Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL
       4             :  * Purpose:  GDALComputedDataset and GDALComputedRasterBand
       5             :  * Author:   Even Rouault <even dot rouault at spatialys.com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com>
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "gdal_priv.h"
      14             : #include "vrtdataset.h"
      15             : 
      16             : #include <cmath>
      17             : #include <limits>
      18             : 
      19             : /************************************************************************/
      20             : /*                        GDALComputedDataset                           */
      21             : /************************************************************************/
      22             : 
      23        1038 : class GDALComputedDataset final : public GDALDataset
      24             : {
      25             :     friend class GDALComputedRasterBand;
      26             : 
      27             :     const GDALComputedRasterBand::Operation m_op;
      28             :     CPLStringList m_aosOptions{};
      29             :     std::vector<std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser>>
      30             :         m_bandDS{};
      31             :     std::vector<GDALRasterBand *> m_poBands{};
      32             :     VRTDataset m_oVRTDS;
      33             : 
      34             :     void AddSources(GDALComputedRasterBand *poBand);
      35             : 
      36             :     static const char *
      37             :     OperationToFunctionName(GDALComputedRasterBand::Operation op);
      38             : 
      39             :     GDALComputedDataset &operator=(const GDALComputedDataset &) = delete;
      40             :     GDALComputedDataset(GDALComputedDataset &&) = delete;
      41             :     GDALComputedDataset &operator=(GDALComputedDataset &&) = delete;
      42             : 
      43             :   public:
      44             :     GDALComputedDataset(const GDALComputedDataset &other);
      45             : 
      46             :     GDALComputedDataset(GDALComputedRasterBand *poBand, int nXSize, int nYSize,
      47             :                         GDALDataType eDT, int nBlockXSize, int nBlockYSize,
      48             :                         GDALComputedRasterBand::Operation op,
      49             :                         const GDALRasterBand *firstBand, double *pFirstConstant,
      50             :                         const GDALRasterBand *secondBand,
      51             :                         double *pSecondConstant);
      52             : 
      53             :     GDALComputedDataset(GDALComputedRasterBand *poBand, int nXSize, int nYSize,
      54             :                         GDALDataType eDT, int nBlockXSize, int nBlockYSize,
      55             :                         GDALComputedRasterBand::Operation op,
      56             :                         const std::vector<const GDALRasterBand *> &bands,
      57             :                         double constant);
      58             : 
      59             :     ~GDALComputedDataset() override;
      60             : 
      61          80 :     CPLErr GetGeoTransform(double *padfGeoTransform) override
      62             :     {
      63          80 :         return m_oVRTDS.GetGeoTransform(padfGeoTransform);
      64             :     }
      65             : 
      66          80 :     const OGRSpatialReference *GetSpatialRef() const override
      67             :     {
      68          80 :         return m_oVRTDS.GetSpatialRef();
      69             :     }
      70             : 
      71           1 :     char **GetMetadata(const char *pszDomain) override
      72             :     {
      73           1 :         return m_oVRTDS.GetMetadata(pszDomain);
      74             :     }
      75             : 
      76         117 :     const char *GetMetadataItem(const char *pszName,
      77             :                                 const char *pszDomain) override
      78             :     {
      79         117 :         return m_oVRTDS.GetMetadataItem(pszName, pszDomain);
      80             :     }
      81             : 
      82          46 :     void *GetInternalHandle(const char *pszHandleName) override
      83             :     {
      84          46 :         if (pszHandleName && EQUAL(pszHandleName, "VRT_DATASET"))
      85          23 :             return &m_oVRTDS;
      86          23 :         return nullptr;
      87             :     }
      88             : };
      89             : 
      90             : /************************************************************************/
      91             : /*                        IsComparisonOperator()                        */
      92             : /************************************************************************/
      93             : 
      94         367 : static bool IsComparisonOperator(GDALComputedRasterBand::Operation op)
      95             : {
      96         367 :     switch (op)
      97             :     {
      98         192 :         case GDALComputedRasterBand::Operation::OP_GT:
      99             :         case GDALComputedRasterBand::Operation::OP_GE:
     100             :         case GDALComputedRasterBand::Operation::OP_LT:
     101             :         case GDALComputedRasterBand::Operation::OP_LE:
     102             :         case GDALComputedRasterBand::Operation::OP_EQ:
     103             :         case GDALComputedRasterBand::Operation::OP_NE:
     104         192 :             return true;
     105         175 :         default:
     106         175 :             break;
     107             :     }
     108         175 :     return false;
     109             : }
     110             : 
     111             : /************************************************************************/
     112             : /*                        GDALComputedDataset()                         */
     113             : /************************************************************************/
     114             : 
     115         279 : GDALComputedDataset::GDALComputedDataset(const GDALComputedDataset &other)
     116         279 :     : GDALDataset(), m_op(other.m_op), m_aosOptions(other.m_aosOptions),
     117         279 :       m_poBands(other.m_poBands),
     118             :       m_oVRTDS(other.GetRasterXSize(), other.GetRasterYSize(),
     119         279 :                other.m_oVRTDS.GetBlockXSize(), other.m_oVRTDS.GetBlockYSize())
     120             : {
     121         279 :     nRasterXSize = other.nRasterXSize;
     122         279 :     nRasterYSize = other.nRasterYSize;
     123             : 
     124             :     auto poBand = new GDALComputedRasterBand(
     125             :         const_cast<const GDALComputedRasterBand &>(
     126         279 :             *cpl::down_cast<GDALComputedRasterBand *>(
     127             :                 const_cast<GDALComputedDataset &>(other).GetRasterBand(1))),
     128         279 :         true);
     129         279 :     SetBand(1, poBand);
     130             : 
     131             :     double adfGT[6];
     132         279 :     if (const_cast<VRTDataset &>(other.m_oVRTDS).GetGeoTransform(adfGT) ==
     133             :         CE_None)
     134             :     {
     135         226 :         m_oVRTDS.SetGeoTransform(adfGT);
     136             :     }
     137             : 
     138         279 :     if (const auto *poSRS =
     139         279 :             const_cast<VRTDataset &>(other.m_oVRTDS).GetSpatialRef())
     140             :     {
     141         226 :         m_oVRTDS.SetSpatialRef(poSRS);
     142             :     }
     143             : 
     144         279 :     m_oVRTDS.AddBand(other.m_oVRTDS.GetRasterBand(1)->GetRasterDataType(),
     145             :                      m_aosOptions.List());
     146             : 
     147         279 :     AddSources(poBand);
     148         279 : }
     149             : 
     150             : /************************************************************************/
     151             : /*                        GDALComputedDataset()                         */
     152             : /************************************************************************/
     153             : 
     154         203 : GDALComputedDataset::GDALComputedDataset(
     155             :     GDALComputedRasterBand *poBand, int nXSize, int nYSize, GDALDataType eDT,
     156             :     int nBlockXSize, int nBlockYSize, GDALComputedRasterBand::Operation op,
     157             :     const GDALRasterBand *firstBand, double *pFirstConstant,
     158         203 :     const GDALRasterBand *secondBand, double *pSecondConstant)
     159         203 :     : m_op(op), m_oVRTDS(nXSize, nYSize, nBlockXSize, nBlockYSize)
     160             : {
     161         203 :     CPLAssert(firstBand != nullptr || secondBand != nullptr);
     162         203 :     if (firstBand)
     163         189 :         m_poBands.push_back(const_cast<GDALRasterBand *>(firstBand));
     164         203 :     if (secondBand)
     165          69 :         m_poBands.push_back(const_cast<GDALRasterBand *>(secondBand));
     166             : 
     167         203 :     nRasterXSize = nXSize;
     168         203 :     nRasterYSize = nYSize;
     169             : 
     170         203 :     if (auto poSrcDS = m_poBands.front()->GetDataset())
     171             :     {
     172             :         double adfGT[6];
     173         203 :         if (poSrcDS->GetGeoTransform(adfGT) == CE_None)
     174             :         {
     175         121 :             m_oVRTDS.SetGeoTransform(adfGT);
     176             :         }
     177             : 
     178         203 :         if (const auto *poSRS = poSrcDS->GetSpatialRef())
     179             :         {
     180         121 :             m_oVRTDS.SetSpatialRef(poSRS);
     181             :         }
     182             :     }
     183             : 
     184         203 :     if (op == GDALComputedRasterBand::Operation::OP_CAST)
     185             :     {
     186             : #ifdef DEBUG
     187             :         // Just for code coverage...
     188          24 :         CPL_IGNORE_RET_VAL(GDALComputedDataset::OperationToFunctionName(op));
     189             : #endif
     190          24 :         m_aosOptions.SetNameValue("subclass", "VRTSourcedRasterBand");
     191             :     }
     192             :     else
     193             :     {
     194         179 :         m_aosOptions.SetNameValue("subclass", "VRTDerivedRasterBand");
     195         179 :         if (IsComparisonOperator(op))
     196             :         {
     197          90 :             m_aosOptions.SetNameValue("PixelFunctionType", "expression");
     198          90 :             if (firstBand && secondBand)
     199             :             {
     200          30 :                 m_aosOptions.SetNameValue(
     201             :                     "_PIXELFN_ARG_expression",
     202             :                     CPLSPrintf(
     203             :                         "source1 %s source2",
     204          30 :                         GDALComputedDataset::OperationToFunctionName(op)));
     205             :             }
     206          60 :             else if (firstBand && pSecondConstant)
     207             :             {
     208          48 :                 m_aosOptions.SetNameValue(
     209             :                     "_PIXELFN_ARG_expression",
     210             :                     CPLSPrintf("source1 %s %.17g",
     211             :                                GDALComputedDataset::OperationToFunctionName(op),
     212          48 :                                *pSecondConstant));
     213             :             }
     214          12 :             else if (pFirstConstant && secondBand)
     215             :             {
     216          12 :                 m_aosOptions.SetNameValue(
     217             :                     "_PIXELFN_ARG_expression",
     218             :                     CPLSPrintf(
     219             :                         "%.17g %s source1", *pFirstConstant,
     220          12 :                         GDALComputedDataset::OperationToFunctionName(op)));
     221             :             }
     222             :             else
     223             :             {
     224           0 :                 CPLAssert(false);
     225             :             }
     226             :         }
     227          89 :         else if (op == GDALComputedRasterBand::Operation::OP_SUBTRACT &&
     228             :                  pSecondConstant)
     229             :         {
     230           2 :             m_aosOptions.SetNameValue("PixelFunctionType", "sum");
     231           2 :             m_aosOptions.SetNameValue("_PIXELFN_ARG_k",
     232           2 :                                       CPLSPrintf("%.17g", -(*pSecondConstant)));
     233             :         }
     234          87 :         else if (op == GDALComputedRasterBand::Operation::OP_DIVIDE)
     235             :         {
     236           6 :             if (pSecondConstant)
     237             :             {
     238           2 :                 m_aosOptions.SetNameValue("PixelFunctionType", "mul");
     239             :                 m_aosOptions.SetNameValue(
     240             :                     "_PIXELFN_ARG_k",
     241           2 :                     CPLSPrintf("%.17g", 1.0 / (*pSecondConstant)));
     242             :             }
     243           4 :             else if (pFirstConstant)
     244             :             {
     245           2 :                 m_aosOptions.SetNameValue("PixelFunctionType", "inv");
     246             :                 m_aosOptions.SetNameValue("_PIXELFN_ARG_k",
     247           2 :                                           CPLSPrintf("%.17g", *pFirstConstant));
     248             :             }
     249             :             else
     250             :             {
     251           2 :                 m_aosOptions.SetNameValue("PixelFunctionType", "div");
     252             :             }
     253             :         }
     254             :         else
     255             :         {
     256             :             m_aosOptions.SetNameValue("PixelFunctionType",
     257          81 :                                       OperationToFunctionName(op));
     258          81 :             if (pSecondConstant)
     259             :                 m_aosOptions.SetNameValue(
     260          58 :                     "_PIXELFN_ARG_k", CPLSPrintf("%.17g", *pSecondConstant));
     261             :         }
     262             :     }
     263         203 :     m_aosOptions.SetNameValue("_PIXELFN_ARG_propagateNoData", "true");
     264         203 :     m_oVRTDS.AddBand(eDT, m_aosOptions.List());
     265             : 
     266         203 :     SetBand(1, poBand);
     267             : 
     268         203 :     AddSources(poBand);
     269         203 : }
     270             : 
     271             : /************************************************************************/
     272             : /*                        GDALComputedDataset()                         */
     273             : /************************************************************************/
     274             : 
     275          37 : GDALComputedDataset::GDALComputedDataset(
     276             :     GDALComputedRasterBand *poBand, int nXSize, int nYSize, GDALDataType eDT,
     277             :     int nBlockXSize, int nBlockYSize, GDALComputedRasterBand::Operation op,
     278          37 :     const std::vector<const GDALRasterBand *> &bands, double constant)
     279          37 :     : m_op(op), m_oVRTDS(nXSize, nYSize, nBlockXSize, nBlockYSize)
     280             : {
     281         132 :     for (const GDALRasterBand *poIterBand : bands)
     282          95 :         m_poBands.push_back(const_cast<GDALRasterBand *>(poIterBand));
     283             : 
     284          37 :     nRasterXSize = nXSize;
     285          37 :     nRasterYSize = nYSize;
     286             : 
     287          37 :     if (auto poSrcDS = m_poBands.front()->GetDataset())
     288             :     {
     289             :         double adfGT[6];
     290          37 :         if (poSrcDS->GetGeoTransform(adfGT) == CE_None)
     291             :         {
     292          16 :             m_oVRTDS.SetGeoTransform(adfGT);
     293             :         }
     294             : 
     295          37 :         if (const auto *poSRS = poSrcDS->GetSpatialRef())
     296             :         {
     297          16 :             m_oVRTDS.SetSpatialRef(poSRS);
     298             :         }
     299             :     }
     300             : 
     301          37 :     m_aosOptions.SetNameValue("subclass", "VRTDerivedRasterBand");
     302          37 :     if (op == GDALComputedRasterBand::Operation::OP_TERNARY)
     303             :     {
     304          18 :         m_aosOptions.SetNameValue("PixelFunctionType", "expression");
     305             :         m_aosOptions.SetNameValue("_PIXELFN_ARG_expression",
     306          18 :                                   "source1 ? source2 : source3");
     307             :     }
     308             :     else
     309             :     {
     310             :         m_aosOptions.SetNameValue("PixelFunctionType",
     311          19 :                                   OperationToFunctionName(op));
     312          19 :         if (!std::isnan(constant))
     313             :         {
     314             :             m_aosOptions.SetNameValue("_PIXELFN_ARG_k",
     315           8 :                                       CPLSPrintf("%.17g", constant));
     316             :         }
     317          19 :         m_aosOptions.SetNameValue("_PIXELFN_ARG_propagateNoData", "true");
     318             :     }
     319          37 :     m_oVRTDS.AddBand(eDT, m_aosOptions.List());
     320             : 
     321          37 :     SetBand(1, poBand);
     322             : 
     323          37 :     AddSources(poBand);
     324          37 : }
     325             : 
     326             : /************************************************************************/
     327             : /*                       ~GDALComputedDataset()                         */
     328             : /************************************************************************/
     329             : 
     330             : GDALComputedDataset::~GDALComputedDataset() = default;
     331             : 
     332             : /************************************************************************/
     333             : /*                       HaveAllBandsSameNoDataValue()                  */
     334             : /************************************************************************/
     335             : 
     336         611 : static bool HaveAllBandsSameNoDataValue(GDALRasterBand **apoBands,
     337             :                                         size_t nBands, bool &hasAtLeastOneNDV,
     338             :                                         double &singleNDV)
     339             : {
     340         611 :     hasAtLeastOneNDV = false;
     341         611 :     singleNDV = 0;
     342             : 
     343         611 :     int bFirstBandHasNoData = false;
     344        1501 :     for (size_t i = 0; i < nBands; ++i)
     345             :     {
     346         896 :         int bHasNoData = false;
     347         896 :         const double dfNoData = apoBands[i]->GetNoDataValue(&bHasNoData);
     348         896 :         if (bHasNoData)
     349          18 :             hasAtLeastOneNDV = true;
     350         896 :         if (i == 0)
     351             :         {
     352         611 :             bFirstBandHasNoData = bHasNoData;
     353         611 :             singleNDV = dfNoData;
     354             :         }
     355         285 :         else if (bHasNoData != bFirstBandHasNoData)
     356             :         {
     357           6 :             return false;
     358             :         }
     359         291 :         else if (bFirstBandHasNoData &&
     360           8 :                  !((std::isnan(singleNDV) && std::isnan(dfNoData)) ||
     361           6 :                    (singleNDV == dfNoData)))
     362             :         {
     363           4 :             return false;
     364             :         }
     365             :     }
     366         605 :     return true;
     367             : }
     368             : 
     369             : /************************************************************************/
     370             : /*                  GDALComputedDataset::AddSources()                   */
     371             : /************************************************************************/
     372             : 
     373         519 : void GDALComputedDataset::AddSources(GDALComputedRasterBand *poBand)
     374             : {
     375             :     auto poSourcedRasterBand =
     376         519 :         cpl::down_cast<VRTSourcedRasterBand *>(m_oVRTDS.GetRasterBand(1));
     377             : 
     378         519 :     bool hasAtLeastOneNDV = false;
     379         519 :     double singleNDV = 0;
     380         519 :     const bool bSameNDV = HaveAllBandsSameNoDataValue(
     381             :         m_poBands.data(), m_poBands.size(), hasAtLeastOneNDV, singleNDV);
     382             : 
     383             :     // For inputs that are instances of GDALComputedDataset, clone them
     384             :     // to make sure we do not depend on temporary instances,
     385             :     // such as "a + b + c", which is evaluated as "(a + b) + c", and the
     386             :     // temporary band/dataset corresponding to a + b will go out of scope
     387             :     // quickly.
     388        1210 :     for (GDALRasterBand *&band : m_poBands)
     389             :     {
     390         691 :         auto poDS = band->GetDataset();
     391         691 :         if (auto poComputedDS = dynamic_cast<GDALComputedDataset *>(poDS))
     392             :         {
     393             :             auto poComputedDSNew =
     394         279 :                 std::make_unique<GDALComputedDataset>(*poComputedDS);
     395         279 :             band = poComputedDSNew->GetRasterBand(1);
     396         279 :             m_bandDS.emplace_back(poComputedDSNew.release());
     397             :         }
     398             : 
     399         691 :         int bHasNoData = false;
     400         691 :         const double dfNoData = band->GetNoDataValue(&bHasNoData);
     401         691 :         if (bHasNoData)
     402             :         {
     403           9 :             poSourcedRasterBand->AddComplexSource(band, -1, -1, -1, -1, -1, -1,
     404             :                                                   -1, -1, 0, 1, dfNoData);
     405             :         }
     406             :         else
     407             :         {
     408         682 :             poSourcedRasterBand->AddSimpleSource(band);
     409             :         }
     410         691 :         poSourcedRasterBand->papoSources[poSourcedRasterBand->nSources - 1]
     411         691 :             ->SetName(CPLSPrintf("source%d", poSourcedRasterBand->nSources));
     412             :     }
     413         519 :     if (hasAtLeastOneNDV)
     414             :     {
     415           5 :         poBand->m_bHasNoData = true;
     416           5 :         if (bSameNDV)
     417             :         {
     418           2 :             poBand->m_dfNoDataValue = singleNDV;
     419             :         }
     420             :         else
     421             :         {
     422           3 :             poBand->m_dfNoDataValue = std::numeric_limits<double>::quiet_NaN();
     423             :         }
     424           5 :         poSourcedRasterBand->SetNoDataValue(poBand->m_dfNoDataValue);
     425             :     }
     426         519 : }
     427             : 
     428             : /************************************************************************/
     429             : /*                       OperationToFunctionName()                      */
     430             : /************************************************************************/
     431             : 
     432         214 : /* static */ const char *GDALComputedDataset::OperationToFunctionName(
     433             :     GDALComputedRasterBand::Operation op)
     434             : {
     435         214 :     const char *ret = "";
     436         214 :     switch (op)
     437             :     {
     438          39 :         case GDALComputedRasterBand::Operation::OP_ADD:
     439          39 :             ret = "sum";
     440          39 :             break;
     441           3 :         case GDALComputedRasterBand::Operation::OP_SUBTRACT:
     442           3 :             ret = "diff";
     443           3 :             break;
     444          38 :         case GDALComputedRasterBand::Operation::OP_MULTIPLY:
     445          38 :             ret = "mul";
     446          38 :             break;
     447           0 :         case GDALComputedRasterBand::Operation::OP_DIVIDE:
     448           0 :             ret = "div";
     449           0 :             break;
     450           9 :         case GDALComputedRasterBand::Operation::OP_MIN:
     451           9 :             ret = "min";
     452           9 :             break;
     453           9 :         case GDALComputedRasterBand::Operation::OP_MAX:
     454           9 :             ret = "max";
     455           9 :             break;
     456           2 :         case GDALComputedRasterBand::Operation::OP_MEAN:
     457           2 :             ret = "mean";
     458           2 :             break;
     459          13 :         case GDALComputedRasterBand::Operation::OP_GT:
     460          13 :             ret = ">";
     461          13 :             break;
     462          16 :         case GDALComputedRasterBand::Operation::OP_GE:
     463          16 :             ret = ">=";
     464          16 :             break;
     465          13 :         case GDALComputedRasterBand::Operation::OP_LT:
     466          13 :             ret = "<";
     467          13 :             break;
     468          14 :         case GDALComputedRasterBand::Operation::OP_LE:
     469          14 :             ret = "<=";
     470          14 :             break;
     471          18 :         case GDALComputedRasterBand::Operation::OP_EQ:
     472          18 :             ret = "==";
     473          18 :             break;
     474          16 :         case GDALComputedRasterBand::Operation::OP_NE:
     475          16 :             ret = "!=";
     476          16 :             break;
     477          24 :         case GDALComputedRasterBand::Operation::OP_CAST:
     478             :         case GDALComputedRasterBand::Operation::OP_TERNARY:
     479          24 :             break;
     480             :     }
     481         214 :     return ret;
     482             : }
     483             : 
     484             : /************************************************************************/
     485             : /*                       GDALComputedRasterBand()                       */
     486             : /************************************************************************/
     487             : 
     488         279 : GDALComputedRasterBand::GDALComputedRasterBand(
     489         279 :     const GDALComputedRasterBand &other, bool)
     490         279 :     : GDALRasterBand()
     491             : {
     492         279 :     nRasterXSize = other.nRasterXSize;
     493         279 :     nRasterYSize = other.nRasterYSize;
     494         279 :     eDataType = other.eDataType;
     495         279 :     nBlockXSize = other.nBlockXSize;
     496         279 :     nBlockYSize = other.nBlockYSize;
     497         279 : }
     498             : 
     499             : //! @cond Doxygen_Suppress
     500             : 
     501             : /************************************************************************/
     502             : /*                       GDALComputedRasterBand()                       */
     503             : /************************************************************************/
     504             : 
     505          37 : GDALComputedRasterBand::GDALComputedRasterBand(
     506             :     Operation op, const std::vector<const GDALRasterBand *> &bands,
     507          37 :     double constant)
     508             : {
     509          37 :     CPLAssert(op == Operation::OP_ADD || op == Operation::OP_MIN ||
     510             :               op == Operation::OP_MAX || op == Operation::OP_MEAN ||
     511             :               op == Operation::OP_TERNARY);
     512             : 
     513          37 :     CPLAssert(!bands.empty());
     514          37 :     nRasterXSize = bands[0]->GetXSize();
     515          37 :     nRasterYSize = bands[0]->GetYSize();
     516          37 :     eDataType = bands[0]->GetRasterDataType();
     517          95 :     for (size_t i = 1; i < bands.size(); ++i)
     518             :     {
     519          58 :         eDataType = GDALDataTypeUnion(eDataType, bands[i]->GetRasterDataType());
     520             :     }
     521             : 
     522          37 :     bool hasAtLeastOneNDV = false;
     523          37 :     double singleNDV = 0;
     524             :     const bool bSameNDV =
     525          37 :         HaveAllBandsSameNoDataValue(const_cast<GDALRasterBand **>(bands.data()),
     526             :                                     bands.size(), hasAtLeastOneNDV, singleNDV);
     527             : 
     528          37 :     if (!bSameNDV)
     529             :     {
     530           0 :         eDataType = eDataType == GDT_Float64 ? GDT_Float64 : GDT_Float32;
     531             :     }
     532          37 :     else if (op == Operation::OP_TERNARY)
     533             :     {
     534          18 :         CPLAssert(bands.size() == 3);
     535          18 :         eDataType = GDALDataTypeUnion(bands[1]->GetRasterDataType(),
     536          18 :                                       bands[2]->GetRasterDataType());
     537             :     }
     538          19 :     else if (!std::isnan(constant) && eDataType != GDT_Float64)
     539             :     {
     540           4 :         if (op == Operation::OP_MIN || op == Operation::OP_MAX)
     541             :         {
     542           4 :             eDataType = GDALDataTypeUnionWithValue(eDataType, constant, false);
     543             :         }
     544             :         else
     545             :         {
     546           0 :             eDataType = (static_cast<float>(constant) == constant)
     547           0 :                             ? GDT_Float32
     548             :                             : GDT_Float64;
     549             :         }
     550             :     }
     551          37 :     bands[0]->GetBlockSize(&nBlockXSize, &nBlockYSize);
     552             :     auto l_poDS = std::make_unique<GDALComputedDataset>(
     553          37 :         this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
     554          74 :         op, bands, constant);
     555          37 :     m_poOwningDS.reset(l_poDS.release());
     556          37 : }
     557             : 
     558             : /************************************************************************/
     559             : /*                       GDALComputedRasterBand()                       */
     560             : /************************************************************************/
     561             : 
     562          55 : GDALComputedRasterBand::GDALComputedRasterBand(Operation op,
     563             :                                                const GDALRasterBand &firstBand,
     564          55 :                                                const GDALRasterBand &secondBand)
     565             : {
     566          55 :     nRasterXSize = firstBand.GetXSize();
     567          55 :     nRasterYSize = firstBand.GetYSize();
     568             : 
     569          55 :     bool hasAtLeastOneNDV = false;
     570          55 :     double singleNDV = 0;
     571             :     GDALRasterBand *apoBands[] = {const_cast<GDALRasterBand *>(&firstBand),
     572          55 :                                   const_cast<GDALRasterBand *>(&secondBand)};
     573             :     const bool bSameNDV =
     574          55 :         HaveAllBandsSameNoDataValue(apoBands, 2, hasAtLeastOneNDV, singleNDV);
     575             : 
     576          55 :     const auto firstDT = firstBand.GetRasterDataType();
     577          55 :     const auto secondDT = secondBand.GetRasterDataType();
     578          55 :     if (!bSameNDV)
     579           6 :         eDataType = (firstDT == GDT_Float64 || secondDT == GDT_Float64)
     580           6 :                         ? GDT_Float64
     581             :                         : GDT_Float32;
     582          52 :     else if (IsComparisonOperator(op))
     583          30 :         eDataType = GDT_Byte;
     584          22 :     else if (op == Operation::OP_ADD && firstDT == GDT_Byte &&
     585             :              secondDT == GDT_Byte)
     586           2 :         eDataType = GDT_UInt16;
     587          20 :     else if (firstDT == GDT_Float32 && secondDT == GDT_Float32)
     588           1 :         eDataType = GDT_Float32;
     589          19 :     else if ((op == Operation::OP_MIN || op == Operation::OP_MAX) &&
     590             :              firstDT == secondDT)
     591           1 :         eDataType = firstDT;
     592             :     else
     593          18 :         eDataType = GDT_Float64;
     594          55 :     firstBand.GetBlockSize(&nBlockXSize, &nBlockYSize);
     595             :     auto l_poDS = std::make_unique<GDALComputedDataset>(
     596          55 :         this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
     597         110 :         op, &firstBand, nullptr, &secondBand, nullptr);
     598          55 :     m_poOwningDS.reset(l_poDS.release());
     599          55 : }
     600             : 
     601             : /************************************************************************/
     602             : /*                       GDALComputedRasterBand()                       */
     603             : /************************************************************************/
     604             : 
     605          14 : GDALComputedRasterBand::GDALComputedRasterBand(Operation op, double constant,
     606          14 :                                                const GDALRasterBand &band)
     607             : {
     608          14 :     CPLAssert(op == Operation::OP_DIVIDE || IsComparisonOperator(op));
     609             : 
     610          14 :     nRasterXSize = band.GetXSize();
     611          14 :     nRasterYSize = band.GetYSize();
     612          14 :     const auto firstDT = band.GetRasterDataType();
     613          14 :     if (IsComparisonOperator(op))
     614          12 :         eDataType = GDT_Byte;
     615           2 :     else if (firstDT == GDT_Float32 && static_cast<float>(constant) == constant)
     616           0 :         eDataType = GDT_Float32;
     617             :     else
     618           2 :         eDataType = GDT_Float64;
     619          14 :     band.GetBlockSize(&nBlockXSize, &nBlockYSize);
     620             :     auto l_poDS = std::make_unique<GDALComputedDataset>(
     621          14 :         this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
     622          28 :         op, nullptr, &constant, &band, nullptr);
     623          14 :     m_poOwningDS.reset(l_poDS.release());
     624          14 : }
     625             : 
     626             : /************************************************************************/
     627             : /*                       GDALComputedRasterBand()                       */
     628             : /************************************************************************/
     629             : 
     630         110 : GDALComputedRasterBand::GDALComputedRasterBand(Operation op,
     631             :                                                const GDALRasterBand &band,
     632         110 :                                                double constant)
     633             : {
     634         110 :     nRasterXSize = band.GetXSize();
     635         110 :     nRasterYSize = band.GetYSize();
     636         110 :     const auto firstDT = band.GetRasterDataType();
     637         110 :     if (IsComparisonOperator(op))
     638          48 :         eDataType = GDT_Byte;
     639          62 :     else if (op == Operation::OP_ADD && firstDT == GDT_Byte &&
     640          10 :              constant >= -128 && constant <= 127 &&
     641          10 :              std::floor(constant) == constant)
     642          10 :         eDataType = GDT_Byte;
     643          52 :     else if (firstDT == GDT_Float32 && static_cast<float>(constant) == constant)
     644           8 :         eDataType = GDT_Float32;
     645             :     else
     646          44 :         eDataType = GDT_Float64;
     647         110 :     band.GetBlockSize(&nBlockXSize, &nBlockYSize);
     648             :     auto l_poDS = std::make_unique<GDALComputedDataset>(
     649         110 :         this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
     650         220 :         op, &band, nullptr, nullptr, &constant);
     651         110 :     m_poOwningDS.reset(l_poDS.release());
     652         110 : }
     653             : 
     654             : /************************************************************************/
     655             : /*                       GDALComputedRasterBand()                       */
     656             : /************************************************************************/
     657             : 
     658          24 : GDALComputedRasterBand::GDALComputedRasterBand(Operation op,
     659             :                                                const GDALRasterBand &band,
     660          24 :                                                GDALDataType dt)
     661             : {
     662          24 :     CPLAssert(op == Operation::OP_CAST);
     663          24 :     nRasterXSize = band.GetXSize();
     664          24 :     nRasterYSize = band.GetYSize();
     665          24 :     eDataType = dt;
     666          24 :     band.GetBlockSize(&nBlockXSize, &nBlockYSize);
     667             :     auto l_poDS = std::make_unique<GDALComputedDataset>(
     668          24 :         this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
     669          48 :         op, &band, nullptr, nullptr, nullptr);
     670          24 :     m_poOwningDS.reset(l_poDS.release());
     671          24 : }
     672             : 
     673             : //! @endcond
     674             : 
     675             : /************************************************************************/
     676             : /*                      ~GDALComputedRasterBand()                       */
     677             : /************************************************************************/
     678             : 
     679         929 : GDALComputedRasterBand::~GDALComputedRasterBand()
     680             : {
     681         519 :     if (m_poOwningDS)
     682         240 :         cpl::down_cast<GDALComputedDataset *>(m_poOwningDS.get())->nBands = 0;
     683         519 :     poDS = nullptr;
     684         929 : }
     685             : 
     686             : /************************************************************************/
     687             : /*                  GDALComputedRasterBand::GetNoDataValue()            */
     688             : /************************************************************************/
     689             : 
     690         849 : double GDALComputedRasterBand::GetNoDataValue(int *pbHasNoData)
     691             : {
     692         849 :     if (pbHasNoData)
     693         849 :         *pbHasNoData = m_bHasNoData;
     694         849 :     return m_dfNoDataValue;
     695             : }
     696             : 
     697             : /************************************************************************/
     698             : /*                    GDALComputedRasterBandRelease()                   */
     699             : /************************************************************************/
     700             : 
     701             : /** Release a GDALComputedRasterBandH
     702             :  *
     703             :  * @since 3.12
     704             :  */
     705         131 : void GDALComputedRasterBandRelease(GDALComputedRasterBandH hBand)
     706             : {
     707         131 :     delete GDALComputedRasterBand::FromHandle(hBand);
     708         131 : }
     709             : 
     710             : /************************************************************************/
     711             : /*                           IReadBlock()                               */
     712             : /************************************************************************/
     713             : 
     714         168 : CPLErr GDALComputedRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
     715             :                                           void *pData)
     716             : {
     717         168 :     auto l_poDS = cpl::down_cast<GDALComputedDataset *>(poDS);
     718         168 :     return l_poDS->m_oVRTDS.GetRasterBand(1)->ReadBlock(nBlockXOff, nBlockYOff,
     719         168 :                                                         pData);
     720             : }
     721             : 
     722             : /************************************************************************/
     723             : /*                           IRasterIO()                                */
     724             : /************************************************************************/
     725             : 
     726         151 : CPLErr GDALComputedRasterBand::IRasterIO(
     727             :     GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
     728             :     void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
     729             :     GSpacing nPixelSpace, GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg)
     730             : {
     731         151 :     auto l_poDS = cpl::down_cast<GDALComputedDataset *>(poDS);
     732         151 :     return l_poDS->m_oVRTDS.GetRasterBand(1)->RasterIO(
     733             :         eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
     734         151 :         eBufType, nPixelSpace, nLineSpace, psExtraArg);
     735             : }