LCOV - code coverage report
Current view: top level - frmts/vrt - vrtprocesseddataset.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 691 842 82.1 %
Date: 2025-01-18 12:42:00 Functions: 23 24 95.8 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  Virtual GDAL Datasets
       4             :  * Purpose:  Implementation of VRTProcessedDataset.
       5             :  * Author:   Even Rouault <even.rouault at spatialys.com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2024, Even Rouault <even.rouault at spatialys.com>
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "cpl_minixml.h"
      14             : #include "cpl_string.h"
      15             : #include "gdal_utils.h"
      16             : #include "vrtdataset.h"
      17             : 
      18             : #include <algorithm>
      19             : #include <limits>
      20             : #include <map>
      21             : #include <vector>
      22             : 
      23             : /************************************************************************/
      24             : /*                        VRTProcessedDatasetFunc                       */
      25             : /************************************************************************/
      26             : 
      27             : //! Structure holding information for a VRTProcessedDataset function.
      28             : struct VRTProcessedDatasetFunc
      29             : {
      30             :     //! Processing function name
      31             :     std::string osFuncName{};
      32             : 
      33             :     //! User data to provide to pfnInit, pfnFree, pfnProcess callbacks.
      34             :     void *pUserData = nullptr;
      35             : 
      36             :     //! Whether XML metadata has been specified
      37             :     bool bMetadataSpecified = false;
      38             : 
      39             :     //! Map of (constant argument name, constant value)
      40             :     std::map<std::string, std::string> oMapConstantArguments{};
      41             : 
      42             :     //! Set of builtin argument names (e.g "offset", "scale", "nodata")
      43             :     std::set<std::string> oSetBuiltinArguments{};
      44             : 
      45             :     //! Arguments defined in the VRT
      46             :     struct OtherArgument
      47             :     {
      48             :         std::string osType{};
      49             :         bool bRequired = false;
      50             :     };
      51             : 
      52             :     std::map<std::string, OtherArgument> oOtherArguments{};
      53             : 
      54             :     //! Requested input data type.
      55             :     GDALDataType eRequestedInputDT = GDT_Unknown;
      56             : 
      57             :     //! List of supported input datatypes. Empty if no restriction.
      58             :     std::vector<GDALDataType> aeSupportedInputDT{};
      59             : 
      60             :     //! List of supported input band counts. Empty if no restriction.
      61             :     std::vector<int> anSupportedInputBandCount{};
      62             : 
      63             :     //! Optional initialization function
      64             :     GDALVRTProcessedDatasetFuncInit pfnInit = nullptr;
      65             : 
      66             :     //! Optional free function
      67             :     GDALVRTProcessedDatasetFuncFree pfnFree = nullptr;
      68             : 
      69             :     //! Required processing function
      70             :     GDALVRTProcessedDatasetFuncProcess pfnProcess = nullptr;
      71             : };
      72             : 
      73             : /************************************************************************/
      74             : /*                      GetGlobalMapProcessedDatasetFunc()              */
      75             : /************************************************************************/
      76             : 
      77             : /** Return the registry of VRTProcessedDatasetFunc functions */
      78             : static std::map<std::string, VRTProcessedDatasetFunc> &
      79        7082 : GetGlobalMapProcessedDatasetFunc()
      80             : {
      81        7082 :     static std::map<std::string, VRTProcessedDatasetFunc> goMap;
      82        7082 :     return goMap;
      83             : }
      84             : 
      85             : /************************************************************************/
      86             : /*                            Step::~Step()                             */
      87             : /************************************************************************/
      88             : 
      89             : /*! @cond Doxygen_Suppress */
      90             : 
      91             : /** Step destructor */
      92         128 : VRTProcessedDataset::Step::~Step()
      93             : {
      94         128 :     deinit();
      95         128 : }
      96             : 
      97             : /************************************************************************/
      98             : /*                           Step::deinit()                             */
      99             : /************************************************************************/
     100             : 
     101             : /** Free pWorkingData */
     102         128 : void VRTProcessedDataset::Step::deinit()
     103             : {
     104         128 :     if (pWorkingData)
     105             :     {
     106          48 :         const auto &oMapFunctions = GetGlobalMapProcessedDatasetFunc();
     107          48 :         const auto oIterFunc = oMapFunctions.find(osAlgorithm);
     108          48 :         if (oIterFunc != oMapFunctions.end())
     109             :         {
     110          48 :             if (oIterFunc->second.pfnFree)
     111             :             {
     112          96 :                 oIterFunc->second.pfnFree(osAlgorithm.c_str(),
     113          48 :                                           oIterFunc->second.pUserData,
     114             :                                           pWorkingData);
     115             :             }
     116             :         }
     117             :         else
     118             :         {
     119           0 :             CPLAssert(false);
     120             :         }
     121          48 :         pWorkingData = nullptr;
     122             :     }
     123         128 : }
     124             : 
     125             : /************************************************************************/
     126             : /*                        Step::Step(Step&& other)                      */
     127             : /************************************************************************/
     128             : 
     129             : /** Move constructor */
     130          51 : VRTProcessedDataset::Step::Step(Step &&other)
     131          51 :     : osAlgorithm(std::move(other.osAlgorithm)),
     132         102 :       aosArguments(std::move(other.aosArguments)), eInDT(other.eInDT),
     133          51 :       eOutDT(other.eOutDT), nInBands(other.nInBands),
     134          51 :       nOutBands(other.nOutBands), adfInNoData(other.adfInNoData),
     135          51 :       adfOutNoData(other.adfOutNoData), pWorkingData(other.pWorkingData)
     136             : {
     137          51 :     other.pWorkingData = nullptr;
     138          51 : }
     139             : 
     140             : /************************************************************************/
     141             : /*                      Step operator=(Step&& other)                    */
     142             : /************************************************************************/
     143             : 
     144             : /** Move assignment operator */
     145           0 : VRTProcessedDataset::Step &VRTProcessedDataset::Step::operator=(Step &&other)
     146             : {
     147           0 :     if (&other != this)
     148             :     {
     149           0 :         deinit();
     150           0 :         osAlgorithm = std::move(other.osAlgorithm);
     151           0 :         aosArguments = std::move(other.aosArguments);
     152           0 :         eInDT = other.eInDT;
     153           0 :         eOutDT = other.eOutDT;
     154           0 :         nInBands = other.nInBands;
     155           0 :         nOutBands = other.nOutBands;
     156           0 :         adfInNoData = std::move(other.adfInNoData);
     157           0 :         adfOutNoData = std::move(other.adfOutNoData);
     158           0 :         std::swap(pWorkingData, other.pWorkingData);
     159             :     }
     160           0 :     return *this;
     161             : }
     162             : 
     163             : /************************************************************************/
     164             : /*                        VRTProcessedDataset()                         */
     165             : /************************************************************************/
     166             : 
     167             : /** Constructor */
     168          87 : VRTProcessedDataset::VRTProcessedDataset(int nXSize, int nYSize)
     169          87 :     : VRTDataset(nXSize, nYSize)
     170             : {
     171          87 : }
     172             : 
     173             : /************************************************************************/
     174             : /*                       ~VRTProcessedDataset()                         */
     175             : /************************************************************************/
     176             : 
     177         174 : VRTProcessedDataset::~VRTProcessedDataset()
     178             : 
     179             : {
     180          87 :     VRTProcessedDataset::FlushCache(true);
     181          87 :     VRTProcessedDataset::CloseDependentDatasets();
     182         174 : }
     183             : 
     184             : /************************************************************************/
     185             : /*                              XMLInit()                               */
     186             : /************************************************************************/
     187             : 
     188             : /** Instantiate object from XML tree */
     189          84 : CPLErr VRTProcessedDataset::XMLInit(const CPLXMLNode *psTree,
     190             :                                     const char *pszVRTPathIn)
     191             : 
     192             : {
     193          84 :     if (Init(psTree, pszVRTPathIn, nullptr, nullptr, -1) != CE_None)
     194          41 :         return CE_Failure;
     195             : 
     196          43 :     const auto poSrcFirstBand = m_poSrcDS->GetRasterBand(1);
     197          43 :     const int nOvrCount = poSrcFirstBand->GetOverviewCount();
     198          46 :     for (int i = 0; i < nOvrCount; ++i)
     199             :     {
     200           3 :         auto poOvrDS = std::make_unique<VRTProcessedDataset>(0, 0);
     201           3 :         if (poOvrDS->Init(psTree, pszVRTPathIn, this, m_poSrcDS.get(), i) !=
     202             :             CE_None)
     203           0 :             break;
     204           3 :         m_apoOverviewDatasets.emplace_back(std::move(poOvrDS));
     205             :     }
     206             : 
     207          43 :     return CE_None;
     208             : }
     209             : 
     210          72 : static bool HasScaleOffset(GDALDataset &oSrcDS)
     211             : {
     212         281 :     for (int i = 1; i <= oSrcDS.GetRasterCount(); i++)
     213             :     {
     214             :         int pbSuccess;
     215         211 :         GDALRasterBand &oBand = *oSrcDS.GetRasterBand(i);
     216         211 :         double scale = oBand.GetScale(&pbSuccess);
     217         211 :         if (pbSuccess && scale != 1)
     218             :         {
     219           2 :             return true;
     220             :         }
     221         209 :         double offset = oBand.GetOffset(&pbSuccess);
     222         209 :         if (pbSuccess && offset != 0)
     223             :         {
     224           0 :             return true;
     225             :         }
     226             :     }
     227             : 
     228          70 :     return false;
     229             : }
     230             : 
     231             : /** Instantiate object from XML tree */
     232          87 : CPLErr VRTProcessedDataset::Init(const CPLXMLNode *psTree,
     233             :                                  const char *pszVRTPathIn,
     234             :                                  const VRTProcessedDataset *poParentDS,
     235             :                                  GDALDataset *poParentSrcDS, int iOvrLevel)
     236             : 
     237             : {
     238          87 :     const CPLXMLNode *psInput = CPLGetXMLNode(psTree, "Input");
     239          87 :     if (!psInput)
     240             :     {
     241           1 :         CPLError(CE_Failure, CPLE_AppDefined, "Input element missing");
     242           1 :         return CE_Failure;
     243             :     }
     244             : 
     245          86 :     if (pszVRTPathIn)
     246          12 :         m_osVRTPath = pszVRTPathIn;
     247             : 
     248          86 :     if (poParentSrcDS)
     249             :     {
     250           3 :         m_poSrcDS.reset(
     251             :             GDALCreateOverviewDataset(poParentSrcDS, iOvrLevel, true));
     252             :     }
     253          83 :     else if (const CPLXMLNode *psSourceFileNameNode =
     254          83 :                  CPLGetXMLNode(psInput, "SourceFilename"))
     255             :     {
     256          81 :         const bool bRelativeToVRT = CPL_TO_BOOL(
     257             :             atoi(CPLGetXMLValue(psSourceFileNameNode, "relativetoVRT", "0")));
     258             :         const std::string osFilename = VRTDataset::BuildSourceFilename(
     259             :             CPLGetXMLValue(psInput, "SourceFilename", ""), pszVRTPathIn,
     260         162 :             bRelativeToVRT);
     261          81 :         m_poSrcDS.reset(GDALDataset::Open(
     262             :             osFilename.c_str(), GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, nullptr,
     263             :             nullptr, nullptr));
     264             :     }
     265           2 :     else if (const CPLXMLNode *psVRTDataset =
     266           2 :                  CPLGetXMLNode(psInput, "VRTDataset"))
     267             :     {
     268           1 :         CPLXMLNode sVRTDatasetTmp = *psVRTDataset;
     269           1 :         sVRTDatasetTmp.psNext = nullptr;
     270           1 :         char *pszXML = CPLSerializeXMLTree(&sVRTDatasetTmp);
     271           1 :         m_poSrcDS.reset(VRTDataset::OpenXML(pszXML, pszVRTPathIn, GA_ReadOnly));
     272           1 :         CPLFree(pszXML);
     273             :     }
     274             :     else
     275             :     {
     276           1 :         CPLError(
     277             :             CE_Failure, CPLE_AppDefined,
     278             :             "Input element should have a SourceFilename or VRTDataset element");
     279           1 :         return CE_Failure;
     280             :     }
     281             : 
     282          85 :     if (!m_poSrcDS)
     283           2 :         return CE_Failure;
     284             : 
     285          83 :     const char *pszUnscale = CPLGetXMLValue(psInput, "unscale", "AUTO");
     286          83 :     bool bUnscale = false;
     287          83 :     if (EQUAL(pszUnscale, "AUTO"))
     288             :     {
     289          72 :         if (HasScaleOffset(*m_poSrcDS))
     290             :         {
     291           2 :             bUnscale = true;
     292             :         }
     293             :     }
     294          11 :     else if (EQUAL(pszUnscale, "YES") || EQUAL(pszUnscale, "ON") ||
     295          11 :              EQUAL(pszUnscale, "TRUE") || EQUAL(pszUnscale, "1"))
     296             :     {
     297           6 :         bUnscale = true;
     298             :     }
     299           5 :     else if (!(EQUAL(pszUnscale, "NO") || EQUAL(pszUnscale, "OFF") ||
     300           5 :                EQUAL(pszUnscale, "FALSE") || EQUAL(pszUnscale, "0")))
     301             :     {
     302           1 :         CPLError(CE_Failure, CPLE_AppDefined, "Invalid value of 'unscale'");
     303           1 :         return CE_Failure;
     304             :     }
     305             : 
     306          82 :     if (bUnscale)
     307             :     {
     308           8 :         CPLStringList oArgs;
     309           8 :         oArgs.AddString("-unscale");
     310           8 :         oArgs.AddString("-ot");
     311           8 :         oArgs.AddString("Float64");
     312           8 :         oArgs.AddString("-of");
     313           8 :         oArgs.AddString("VRT");
     314           8 :         oArgs.AddString("-a_nodata");
     315           8 :         oArgs.AddString("nan");
     316           8 :         auto *poArgs = GDALTranslateOptionsNew(oArgs.List(), nullptr);
     317             :         int pbUsageError;
     318           8 :         CPLAssert(poArgs);
     319           8 :         m_poVRTSrcDS.reset(m_poSrcDS.release());
     320             :         // https://trac.cppcheck.net/ticket/11325
     321             :         // cppcheck-suppress accessMoved
     322           8 :         m_poSrcDS.reset(GDALDataset::FromHandle(
     323           8 :             GDALTranslate("", m_poVRTSrcDS.get(), poArgs, &pbUsageError)));
     324           8 :         GDALTranslateOptionsFree(poArgs);
     325             : 
     326           8 :         if (pbUsageError || !m_poSrcDS)
     327             :         {
     328           0 :             return CE_Failure;
     329             :         }
     330             :     }
     331             : 
     332          82 :     if (nRasterXSize == 0 && nRasterYSize == 0)
     333             :     {
     334          80 :         nRasterXSize = m_poSrcDS->GetRasterXSize();
     335          80 :         nRasterYSize = m_poSrcDS->GetRasterYSize();
     336             :     }
     337           2 :     else if (nRasterXSize != m_poSrcDS->GetRasterXSize() ||
     338           0 :              nRasterYSize != m_poSrcDS->GetRasterYSize())
     339             :     {
     340           2 :         CPLError(CE_Failure, CPLE_AppDefined,
     341             :                  "Inconsistent declared VRT dimensions with input dataset");
     342           2 :         return CE_Failure;
     343             :     }
     344             : 
     345          80 :     if (m_poSrcDS->GetRasterCount() == 0)
     346           0 :         return CE_Failure;
     347             : 
     348             :     // Inherit SRS from source if not explicitly defined in VRT
     349          80 :     if (!CPLGetXMLNode(psTree, "SRS"))
     350             :     {
     351          80 :         const OGRSpatialReference *poSRS = m_poSrcDS->GetSpatialRef();
     352          80 :         if (poSRS)
     353             :         {
     354          10 :             m_poSRS.reset(poSRS->Clone());
     355             :         }
     356             :     }
     357             : 
     358             :     // Inherit GeoTransform from source if not explicitly defined in VRT
     359          80 :     if (iOvrLevel < 0 && !CPLGetXMLNode(psTree, "GeoTransform"))
     360             :     {
     361          77 :         if (m_poSrcDS->GetGeoTransform(m_adfGeoTransform) == CE_None)
     362          40 :             m_bGeoTransformSet = true;
     363             :     }
     364             : 
     365             :     /* -------------------------------------------------------------------- */
     366             :     /*      Initialize blocksize before calling sub-init so that the        */
     367             :     /*      band initializers can get it from the dataset object when       */
     368             :     /*      they are created.                                               */
     369             :     /* -------------------------------------------------------------------- */
     370             : 
     371          80 :     const auto poSrcFirstBand = m_poSrcDS->GetRasterBand(1);
     372          80 :     poSrcFirstBand->GetBlockSize(&m_nBlockXSize, &m_nBlockYSize);
     373          80 :     bool bUserBlockSize = false;
     374          80 :     if (const char *pszBlockXSize =
     375          80 :             CPLGetXMLValue(psTree, "BlockXSize", nullptr))
     376             :     {
     377           0 :         bUserBlockSize = true;
     378           0 :         m_nBlockXSize = atoi(pszBlockXSize);
     379           0 :         if (m_nBlockXSize <= 1)
     380             :         {
     381           0 :             CPLError(CE_Failure, CPLE_AppDefined, "Invalid BlockXSize");
     382           0 :             return CE_Failure;
     383             :         }
     384             :     }
     385          80 :     if (const char *pszBlockYSize =
     386          80 :             CPLGetXMLValue(psTree, "BlockYSize", nullptr))
     387             :     {
     388           0 :         bUserBlockSize = true;
     389           0 :         m_nBlockYSize = atoi(pszBlockYSize);
     390           0 :         if (m_nBlockYSize <= 1)
     391             :         {
     392           0 :             CPLError(CE_Failure, CPLE_AppDefined, "Invalid BlockYSize");
     393           0 :             return CE_Failure;
     394             :         }
     395             :     }
     396             : 
     397             :     // Initialize all the general VRT stuff.
     398          80 :     if (VRTDataset::XMLInit(psTree, pszVRTPathIn) != CE_None)
     399             :     {
     400           0 :         return CE_Failure;
     401             :     }
     402             : 
     403             :     // Use geotransform from parent for overviews
     404          80 :     if (iOvrLevel >= 0 && poParentDS->m_bGeoTransformSet)
     405             :     {
     406           1 :         m_bGeoTransformSet = true;
     407           1 :         m_adfGeoTransform[0] = poParentDS->m_adfGeoTransform[0];
     408           1 :         m_adfGeoTransform[1] = poParentDS->m_adfGeoTransform[1];
     409           1 :         m_adfGeoTransform[2] = poParentDS->m_adfGeoTransform[2];
     410           1 :         m_adfGeoTransform[3] = poParentDS->m_adfGeoTransform[3];
     411           1 :         m_adfGeoTransform[4] = poParentDS->m_adfGeoTransform[4];
     412           1 :         m_adfGeoTransform[5] = poParentDS->m_adfGeoTransform[5];
     413             : 
     414           1 :         m_adfGeoTransform[1] *=
     415           1 :             static_cast<double>(poParentDS->GetRasterXSize()) / nRasterXSize;
     416           1 :         m_adfGeoTransform[2] *=
     417           1 :             static_cast<double>(poParentDS->GetRasterYSize()) / nRasterYSize;
     418           1 :         m_adfGeoTransform[4] *=
     419           1 :             static_cast<double>(poParentDS->GetRasterXSize()) / nRasterXSize;
     420           1 :         m_adfGeoTransform[5] *=
     421           1 :             static_cast<double>(poParentDS->GetRasterYSize()) / nRasterYSize;
     422             :     }
     423             : 
     424          80 :     const CPLXMLNode *psOutputBands = CPLGetXMLNode(psTree, "OutputBands");
     425          80 :     if (psOutputBands)
     426             :     {
     427          14 :         if (const char *pszCount =
     428          14 :                 CPLGetXMLValue(psOutputBands, "count", nullptr))
     429             :         {
     430          14 :             if (EQUAL(pszCount, "FROM_LAST_STEP"))
     431             :             {
     432           8 :                 m_outputBandCountProvenance = ValueProvenance::FROM_LAST_STEP;
     433             :             }
     434           6 :             else if (!EQUAL(pszCount, "FROM_SOURCE"))
     435             :             {
     436           4 :                 if (CPLGetValueType(pszCount) == CPL_VALUE_INTEGER)
     437             :                 {
     438           3 :                     m_outputBandCountProvenance =
     439             :                         ValueProvenance::USER_PROVIDED;
     440           3 :                     m_outputBandCountValue = atoi(pszCount);
     441           3 :                     if (!GDALCheckBandCount(m_outputBandCountValue,
     442             :                                             /* bIsZeroAllowed = */ false))
     443             :                     {
     444             :                         // Error emitted by above function
     445           1 :                         return CE_Failure;
     446             :                     }
     447             :                 }
     448             :                 else
     449             :                 {
     450           1 :                     CPLError(CE_Failure, CPLE_AppDefined,
     451             :                              "Invalid value for OutputBands.count");
     452           1 :                     return CE_Failure;
     453             :                 }
     454             :             }
     455             :         }
     456             : 
     457          12 :         if (const char *pszType =
     458          12 :                 CPLGetXMLValue(psOutputBands, "dataType", nullptr))
     459             :         {
     460          12 :             if (EQUAL(pszType, "FROM_LAST_STEP"))
     461             :             {
     462           4 :                 m_outputBandDataTypeProvenance =
     463             :                     ValueProvenance::FROM_LAST_STEP;
     464             :             }
     465           8 :             else if (!EQUAL(pszType, "FROM_SOURCE"))
     466             :             {
     467           6 :                 m_outputBandDataTypeProvenance = ValueProvenance::USER_PROVIDED;
     468           6 :                 m_outputBandDataTypeValue = GDALGetDataTypeByName(pszType);
     469           6 :                 if (m_outputBandDataTypeValue == GDT_Unknown)
     470             :                 {
     471           1 :                     CPLError(CE_Failure, CPLE_AppDefined,
     472             :                              "Invalid value for OutputBands.dataType");
     473           1 :                     return CE_Failure;
     474             :                 }
     475             :             }
     476             :         }
     477             :     }
     478          66 :     else if (CPLGetXMLNode(psTree, "VRTRasterBand"))
     479             :     {
     480           5 :         m_outputBandCountProvenance = ValueProvenance::FROM_VRTRASTERBAND;
     481           5 :         m_outputBandDataTypeProvenance = ValueProvenance::FROM_VRTRASTERBAND;
     482             :     }
     483             : 
     484          77 :     int nOutputBandCount = 0;
     485          77 :     switch (m_outputBandCountProvenance)
     486             :     {
     487           1 :         case ValueProvenance::USER_PROVIDED:
     488           1 :             nOutputBandCount = m_outputBandCountValue;
     489           1 :             break;
     490          63 :         case ValueProvenance::FROM_SOURCE:
     491          63 :             nOutputBandCount = m_poSrcDS->GetRasterCount();
     492          63 :             break;
     493           5 :         case ValueProvenance::FROM_VRTRASTERBAND:
     494           5 :             nOutputBandCount = nBands;
     495           5 :             break;
     496           8 :         case ValueProvenance::FROM_LAST_STEP:
     497           8 :             break;
     498             :     }
     499             : 
     500             :     const CPLXMLNode *psProcessingSteps =
     501          77 :         CPLGetXMLNode(psTree, "ProcessingSteps");
     502          77 :     if (!psProcessingSteps)
     503             :     {
     504           1 :         CPLError(CE_Failure, CPLE_AppDefined,
     505             :                  "ProcessingSteps element missing");
     506           1 :         return CE_Failure;
     507             :     }
     508             : 
     509          76 :     const auto eInDT = poSrcFirstBand->GetRasterDataType();
     510         219 :     for (int i = 1; i < m_poSrcDS->GetRasterCount(); ++i)
     511             :     {
     512         143 :         const auto eDT = m_poSrcDS->GetRasterBand(i + 1)->GetRasterDataType();
     513         143 :         if (eDT != eInDT)
     514             :         {
     515           0 :             CPLError(CE_Warning, CPLE_AppDefined,
     516             :                      "Not all bands of the input dataset have the same data "
     517             :                      "type. The data type of the first band will be used as "
     518             :                      "the reference one.");
     519           0 :             break;
     520             :         }
     521             :     }
     522             : 
     523          76 :     GDALDataType eCurrentDT = eInDT;
     524          76 :     int nCurrentBandCount = m_poSrcDS->GetRasterCount();
     525             : 
     526         152 :     std::vector<double> adfNoData;
     527         295 :     for (int i = 1; i <= nCurrentBandCount; ++i)
     528             :     {
     529         219 :         int bHasVal = FALSE;
     530             :         const double dfVal =
     531         219 :             m_poSrcDS->GetRasterBand(i)->GetNoDataValue(&bHasVal);
     532             :         adfNoData.emplace_back(
     533         219 :             bHasVal ? dfVal : std::numeric_limits<double>::quiet_NaN());
     534             :     }
     535             : 
     536          76 :     int nStepCount = 0;
     537         153 :     for (const CPLXMLNode *psStep = psProcessingSteps->psChild; psStep;
     538          77 :          psStep = psStep->psNext)
     539             :     {
     540          77 :         if (psStep->eType == CXT_Element &&
     541          77 :             strcmp(psStep->pszValue, "Step") == 0)
     542             :         {
     543          77 :             ++nStepCount;
     544             :         }
     545             :     }
     546             : 
     547          76 :     int iStep = 0;
     548         124 :     for (const CPLXMLNode *psStep = psProcessingSteps->psChild; psStep;
     549          48 :          psStep = psStep->psNext)
     550             :     {
     551          77 :         if (psStep->eType == CXT_Element &&
     552          77 :             strcmp(psStep->pszValue, "Step") == 0)
     553             :         {
     554          77 :             ++iStep;
     555          77 :             const bool bIsFinalStep = (iStep == nStepCount);
     556          77 :             std::vector<double> adfOutNoData;
     557          77 :             if (bIsFinalStep)
     558             :             {
     559             :                 // Initialize adfOutNoData with nodata value of *output* bands
     560             :                 // for final step
     561          75 :                 if (m_outputBandCountProvenance ==
     562             :                     ValueProvenance::FROM_VRTRASTERBAND)
     563             :                 {
     564          13 :                     for (int i = 1; i <= nBands; ++i)
     565             :                     {
     566           8 :                         int bHasVal = FALSE;
     567             :                         const double dfVal =
     568           8 :                             GetRasterBand(i)->GetNoDataValue(&bHasVal);
     569             :                         adfOutNoData.emplace_back(
     570           8 :                             bHasVal ? dfVal
     571           8 :                                     : std::numeric_limits<double>::quiet_NaN());
     572             :                     }
     573             :                 }
     574          70 :                 else if (m_outputBandCountProvenance ==
     575             :                          ValueProvenance::FROM_SOURCE)
     576             :                 {
     577         213 :                     for (int i = 1; i <= m_poSrcDS->GetRasterCount(); ++i)
     578             :                     {
     579         152 :                         int bHasVal = FALSE;
     580             :                         const double dfVal =
     581         152 :                             m_poSrcDS->GetRasterBand(i)->GetNoDataValue(
     582         152 :                                 &bHasVal);
     583             :                         adfOutNoData.emplace_back(
     584         182 :                             bHasVal ? dfVal
     585         182 :                                     : std::numeric_limits<double>::quiet_NaN());
     586             :                     }
     587             :                 }
     588           9 :                 else if (m_outputBandCountProvenance ==
     589             :                          ValueProvenance::USER_PROVIDED)
     590             :                 {
     591           1 :                     adfOutNoData.resize(
     592           1 :                         m_outputBandCountValue,
     593           1 :                         std::numeric_limits<double>::quiet_NaN());
     594             :                 }
     595             :             }
     596          77 :             if (!ParseStep(psStep, bIsFinalStep, eCurrentDT, nCurrentBandCount,
     597             :                            adfNoData, adfOutNoData))
     598          29 :                 return CE_Failure;
     599          48 :             adfNoData = std::move(adfOutNoData);
     600             :         }
     601             :     }
     602             : 
     603          47 :     if (m_aoSteps.empty())
     604             :     {
     605           1 :         CPLError(CE_Failure, CPLE_AppDefined,
     606             :                  "At least one step should be defined");
     607           1 :         return CE_Failure;
     608             :     }
     609             : 
     610          46 :     int nLargestInDTSizeTimesBand = 1;
     611          46 :     int nLargestOutDTSizeTimesBand = 1;
     612          94 :     for (const auto &oStep : m_aoSteps)
     613             :     {
     614             :         const int nInDTSizeTimesBand =
     615          48 :             GDALGetDataTypeSizeBytes(oStep.eInDT) * oStep.nInBands;
     616          48 :         nLargestInDTSizeTimesBand =
     617          48 :             std::max(nLargestInDTSizeTimesBand, nInDTSizeTimesBand);
     618             :         const int nOutDTSizeTimesBand =
     619          48 :             GDALGetDataTypeSizeBytes(oStep.eOutDT) * oStep.nOutBands;
     620          48 :         nLargestOutDTSizeTimesBand =
     621          48 :             std::max(nLargestOutDTSizeTimesBand, nOutDTSizeTimesBand);
     622             :     }
     623          46 :     m_nWorkingBytesPerPixel =
     624          46 :         nLargestInDTSizeTimesBand + nLargestOutDTSizeTimesBand;
     625             : 
     626             :     // Use only up to 40% of RAM to acquire source bands and generate the output
     627             :     // buffer.
     628          46 :     m_nAllowedRAMUsage = CPLGetUsablePhysicalRAM() / 10 * 4;
     629             :     // Only for tests now
     630          46 :     const char *pszMAX_RAM = "VRT_PROCESSED_DATASET_ALLOWED_RAM_USAGE";
     631          46 :     if (const char *pszVal = CPLGetConfigOption(pszMAX_RAM, nullptr))
     632             :     {
     633           3 :         CPL_IGNORE_RET_VAL(
     634           3 :             CPLParseMemorySize(pszVal, &m_nAllowedRAMUsage, nullptr));
     635             :     }
     636             : 
     637          46 :     if (m_nAllowedRAMUsage > 0)
     638             :     {
     639          46 :         bool bBlockSizeModified = false;
     640          78 :         while ((m_nBlockXSize >= 2 || m_nBlockYSize >= 2) &&
     641          71 :                static_cast<GIntBig>(m_nBlockXSize) * m_nBlockYSize >
     642          71 :                    m_nAllowedRAMUsage / m_nWorkingBytesPerPixel)
     643             :         {
     644          32 :             if ((m_nBlockXSize == nRasterXSize ||
     645          30 :                  m_nBlockYSize >= m_nBlockXSize) &&
     646          17 :                 m_nBlockYSize >= 2)
     647             :             {
     648          16 :                 m_nBlockYSize /= 2;
     649             :             }
     650             :             else
     651             :             {
     652          16 :                 m_nBlockXSize /= 2;
     653             :             }
     654          32 :             bBlockSizeModified = true;
     655             :         }
     656          46 :         if (bBlockSizeModified)
     657             :         {
     658           3 :             if (bUserBlockSize)
     659             :             {
     660           0 :                 CPLError(
     661             :                     CE_Warning, CPLE_AppDefined,
     662             :                     "Reducing block size to %d x %d to avoid consuming too "
     663             :                     "much RAM",
     664             :                     m_nBlockXSize, m_nBlockYSize);
     665             :             }
     666             :             else
     667             :             {
     668           3 :                 CPLDebug(
     669             :                     "VRT",
     670             :                     "Reducing block size to %d x %d to avoid consuming too "
     671             :                     "much RAM",
     672             :                     m_nBlockXSize, m_nBlockYSize);
     673             :             }
     674             :         }
     675             :     }
     676             : 
     677          46 :     if (m_outputBandCountProvenance == ValueProvenance::FROM_LAST_STEP)
     678             :     {
     679           7 :         nOutputBandCount = nCurrentBandCount;
     680             :     }
     681          39 :     else if (nOutputBandCount != nCurrentBandCount)
     682             :     {
     683             :         // Should not happen frequently as pixel init functions are expected
     684             :         // to validate that they can accept the number of output bands provided
     685             :         // to them
     686           0 :         CPLError(
     687             :             CE_Failure, CPLE_AppDefined,
     688             :             "Number of output bands of last step (%d) is not consistent with "
     689             :             "number of VRTProcessedRasterBand's (%d)",
     690             :             nCurrentBandCount, nBands);
     691           0 :         return CE_Failure;
     692             :     }
     693             : 
     694          46 :     if (m_outputBandDataTypeProvenance == ValueProvenance::FROM_LAST_STEP)
     695             :     {
     696           3 :         m_outputBandDataTypeValue = eCurrentDT;
     697             :     }
     698             : 
     699          53 :     if (nBands != 0 &&
     700           7 :         (nBands != nOutputBandCount ||
     701           6 :          (m_outputBandDataTypeProvenance == ValueProvenance::FROM_LAST_STEP &&
     702           1 :           m_outputBandDataTypeValue != papoBands[0]->GetRasterDataType())))
     703             :     {
     704           2 :         for (int i = 0; i < nBands; ++i)
     705           1 :             delete papoBands[i];
     706           1 :         CPLFree(papoBands);
     707           1 :         papoBands = nullptr;
     708           1 :         nBands = 0;
     709             :     }
     710             : 
     711         176 :     const auto GetOutputBandType = [this, eCurrentDT](GDALDataType eSourceDT)
     712             :     {
     713          87 :         if (m_outputBandDataTypeProvenance == ValueProvenance::FROM_LAST_STEP)
     714           3 :             return eCurrentDT;
     715          84 :         else if (m_outputBandDataTypeProvenance ==
     716             :                  ValueProvenance::USER_PROVIDED)
     717           5 :             return m_outputBandDataTypeValue;
     718             :         else
     719          79 :             return eSourceDT;
     720          46 :     };
     721             : 
     722          46 :     if (m_outputBandCountProvenance == ValueProvenance::FROM_SOURCE)
     723             :     {
     724         112 :         for (int i = 0; i < m_poSrcDS->GetRasterCount(); ++i)
     725             :         {
     726          79 :             const auto poSrcBand = m_poSrcDS->GetRasterBand(i + 1);
     727             :             const GDALDataType eOutputBandType =
     728          79 :                 GetOutputBandType(poSrcBand->GetRasterDataType());
     729             :             auto poBand =
     730          79 :                 new VRTProcessedRasterBand(this, i + 1, eOutputBandType);
     731          79 :             poBand->CopyCommonInfoFrom(poSrcBand);
     732          79 :             SetBand(i + 1, poBand);
     733             :         }
     734             :     }
     735          13 :     else if (m_outputBandCountProvenance != ValueProvenance::FROM_VRTRASTERBAND)
     736             :     {
     737           8 :         const GDALDataType eOutputBandType = GetOutputBandType(eInDT);
     738          35 :         for (int i = 0; i < nOutputBandCount; ++i)
     739             :         {
     740             :             auto poBand =
     741          27 :                 new VRTProcessedRasterBand(this, i + 1, eOutputBandType);
     742          27 :             SetBand(i + 1, poBand);
     743             :         }
     744             :     }
     745             : 
     746          46 :     if (nBands > 1)
     747          37 :         SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE");
     748             : 
     749          46 :     m_oXMLTree.reset(CPLCloneXMLTree(psTree));
     750             : 
     751          46 :     return CE_None;
     752             : }
     753             : 
     754             : /************************************************************************/
     755             : /*                            ParseStep()                               */
     756             : /************************************************************************/
     757             : 
     758             : /** Parse the current Step node and create a corresponding entry in m_aoSteps.
     759             :  *
     760             :  * @param psStep Step node
     761             :  * @param bIsFinalStep Whether this is the final step.
     762             :  * @param[in,out] eCurrentDT Input data type for this step.
     763             :  *                           Updated to output data type at end of method.
     764             :  * @param[in,out] nCurrentBandCount Input band count for this step.
     765             :  *                                  Updated to output band cout at end of
     766             :  *                                  method.
     767             :  * @param adfInNoData Input nodata values
     768             :  * @param[in,out] adfOutNoData Output nodata values, to be filled by this
     769             :  *                             method. When bIsFinalStep, this is also an
     770             :  *                             input parameter.
     771             :  * @return true on success.
     772             :  */
     773          77 : bool VRTProcessedDataset::ParseStep(const CPLXMLNode *psStep, bool bIsFinalStep,
     774             :                                     GDALDataType &eCurrentDT,
     775             :                                     int &nCurrentBandCount,
     776             :                                     std::vector<double> &adfInNoData,
     777             :                                     std::vector<double> &adfOutNoData)
     778             : {
     779          77 :     const char *pszStepName = CPLGetXMLValue(
     780          77 :         psStep, "name", CPLSPrintf("nr %d", 1 + int(m_aoSteps.size())));
     781          77 :     const char *pszAlgorithm = CPLGetXMLValue(psStep, "Algorithm", nullptr);
     782          77 :     if (!pszAlgorithm)
     783             :     {
     784           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     785             :                  "Step '%s' lacks a Algorithm element", pszStepName);
     786           0 :         return false;
     787             :     }
     788             : 
     789          77 :     const auto &oMapFunctions = GetGlobalMapProcessedDatasetFunc();
     790          77 :     const auto oIterFunc = oMapFunctions.find(pszAlgorithm);
     791          77 :     if (oIterFunc == oMapFunctions.end())
     792             :     {
     793           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     794             :                  "Step '%s' uses unregistered algorithm '%s'", pszStepName,
     795             :                  pszAlgorithm);
     796           0 :         return false;
     797             :     }
     798             : 
     799          77 :     const auto &oFunc = oIterFunc->second;
     800             : 
     801          77 :     if (!oFunc.aeSupportedInputDT.empty())
     802             :     {
     803           0 :         if (std::find(oFunc.aeSupportedInputDT.begin(),
     804             :                       oFunc.aeSupportedInputDT.end(),
     805           0 :                       eCurrentDT) == oFunc.aeSupportedInputDT.end())
     806             :         {
     807           0 :             CPLError(CE_Failure, CPLE_AppDefined,
     808             :                      "Step '%s' (using algorithm '%s') does not "
     809             :                      "support input data type = '%s'",
     810             :                      pszStepName, pszAlgorithm,
     811             :                      GDALGetDataTypeName(eCurrentDT));
     812           0 :             return false;
     813             :         }
     814             :     }
     815             : 
     816          77 :     if (!oFunc.anSupportedInputBandCount.empty())
     817             :     {
     818           0 :         if (std::find(oFunc.anSupportedInputBandCount.begin(),
     819             :                       oFunc.anSupportedInputBandCount.end(),
     820           0 :                       nCurrentBandCount) ==
     821           0 :             oFunc.anSupportedInputBandCount.end())
     822             :         {
     823           0 :             CPLError(CE_Failure, CPLE_AppDefined,
     824             :                      "Step '%s' (using algorithm '%s') does not "
     825             :                      "support input band count = %d",
     826             :                      pszStepName, pszAlgorithm, nCurrentBandCount);
     827           0 :             return false;
     828             :         }
     829             :     }
     830             : 
     831         154 :     Step oStep;
     832          77 :     oStep.osAlgorithm = pszAlgorithm;
     833         154 :     oStep.eInDT = oFunc.eRequestedInputDT != GDT_Unknown
     834          77 :                       ? oFunc.eRequestedInputDT
     835             :                       : eCurrentDT;
     836          77 :     oStep.nInBands = nCurrentBandCount;
     837             : 
     838             :     // Unless modified by pfnInit...
     839          77 :     oStep.eOutDT = oStep.eInDT;
     840             : 
     841          77 :     oStep.adfInNoData = adfInNoData;
     842          77 :     oStep.adfOutNoData = bIsFinalStep ? adfOutNoData : adfInNoData;
     843             : 
     844             :     // Deal with constant arguments
     845          77 :     for (const auto &nameValuePair : oFunc.oMapConstantArguments)
     846             :     {
     847             :         oStep.aosArguments.AddNameValue(nameValuePair.first.c_str(),
     848           0 :                                         nameValuePair.second.c_str());
     849             :     }
     850             : 
     851             :     // Deal with built-in arguments
     852          77 :     if (oFunc.oSetBuiltinArguments.find("nodata") !=
     853         154 :         oFunc.oSetBuiltinArguments.end())
     854             :     {
     855           0 :         int bHasVal = false;
     856           0 :         const auto poSrcFirstBand = m_poSrcDS->GetRasterBand(1);
     857           0 :         const double dfVal = poSrcFirstBand->GetNoDataValue(&bHasVal);
     858           0 :         if (bHasVal)
     859             :         {
     860             :             oStep.aosArguments.AddNameValue("nodata",
     861           0 :                                             CPLSPrintf("%.17g", dfVal));
     862             :         }
     863             :     }
     864             : 
     865          77 :     if (oFunc.oSetBuiltinArguments.find("offset_{band}") !=
     866         154 :         oFunc.oSetBuiltinArguments.end())
     867             :     {
     868           0 :         for (int i = 1; i <= m_poSrcDS->GetRasterCount(); ++i)
     869             :         {
     870           0 :             int bHasVal = false;
     871           0 :             const double dfVal = GetRasterBand(i)->GetOffset(&bHasVal);
     872             :             oStep.aosArguments.AddNameValue(
     873             :                 CPLSPrintf("offset_%d", i),
     874           0 :                 CPLSPrintf("%.17g", bHasVal ? dfVal : 0.0));
     875             :         }
     876             :     }
     877             : 
     878          77 :     if (oFunc.oSetBuiltinArguments.find("scale_{band}") !=
     879         154 :         oFunc.oSetBuiltinArguments.end())
     880             :     {
     881           0 :         for (int i = 1; i <= m_poSrcDS->GetRasterCount(); ++i)
     882             :         {
     883           0 :             int bHasVal = false;
     884           0 :             const double dfVal = GetRasterBand(i)->GetScale(&bHasVal);
     885             :             oStep.aosArguments.AddNameValue(
     886             :                 CPLSPrintf("scale_%d", i),
     887           0 :                 CPLSPrintf("%.17g", bHasVal ? dfVal : 1.0));
     888             :         }
     889             :     }
     890             : 
     891             :     // Parse arguments specified in VRT
     892         154 :     std::set<std::string> oFoundArguments;
     893             : 
     894         421 :     for (const CPLXMLNode *psStepChild = psStep->psChild; psStepChild;
     895         344 :          psStepChild = psStepChild->psNext)
     896             :     {
     897         344 :         if (psStepChild->eType == CXT_Element &&
     898         320 :             strcmp(psStepChild->pszValue, "Argument") == 0)
     899             :         {
     900             :             const char *pszParamName =
     901         243 :                 CPLGetXMLValue(psStepChild, "name", nullptr);
     902         243 :             if (!pszParamName)
     903             :             {
     904           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
     905             :                          "Step '%s' has a Argument without a name attribute",
     906             :                          pszStepName);
     907           0 :                 return false;
     908             :             }
     909         243 :             const char *pszValue = CPLGetXMLValue(psStepChild, nullptr, "");
     910             :             auto oOtherArgIter =
     911         243 :                 oFunc.oOtherArguments.find(CPLString(pszParamName).tolower());
     912         486 :             if (!oFunc.oOtherArguments.empty() &&
     913         486 :                 oOtherArgIter == oFunc.oOtherArguments.end())
     914             :             {
     915             :                 // If we got a parameter name like 'coefficients_1',
     916             :                 // try to fetch the generic 'coefficients_{band}'
     917         296 :                 std::string osParamName(pszParamName);
     918         148 :                 const auto nPos = osParamName.rfind('_');
     919         148 :                 if (nPos != std::string::npos)
     920             :                 {
     921         148 :                     osParamName.resize(nPos + 1);
     922         148 :                     osParamName += "{band}";
     923             :                     oOtherArgIter = oFunc.oOtherArguments.find(
     924         148 :                         CPLString(osParamName).tolower());
     925             :                 }
     926             :             }
     927         243 :             if (oOtherArgIter != oFunc.oOtherArguments.end())
     928             :             {
     929         243 :                 oFoundArguments.insert(oOtherArgIter->first);
     930             : 
     931         243 :                 const std::string &osType = oOtherArgIter->second.osType;
     932         243 :                 if (osType == "boolean")
     933             :                 {
     934           0 :                     if (!EQUAL(pszValue, "true") && !EQUAL(pszValue, "false"))
     935             :                     {
     936           0 :                         CPLError(CE_Failure, CPLE_NotSupported,
     937             :                                  "Step '%s' has a Argument '%s' whose "
     938             :                                  "value '%s' is not a boolean",
     939             :                                  pszStepName, pszParamName, pszValue);
     940           0 :                         return false;
     941             :                     }
     942             :                 }
     943         243 :                 else if (osType == "integer")
     944             :                 {
     945          44 :                     if (CPLGetValueType(pszValue) != CPL_VALUE_INTEGER)
     946             :                     {
     947           0 :                         CPLError(CE_Failure, CPLE_NotSupported,
     948             :                                  "Step '%s' has a Argument '%s' whose "
     949             :                                  "value '%s' is not a integer",
     950             :                                  pszStepName, pszParamName, pszValue);
     951           0 :                         return false;
     952             :                     }
     953             :                 }
     954         199 :                 else if (osType == "double")
     955             :                 {
     956          49 :                     const auto eType = CPLGetValueType(pszValue);
     957          49 :                     if (eType != CPL_VALUE_INTEGER && eType != CPL_VALUE_REAL)
     958             :                     {
     959           0 :                         CPLError(CE_Failure, CPLE_NotSupported,
     960             :                                  "Step '%s' has a Argument '%s' whose "
     961             :                                  "value '%s' is not a double",
     962             :                                  pszStepName, pszParamName, pszValue);
     963           0 :                         return false;
     964             :                     }
     965             :                 }
     966         150 :                 else if (osType == "double_list")
     967             :                 {
     968             :                     const CPLStringList aosTokens(
     969          88 :                         CSLTokenizeString2(pszValue, ",", 0));
     970         405 :                     for (int i = 0; i < aosTokens.size(); ++i)
     971             :                     {
     972         317 :                         const auto eType = CPLGetValueType(aosTokens[i]);
     973         317 :                         if (eType != CPL_VALUE_INTEGER &&
     974             :                             eType != CPL_VALUE_REAL)
     975             :                         {
     976           0 :                             CPLError(CE_Failure, CPLE_NotSupported,
     977             :                                      "Step '%s' has a Argument '%s' "
     978             :                                      "whose value '%s' is not a "
     979             :                                      "comma-separated list of doubles",
     980             :                                      pszStepName, pszParamName, pszValue);
     981           0 :                             return false;
     982             :                         }
     983             :                     }
     984             :                 }
     985          62 :                 else if (osType != "string")
     986             :                 {
     987           0 :                     CPLDebug("VRT", "Unhandled argument type '%s'",
     988             :                              osType.c_str());
     989           0 :                     CPLAssert(0);
     990             :                 }
     991             :             }
     992           0 :             else if (oFunc.bMetadataSpecified &&
     993           0 :                      oFunc.oSetBuiltinArguments.find(
     994           0 :                          CPLString(pszParamName).tolower()) ==
     995           0 :                          oFunc.oSetBuiltinArguments.end() &&
     996           0 :                      oFunc.oMapConstantArguments.find(
     997           0 :                          CPLString(pszParamName).tolower()) ==
     998           0 :                          oFunc.oMapConstantArguments.end())
     999             :             {
    1000           0 :                 CPLError(CE_Warning, CPLE_NotSupported,
    1001             :                          "Step '%s' has a Argument '%s' which is not "
    1002             :                          "supported",
    1003             :                          pszStepName, pszParamName);
    1004             :             }
    1005             : 
    1006         243 :             oStep.aosArguments.AddNameValue(pszParamName, pszValue);
    1007             :         }
    1008             :     }
    1009             : 
    1010             :     // Check that required arguments have been specified
    1011         619 :     for (const auto &oIter : oFunc.oOtherArguments)
    1012             :     {
    1013         693 :         if (oIter.second.bRequired &&
    1014         693 :             oFoundArguments.find(oIter.first) == oFoundArguments.end())
    1015             :         {
    1016           3 :             CPLError(CE_Failure, CPLE_AppDefined,
    1017             :                      "Step '%s' lacks required Argument '%s'", pszStepName,
    1018             :                      oIter.first.c_str());
    1019           3 :             return false;
    1020             :         }
    1021             :     }
    1022             : 
    1023          74 :     if (oFunc.pfnInit)
    1024             :     {
    1025          74 :         double *padfOutNoData = nullptr;
    1026          74 :         if (bIsFinalStep && !adfOutNoData.empty())
    1027             :         {
    1028          64 :             oStep.nOutBands = static_cast<int>(adfOutNoData.size());
    1029          64 :             padfOutNoData = static_cast<double *>(
    1030          64 :                 CPLMalloc(adfOutNoData.size() * sizeof(double)));
    1031          64 :             memcpy(padfOutNoData, adfOutNoData.data(),
    1032          64 :                    adfOutNoData.size() * sizeof(double));
    1033             :         }
    1034             :         else
    1035             :         {
    1036          10 :             oStep.nOutBands = 0;
    1037             :         }
    1038             : 
    1039         148 :         if (oFunc.pfnInit(pszAlgorithm, oFunc.pUserData,
    1040          74 :                           oStep.aosArguments.List(), oStep.nInBands,
    1041             :                           oStep.eInDT, adfInNoData.data(), &(oStep.nOutBands),
    1042             :                           &(oStep.eOutDT), &padfOutNoData, m_osVRTPath.c_str(),
    1043          74 :                           &(oStep.pWorkingData)) != CE_None)
    1044             :         {
    1045          26 :             CPLError(CE_Failure, CPLE_AppDefined,
    1046             :                      "Step '%s' (using algorithm '%s') init() function "
    1047             :                      "failed",
    1048             :                      pszStepName, pszAlgorithm);
    1049          26 :             CPLFree(padfOutNoData);
    1050          26 :             return false;
    1051             :         }
    1052             : 
    1053             :         // Input nodata values may have been modified by pfnInit()
    1054          48 :         oStep.adfInNoData = adfInNoData;
    1055             : 
    1056          48 :         if (padfOutNoData)
    1057             :         {
    1058          88 :             adfOutNoData = std::vector<double>(padfOutNoData,
    1059          44 :                                                padfOutNoData + oStep.nOutBands);
    1060             :         }
    1061             :         else
    1062             :         {
    1063          12 :             adfOutNoData = std::vector<double>(
    1064           8 :                 oStep.nOutBands, std::numeric_limits<double>::quiet_NaN());
    1065             :         }
    1066          48 :         CPLFree(padfOutNoData);
    1067             : 
    1068          48 :         oStep.adfOutNoData = adfOutNoData;
    1069             :     }
    1070             :     else
    1071             :     {
    1072           0 :         oStep.nOutBands = oStep.nInBands;
    1073           0 :         adfOutNoData = oStep.adfOutNoData;
    1074             :     }
    1075             : 
    1076          48 :     eCurrentDT = oStep.eOutDT;
    1077          48 :     nCurrentBandCount = oStep.nOutBands;
    1078             : 
    1079          48 :     m_aoSteps.emplace_back(std::move(oStep));
    1080             : 
    1081          48 :     return true;
    1082             : }
    1083             : 
    1084             : /************************************************************************/
    1085             : /*                           SerializeToXML()                           */
    1086             : /************************************************************************/
    1087             : 
    1088           1 : CPLXMLNode *VRTProcessedDataset::SerializeToXML(const char *pszVRTPathIn)
    1089             : 
    1090             : {
    1091           1 :     CPLXMLNode *psTree = CPLCloneXMLTree(m_oXMLTree.get());
    1092           1 :     if (psTree == nullptr)
    1093           0 :         return psTree;
    1094             : 
    1095             :     /* -------------------------------------------------------------------- */
    1096             :     /*      Remove VRTRasterBand nodes from the original tree and find the  */
    1097             :     /*      last child.                                                     */
    1098             :     /* -------------------------------------------------------------------- */
    1099           1 :     CPLXMLNode *psLastChild = psTree->psChild;
    1100           1 :     CPLXMLNode *psPrevChild = nullptr;
    1101           5 :     while (psLastChild)
    1102             :     {
    1103           5 :         CPLXMLNode *psNextChild = psLastChild->psNext;
    1104           5 :         if (psLastChild->eType == CXT_Element &&
    1105           3 :             strcmp(psLastChild->pszValue, "VRTRasterBand") == 0)
    1106             :         {
    1107           1 :             if (psPrevChild)
    1108           1 :                 psPrevChild->psNext = psNextChild;
    1109             :             else
    1110           0 :                 psTree->psChild = psNextChild;
    1111           1 :             psLastChild->psNext = nullptr;
    1112           1 :             CPLDestroyXMLNode(psLastChild);
    1113           1 :             psLastChild = psPrevChild ? psPrevChild : psTree->psChild;
    1114             :         }
    1115           4 :         else if (!psNextChild)
    1116             :         {
    1117           1 :             break;
    1118             :         }
    1119             :         else
    1120             :         {
    1121           3 :             psPrevChild = psLastChild;
    1122           3 :             psLastChild = psNextChild;
    1123             :         }
    1124             :     }
    1125           1 :     CPLAssert(psLastChild);  // we have at least Input
    1126             : 
    1127             :     /* -------------------------------------------------------------------- */
    1128             :     /*      Serialize bands.                                                */
    1129             :     /* -------------------------------------------------------------------- */
    1130           1 :     bool bHasWarnedAboutRAMUsage = false;
    1131           1 :     size_t nAccRAMUsage = 0;
    1132           2 :     for (int iBand = 0; iBand < nBands; iBand++)
    1133             :     {
    1134             :         CPLXMLNode *psBandTree =
    1135           1 :             static_cast<VRTRasterBand *>(papoBands[iBand])
    1136           2 :                 ->SerializeToXML(pszVRTPathIn, bHasWarnedAboutRAMUsage,
    1137           1 :                                  nAccRAMUsage);
    1138             : 
    1139           1 :         if (psBandTree != nullptr)
    1140             :         {
    1141           1 :             psLastChild->psNext = psBandTree;
    1142           1 :             psLastChild = psBandTree;
    1143             :         }
    1144             :     }
    1145             : 
    1146           1 :     return psTree;
    1147             : }
    1148             : 
    1149             : /************************************************************************/
    1150             : /*                           SerializeToXML()                           */
    1151             : /************************************************************************/
    1152             : 
    1153             : CPLXMLNode *
    1154           1 : VRTProcessedRasterBand::SerializeToXML(const char *pszVRTPathIn,
    1155             :                                        bool &bHasWarnedAboutRAMUsage,
    1156             :                                        size_t &nAccRAMUsage)
    1157             : 
    1158             : {
    1159           1 :     CPLXMLNode *psTree = VRTRasterBand::SerializeToXML(
    1160             :         pszVRTPathIn, bHasWarnedAboutRAMUsage, nAccRAMUsage);
    1161             : 
    1162             :     /* -------------------------------------------------------------------- */
    1163             :     /*      Set subclass.                                                   */
    1164             :     /* -------------------------------------------------------------------- */
    1165           1 :     CPLCreateXMLNode(CPLCreateXMLNode(psTree, CXT_Attribute, "subClass"),
    1166             :                      CXT_Text, "VRTProcessedRasterBand");
    1167             : 
    1168           1 :     return psTree;
    1169             : }
    1170             : 
    1171             : /************************************************************************/
    1172             : /*                            GetBlockSize()                            */
    1173             : /************************************************************************/
    1174             : 
    1175             : /** Return block size */
    1176         119 : void VRTProcessedDataset::GetBlockSize(int *pnBlockXSize,
    1177             :                                        int *pnBlockYSize) const
    1178             : 
    1179             : {
    1180         119 :     *pnBlockXSize = m_nBlockXSize;
    1181         119 :     *pnBlockYSize = m_nBlockYSize;
    1182         119 : }
    1183             : 
    1184             : /************************************************************************/
    1185             : /*                            ProcessRegion()                           */
    1186             : /************************************************************************/
    1187             : 
    1188             : /** Compute pixel values for the specified region.
    1189             :  *
    1190             :  * The output is stored in m_abyInput in a pixel-interleaved way.
    1191             :  */
    1192          60 : bool VRTProcessedDataset::ProcessRegion(int nXOff, int nYOff, int nBufXSize,
    1193             :                                         int nBufYSize,
    1194             :                                         GDALProgressFunc pfnProgress,
    1195             :                                         void *pProgressData)
    1196             : {
    1197             : 
    1198          60 :     CPLAssert(!m_aoSteps.empty());
    1199             : 
    1200          60 :     const size_t nPixelCount = static_cast<size_t>(nBufXSize) * nBufYSize;
    1201             : 
    1202          60 :     const int nFirstBandCount = m_aoSteps.front().nInBands;
    1203          60 :     CPLAssert(nFirstBandCount == m_poSrcDS->GetRasterCount());
    1204          60 :     const GDALDataType eFirstDT = m_aoSteps.front().eInDT;
    1205          60 :     const int nFirstDTSize = GDALGetDataTypeSizeBytes(eFirstDT);
    1206          60 :     auto &abyInput = m_abyInput;
    1207          60 :     auto &abyOutput = m_abyOutput;
    1208             : 
    1209             :     const char *pszInterleave =
    1210          60 :         m_poSrcDS->GetMetadataItem("INTERLEAVE", "IMAGE_STRUCTURE");
    1211          60 :     if (nFirstBandCount > 1 && (!pszInterleave || EQUAL(pszInterleave, "BAND")))
    1212             :     {
    1213             :         // If there are several bands and the source dataset organization
    1214             :         // is apparently band interleaved, then first acquire data in
    1215             :         // a BSQ organization in the abyInput array use in the native
    1216             :         // data type.
    1217             :         // And then transpose it and convert it to the expected data type
    1218             :         // of the first step.
    1219           1 :         const auto eSrcDT = m_poSrcDS->GetRasterBand(1)->GetRasterDataType();
    1220             :         try
    1221             :         {
    1222           2 :             abyInput.resize(nPixelCount * nFirstBandCount *
    1223           1 :                             GDALGetDataTypeSizeBytes(eSrcDT));
    1224             :         }
    1225           0 :         catch (const std::bad_alloc &)
    1226             :         {
    1227           0 :             CPLError(CE_Failure, CPLE_OutOfMemory,
    1228             :                      "Out of memory allocating working buffer");
    1229           0 :             return false;
    1230             :         }
    1231             : 
    1232             :         try
    1233             :         {
    1234           1 :             abyOutput.resize(nPixelCount * nFirstBandCount * nFirstDTSize);
    1235             :         }
    1236           0 :         catch (const std::bad_alloc &)
    1237             :         {
    1238           0 :             CPLError(CE_Failure, CPLE_OutOfMemory,
    1239             :                      "Out of memory allocating working buffer");
    1240           0 :             return false;
    1241             :         }
    1242             : 
    1243             :         GDALRasterIOExtraArg sArg;
    1244           1 :         INIT_RASTERIO_EXTRA_ARG(sArg);
    1245           1 :         sArg.pfnProgress = GDALScaledProgress;
    1246           1 :         sArg.pProgressData =
    1247           1 :             GDALCreateScaledProgress(0, 0.5, pfnProgress, pProgressData);
    1248           1 :         if (sArg.pProgressData == nullptr)
    1249           1 :             sArg.pfnProgress = nullptr;
    1250             : 
    1251           1 :         CPLDebugOnly("VRT", "ProcessRegion(): start RasterIO()");
    1252             :         const bool bOK =
    1253           2 :             m_poSrcDS->RasterIO(GF_Read, nXOff, nYOff, nBufXSize, nBufYSize,
    1254           1 :                                 abyInput.data(), nBufXSize, nBufYSize, eSrcDT,
    1255             :                                 nFirstBandCount, nullptr, 0, 0, 0,
    1256           1 :                                 &sArg) == CE_None;
    1257           1 :         CPLDebugOnly("VRT", "ProcessRegion(): end RasterIO()");
    1258           1 :         GDALDestroyScaledProgress(sArg.pProgressData);
    1259           1 :         if (!bOK)
    1260           0 :             return false;
    1261             : 
    1262           1 :         CPLDebugOnly("VRT", "ProcessRegion(): start GDALTranspose2D()");
    1263           1 :         GDALTranspose2D(abyInput.data(), eSrcDT, abyOutput.data(), eFirstDT,
    1264           1 :                         static_cast<size_t>(nBufXSize) * nBufYSize,
    1265             :                         nFirstBandCount);
    1266           1 :         CPLDebugOnly("VRT", "ProcessRegion(): end GDALTranspose2D()");
    1267             : 
    1268             :         // Swap arrays
    1269           1 :         std::swap(abyInput, abyOutput);
    1270             :     }
    1271             :     else
    1272             :     {
    1273             :         try
    1274             :         {
    1275          59 :             abyInput.resize(nPixelCount * nFirstBandCount * nFirstDTSize);
    1276             :         }
    1277           0 :         catch (const std::bad_alloc &)
    1278             :         {
    1279           0 :             CPLError(CE_Failure, CPLE_OutOfMemory,
    1280             :                      "Out of memory allocating working buffer");
    1281           0 :             return false;
    1282             :         }
    1283             : 
    1284             :         GDALRasterIOExtraArg sArg;
    1285          59 :         INIT_RASTERIO_EXTRA_ARG(sArg);
    1286          59 :         sArg.pfnProgress = GDALScaledProgress;
    1287          59 :         sArg.pProgressData =
    1288          59 :             GDALCreateScaledProgress(0, 0.5, pfnProgress, pProgressData);
    1289          59 :         if (sArg.pProgressData == nullptr)
    1290          59 :             sArg.pfnProgress = nullptr;
    1291             : 
    1292             :         const bool bOK =
    1293         118 :             m_poSrcDS->RasterIO(
    1294          59 :                 GF_Read, nXOff, nYOff, nBufXSize, nBufYSize, abyInput.data(),
    1295             :                 nBufXSize, nBufYSize, eFirstDT, nFirstBandCount, nullptr,
    1296          59 :                 static_cast<GSpacing>(nFirstDTSize) * nFirstBandCount,
    1297          59 :                 static_cast<GSpacing>(nFirstDTSize) * nFirstBandCount *
    1298          59 :                     nBufXSize,
    1299          59 :                 nFirstDTSize, &sArg) == CE_None;
    1300             : 
    1301          59 :         GDALDestroyScaledProgress(sArg.pProgressData);
    1302          59 :         if (!bOK)
    1303           3 :             return false;
    1304             :     }
    1305             : 
    1306          57 :     const double dfSrcXOff = nXOff;
    1307          57 :     const double dfSrcYOff = nYOff;
    1308          57 :     const double dfSrcXSize = nBufXSize;
    1309          57 :     const double dfSrcYSize = nBufYSize;
    1310             : 
    1311             :     double adfSrcGT[6];
    1312          57 :     if (m_poSrcDS->GetGeoTransform(adfSrcGT) != CE_None)
    1313             :     {
    1314          38 :         adfSrcGT[0] = 0;
    1315          38 :         adfSrcGT[1] = 1;
    1316          38 :         adfSrcGT[2] = 0;
    1317          38 :         adfSrcGT[3] = 0;
    1318          38 :         adfSrcGT[4] = 0;
    1319          38 :         adfSrcGT[5] = 1;
    1320             :     }
    1321             : 
    1322          57 :     GDALDataType eLastDT = eFirstDT;
    1323          57 :     const auto &oMapFunctions = GetGlobalMapProcessedDatasetFunc();
    1324             : 
    1325          57 :     int iStep = 0;
    1326         116 :     for (const auto &oStep : m_aoSteps)
    1327             :     {
    1328          59 :         const auto oIterFunc = oMapFunctions.find(oStep.osAlgorithm);
    1329          59 :         CPLAssert(oIterFunc != oMapFunctions.end());
    1330             : 
    1331             :         // Data type adaptation
    1332          59 :         if (eLastDT != oStep.eInDT)
    1333             :         {
    1334             :             try
    1335             :             {
    1336           0 :                 abyOutput.resize(nPixelCount * oStep.nInBands *
    1337           0 :                                  GDALGetDataTypeSizeBytes(oStep.eInDT));
    1338             :             }
    1339           0 :             catch (const std::bad_alloc &)
    1340             :             {
    1341           0 :                 CPLError(CE_Failure, CPLE_OutOfMemory,
    1342             :                          "Out of memory allocating working buffer");
    1343           0 :                 return false;
    1344             :             }
    1345             : 
    1346           0 :             GDALCopyWords64(abyInput.data(), eLastDT,
    1347           0 :                             GDALGetDataTypeSizeBytes(eLastDT), abyOutput.data(),
    1348           0 :                             oStep.eInDT, GDALGetDataTypeSizeBytes(oStep.eInDT),
    1349           0 :                             nPixelCount * oStep.nInBands);
    1350             : 
    1351           0 :             std::swap(abyInput, abyOutput);
    1352             :         }
    1353             : 
    1354             :         try
    1355             :         {
    1356         118 :             abyOutput.resize(nPixelCount * oStep.nOutBands *
    1357          59 :                              GDALGetDataTypeSizeBytes(oStep.eOutDT));
    1358             :         }
    1359           0 :         catch (const std::bad_alloc &)
    1360             :         {
    1361           0 :             CPLError(CE_Failure, CPLE_OutOfMemory,
    1362             :                      "Out of memory allocating working buffer");
    1363           0 :             return false;
    1364             :         }
    1365             : 
    1366          59 :         const auto &oFunc = oIterFunc->second;
    1367         236 :         if (oFunc.pfnProcess(
    1368          59 :                 oStep.osAlgorithm.c_str(), oFunc.pUserData, oStep.pWorkingData,
    1369             :                 oStep.aosArguments.List(), nBufXSize, nBufYSize,
    1370          59 :                 abyInput.data(), abyInput.size(), oStep.eInDT, oStep.nInBands,
    1371          59 :                 oStep.adfInNoData.data(), abyOutput.data(), abyOutput.size(),
    1372          59 :                 oStep.eOutDT, oStep.nOutBands, oStep.adfOutNoData.data(),
    1373             :                 dfSrcXOff, dfSrcYOff, dfSrcXSize, dfSrcYSize, adfSrcGT,
    1374             :                 m_osVRTPath.c_str(),
    1375          59 :                 /*papszExtra=*/nullptr) != CE_None)
    1376             :         {
    1377           0 :             return false;
    1378             :         }
    1379             : 
    1380          59 :         std::swap(abyInput, abyOutput);
    1381          59 :         eLastDT = oStep.eOutDT;
    1382             : 
    1383          59 :         ++iStep;
    1384          59 :         if (pfnProgress &&
    1385           0 :             !pfnProgress(0.5 + 0.5 * iStep / static_cast<int>(m_aoSteps.size()),
    1386             :                          "", pProgressData))
    1387           0 :             return false;
    1388             :     }
    1389             : 
    1390          57 :     return true;
    1391             : }
    1392             : 
    1393             : /************************************************************************/
    1394             : /*                        VRTProcessedRasterBand()                      */
    1395             : /************************************************************************/
    1396             : 
    1397             : /** Constructor */
    1398         119 : VRTProcessedRasterBand::VRTProcessedRasterBand(VRTProcessedDataset *poDSIn,
    1399             :                                                int nBandIn,
    1400         119 :                                                GDALDataType eDataTypeIn)
    1401             : {
    1402         119 :     Initialize(poDSIn->GetRasterXSize(), poDSIn->GetRasterYSize());
    1403             : 
    1404         119 :     poDS = poDSIn;
    1405         119 :     nBand = nBandIn;
    1406         119 :     eAccess = GA_Update;
    1407         119 :     eDataType = eDataTypeIn;
    1408             : 
    1409         119 :     poDSIn->GetBlockSize(&nBlockXSize, &nBlockYSize);
    1410         119 : }
    1411             : 
    1412             : /************************************************************************/
    1413             : /*                           GetOverviewCount()                         */
    1414             : /************************************************************************/
    1415             : 
    1416           6 : int VRTProcessedRasterBand::GetOverviewCount()
    1417             : {
    1418           6 :     auto poVRTDS = cpl::down_cast<VRTProcessedDataset *>(poDS);
    1419           6 :     return static_cast<int>(poVRTDS->m_apoOverviewDatasets.size());
    1420             : }
    1421             : 
    1422             : /************************************************************************/
    1423             : /*                              GetOverview()                           */
    1424             : /************************************************************************/
    1425             : 
    1426          11 : GDALRasterBand *VRTProcessedRasterBand::GetOverview(int iOvr)
    1427             : {
    1428          11 :     auto poVRTDS = cpl::down_cast<VRTProcessedDataset *>(poDS);
    1429          22 :     if (iOvr < 0 ||
    1430          11 :         iOvr >= static_cast<int>(poVRTDS->m_apoOverviewDatasets.size()))
    1431           0 :         return nullptr;
    1432          11 :     return poVRTDS->m_apoOverviewDatasets[iOvr]->GetRasterBand(nBand);
    1433             : }
    1434             : 
    1435             : /************************************************************************/
    1436             : /*                             IReadBlock()                             */
    1437             : /************************************************************************/
    1438             : 
    1439          39 : CPLErr VRTProcessedRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
    1440             :                                           void *pImage)
    1441             : 
    1442             : {
    1443          39 :     auto poVRTDS = cpl::down_cast<VRTProcessedDataset *>(poDS);
    1444             : 
    1445          39 :     int nBufXSize = 0;
    1446          39 :     int nBufYSize = 0;
    1447          39 :     GetActualBlockSize(nBlockXOff, nBlockYOff, &nBufXSize, &nBufYSize);
    1448             : 
    1449          39 :     const int nXPixelOff = nBlockXOff * nBlockXSize;
    1450          39 :     const int nYPixelOff = nBlockYOff * nBlockYSize;
    1451          39 :     if (!poVRTDS->ProcessRegion(nXPixelOff, nYPixelOff, nBufXSize, nBufYSize,
    1452             :                                 nullptr, nullptr))
    1453             :     {
    1454           0 :         return CE_Failure;
    1455             :     }
    1456             : 
    1457          39 :     const int nOutBands = poVRTDS->m_aoSteps.back().nOutBands;
    1458          39 :     CPLAssert(nOutBands == poVRTDS->GetRasterCount());
    1459          39 :     const auto eLastDT = poVRTDS->m_aoSteps.back().eOutDT;
    1460          39 :     const int nLastDTSize = GDALGetDataTypeSizeBytes(eLastDT);
    1461          39 :     const int nDTSize = GDALGetDataTypeSizeBytes(eDataType);
    1462             : 
    1463             :     // Dispatch final output buffer to cached blocks of output bands
    1464         114 :     for (int iDstBand = 0; iDstBand < nOutBands; ++iDstBand)
    1465             :     {
    1466          75 :         GDALRasterBlock *poBlock = nullptr;
    1467             :         GByte *pDst;
    1468          75 :         if (iDstBand + 1 == nBand)
    1469             :         {
    1470          39 :             pDst = static_cast<GByte *>(pImage);
    1471             :         }
    1472             :         else
    1473             :         {
    1474          36 :             auto poOtherBand = poVRTDS->papoBands[iDstBand];
    1475          36 :             poBlock = poOtherBand->TryGetLockedBlockRef(nBlockXOff, nBlockYOff);
    1476          36 :             if (poBlock)
    1477             :             {
    1478           0 :                 poBlock->DropLock();
    1479           0 :                 continue;
    1480             :             }
    1481          72 :             poBlock = poOtherBand->GetLockedBlockRef(
    1482          36 :                 nBlockXOff, nBlockYOff, /* bJustInitialized = */ true);
    1483          36 :             if (!poBlock)
    1484           0 :                 continue;
    1485          36 :             pDst = static_cast<GByte *>(poBlock->GetDataRef());
    1486             :         }
    1487         847 :         for (int iY = 0; iY < nBufYSize; ++iY)
    1488             :         {
    1489        1544 :             GDALCopyWords64(poVRTDS->m_abyInput.data() +
    1490         772 :                                 (iDstBand + static_cast<size_t>(iY) *
    1491         772 :                                                 nBufXSize * nOutBands) *
    1492         772 :                                     nLastDTSize,
    1493             :                             eLastDT, nLastDTSize * nOutBands,
    1494         772 :                             pDst +
    1495         772 :                                 static_cast<size_t>(iY) * nBlockXSize * nDTSize,
    1496             :                             eDataType, nDTSize, nBufXSize);
    1497             :         }
    1498          75 :         if (poBlock)
    1499          36 :             poBlock->DropLock();
    1500             :     }
    1501             : 
    1502          39 :     return CE_None;
    1503             : }
    1504             : 
    1505             : /************************************************************************/
    1506             : /*                VRTProcessedDataset::IRasterIO()                      */
    1507             : /************************************************************************/
    1508             : 
    1509          52 : CPLErr VRTProcessedDataset::IRasterIO(
    1510             :     GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
    1511             :     void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
    1512             :     int nBandCount, BANDMAP_TYPE panBandMap, GSpacing nPixelSpace,
    1513             :     GSpacing nLineSpace, GSpacing nBandSpace, GDALRasterIOExtraArg *psExtraArg)
    1514             : {
    1515             :     // Try to pass the request to the most appropriate overview dataset.
    1516          52 :     if (nBufXSize < nXSize && nBufYSize < nYSize)
    1517             :     {
    1518           2 :         int bTried = FALSE;
    1519           2 :         const CPLErr eErr = TryOverviewRasterIO(
    1520             :             eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
    1521             :             eBufType, nBandCount, panBandMap, nPixelSpace, nLineSpace,
    1522             :             nBandSpace, psExtraArg, &bTried);
    1523           2 :         if (bTried)
    1524           2 :             return eErr;
    1525             :     }
    1526             : 
    1527             :     // Optimize reading of all bands at nominal resolution for BIP-like or
    1528             :     // BSQ-like buffer spacing.
    1529          50 :     if (eRWFlag == GF_Read && nXSize == nBufXSize && nYSize == nBufYSize &&
    1530          49 :         nBandCount == nBands)
    1531             :     {
    1532         474 :         const auto IsSequentialBandMap = [panBandMap, nBandCount]()
    1533             :         {
    1534         236 :             for (int i = 0; i < nBandCount; ++i)
    1535             :             {
    1536         190 :                 if (panBandMap[i] != i + 1)
    1537             :                 {
    1538           2 :                     return false;
    1539             :                 }
    1540             :             }
    1541          46 :             return true;
    1542          48 :         };
    1543             : 
    1544          48 :         const int nBufTypeSize = GDALGetDataTypeSizeBytes(eBufType);
    1545             :         const bool bIsBIPLike =
    1546           9 :             nBandSpace == nBufTypeSize && nPixelSpace == nBandSpace * nBands &&
    1547          57 :             nLineSpace >= nPixelSpace * nBufXSize && IsSequentialBandMap();
    1548          91 :         const bool bIsBSQLike = nPixelSpace == nBufTypeSize &&
    1549          43 :                                 nLineSpace >= nPixelSpace * nBufXSize &&
    1550         134 :                                 nBandSpace >= nLineSpace * nBufYSize &&
    1551          43 :                                 IsSequentialBandMap();
    1552          48 :         if (bIsBIPLike || bIsBSQLike)
    1553             :         {
    1554          46 :             GByte *pabyData = static_cast<GByte *>(pData);
    1555             :             // If acquiring the region of interest in a single time is going
    1556             :             // to consume too much RAM, split in halves.
    1557          46 :             if (m_nAllowedRAMUsage > 0 &&
    1558          46 :                 static_cast<GIntBig>(nBufXSize) * nBufYSize >
    1559          46 :                     m_nAllowedRAMUsage / m_nWorkingBytesPerPixel)
    1560             :             {
    1561          25 :                 if ((nBufXSize == nRasterXSize || nBufYSize >= nBufXSize) &&
    1562             :                     nBufYSize >= 2)
    1563             :                 {
    1564             :                     GDALRasterIOExtraArg sArg;
    1565          12 :                     INIT_RASTERIO_EXTRA_ARG(sArg);
    1566          12 :                     const int nHalfHeight = nBufYSize / 2;
    1567             : 
    1568          12 :                     sArg.pfnProgress = GDALScaledProgress;
    1569          12 :                     sArg.pProgressData = GDALCreateScaledProgress(
    1570             :                         0, 0.5, psExtraArg->pfnProgress,
    1571             :                         psExtraArg->pProgressData);
    1572          12 :                     if (sArg.pProgressData == nullptr)
    1573          12 :                         sArg.pfnProgress = nullptr;
    1574             :                     bool bOK =
    1575          12 :                         IRasterIO(eRWFlag, nXOff, nYOff, nBufXSize, nHalfHeight,
    1576             :                                   pabyData, nBufXSize, nHalfHeight, eBufType,
    1577             :                                   nBandCount, panBandMap, nPixelSpace,
    1578          12 :                                   nLineSpace, nBandSpace, &sArg) == CE_None;
    1579          12 :                     GDALDestroyScaledProgress(sArg.pProgressData);
    1580             : 
    1581          12 :                     if (bOK)
    1582             :                     {
    1583           2 :                         sArg.pfnProgress = GDALScaledProgress;
    1584           2 :                         sArg.pProgressData = GDALCreateScaledProgress(
    1585             :                             0.5, 1, psExtraArg->pfnProgress,
    1586             :                             psExtraArg->pProgressData);
    1587           2 :                         if (sArg.pProgressData == nullptr)
    1588           2 :                             sArg.pfnProgress = nullptr;
    1589           4 :                         bOK = IRasterIO(eRWFlag, nXOff, nYOff + nHalfHeight,
    1590             :                                         nBufXSize, nBufYSize - nHalfHeight,
    1591           2 :                                         pabyData + nHalfHeight * nLineSpace,
    1592             :                                         nBufXSize, nBufYSize - nHalfHeight,
    1593             :                                         eBufType, nBandCount, panBandMap,
    1594             :                                         nPixelSpace, nLineSpace, nBandSpace,
    1595             :                                         &sArg) == CE_None;
    1596           2 :                         GDALDestroyScaledProgress(sArg.pProgressData);
    1597             :                     }
    1598          12 :                     return bOK ? CE_None : CE_Failure;
    1599             :                 }
    1600          13 :                 else if (nBufXSize >= 2)
    1601             :                 {
    1602             :                     GDALRasterIOExtraArg sArg;
    1603          13 :                     INIT_RASTERIO_EXTRA_ARG(sArg);
    1604          13 :                     const int nHalfWidth = nBufXSize / 2;
    1605             : 
    1606          13 :                     sArg.pfnProgress = GDALScaledProgress;
    1607          13 :                     sArg.pProgressData = GDALCreateScaledProgress(
    1608             :                         0, 0.5, psExtraArg->pfnProgress,
    1609             :                         psExtraArg->pProgressData);
    1610          13 :                     if (sArg.pProgressData == nullptr)
    1611          13 :                         sArg.pfnProgress = nullptr;
    1612             :                     bool bOK =
    1613          13 :                         IRasterIO(eRWFlag, nXOff, nYOff, nHalfWidth, nBufYSize,
    1614             :                                   pabyData, nHalfWidth, nBufYSize, eBufType,
    1615             :                                   nBandCount, panBandMap, nPixelSpace,
    1616          13 :                                   nLineSpace, nBandSpace, &sArg) == CE_None;
    1617          13 :                     GDALDestroyScaledProgress(sArg.pProgressData);
    1618             : 
    1619          13 :                     if (bOK)
    1620             :                     {
    1621           3 :                         sArg.pfnProgress = GDALScaledProgress;
    1622           3 :                         sArg.pProgressData = GDALCreateScaledProgress(
    1623             :                             0.5, 1, psExtraArg->pfnProgress,
    1624             :                             psExtraArg->pProgressData);
    1625           3 :                         if (sArg.pProgressData == nullptr)
    1626           3 :                             sArg.pfnProgress = nullptr;
    1627           6 :                         bOK = IRasterIO(eRWFlag, nXOff + nHalfWidth, nYOff,
    1628             :                                         nBufXSize - nHalfWidth, nBufYSize,
    1629           3 :                                         pabyData + nHalfWidth * nPixelSpace,
    1630             :                                         nBufXSize - nHalfWidth, nBufYSize,
    1631             :                                         eBufType, nBandCount, panBandMap,
    1632             :                                         nPixelSpace, nLineSpace, nBandSpace,
    1633             :                                         &sArg) == CE_None;
    1634           3 :                         GDALDestroyScaledProgress(sArg.pProgressData);
    1635             :                     }
    1636          13 :                     return bOK ? CE_None : CE_Failure;
    1637             :                 }
    1638             :             }
    1639             : 
    1640          21 :             if (!ProcessRegion(nXOff, nYOff, nBufXSize, nBufYSize,
    1641             :                                psExtraArg->pfnProgress,
    1642             :                                psExtraArg->pProgressData))
    1643             :             {
    1644           3 :                 return CE_Failure;
    1645             :             }
    1646          18 :             const auto eLastDT = m_aoSteps.back().eOutDT;
    1647          18 :             const int nLastDTSize = GDALGetDataTypeSizeBytes(eLastDT);
    1648          18 :             if (bIsBIPLike)
    1649             :             {
    1650          10 :                 for (int iY = 0; iY < nBufYSize; ++iY)
    1651             :                 {
    1652          21 :                     GDALCopyWords64(
    1653          14 :                         m_abyInput.data() + static_cast<size_t>(iY) * nBands *
    1654           7 :                                                 nBufXSize * nLastDTSize,
    1655           7 :                         eLastDT, nLastDTSize, pabyData + iY * nLineSpace,
    1656             :                         eBufType, GDALGetDataTypeSizeBytes(eBufType),
    1657           7 :                         static_cast<size_t>(nBufXSize) * nBands);
    1658             :                 }
    1659             :             }
    1660             :             else
    1661             :             {
    1662          15 :                 CPLAssert(bIsBSQLike);
    1663          79 :                 for (int iBand = 0; iBand < nBands; ++iBand)
    1664             :                 {
    1665         148 :                     for (int iY = 0; iY < nBufYSize; ++iY)
    1666             :                     {
    1667         168 :                         GDALCopyWords64(
    1668         168 :                             m_abyInput.data() +
    1669          84 :                                 (static_cast<size_t>(iY) * nBands * nBufXSize +
    1670          84 :                                  iBand) *
    1671          84 :                                     nLastDTSize,
    1672          84 :                             eLastDT, nLastDTSize * nBands,
    1673          84 :                             pabyData + iBand * nBandSpace + iY * nLineSpace,
    1674             :                             eBufType, nBufTypeSize, nBufXSize);
    1675             :                     }
    1676             :                 }
    1677             :             }
    1678          18 :             return CE_None;
    1679             :         }
    1680             :     }
    1681             : 
    1682           4 :     return VRTDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, pData,
    1683             :                                  nBufXSize, nBufYSize, eBufType, nBandCount,
    1684             :                                  panBandMap, nPixelSpace, nLineSpace,
    1685           4 :                                  nBandSpace, psExtraArg);
    1686             : }
    1687             : 
    1688             : /*! @endcond */
    1689             : 
    1690             : /************************************************************************/
    1691             : /*                GDALVRTRegisterProcessedDatasetFunc()                 */
    1692             : /************************************************************************/
    1693             : 
    1694             : /** Register a function to be used by VRTProcessedDataset.
    1695             : 
    1696             :  An example of content for pszXMLMetadata is:
    1697             :  \verbatim
    1698             :   <ProcessedDatasetFunctionArgumentsList>
    1699             :      <Argument name='src_nodata' type='double' description='Override input nodata value'/>
    1700             :      <Argument name='dst_nodata' type='double' description='Override output nodata value'/>
    1701             :      <Argument name='replacement_nodata' description='value to substitute to a valid computed value that would be nodata' type='double'/>
    1702             :      <Argument name='dst_intended_datatype' type='string' description='Intented datatype of output (which might be different than the working data type)'/>
    1703             :      <Argument name='coefficients_{band}' description='Comma-separated coefficients for combining bands. First one is constant term' type='double_list' required='true'/>
    1704             :   </ProcessedDatasetFunctionArgumentsList>
    1705             :  \endverbatim
    1706             : 
    1707             :  @param pszFuncName Function name. Must be unique and not null.
    1708             :  @param pUserData User data. May be nullptr. Must remain valid during the
    1709             :                   lifetime of GDAL.
    1710             :  @param pszXMLMetadata XML metadata describing the function arguments. May be
    1711             :                        nullptr if there are no arguments.
    1712             :  @param eRequestedInputDT If the pfnProcess callback only supports a single
    1713             :                           data type, it should be specified in this parameter.
    1714             :                           Otherwise set it to GDT_Unknown.
    1715             :  @param paeSupportedInputDT List of supported input data types. May be nullptr
    1716             :                             if all are supported or if eRequestedInputDT is
    1717             :                             set to a non GDT_Unknown value.
    1718             :  @param nSupportedInputDTSize Size of paeSupportedInputDT
    1719             :  @param panSupportedInputBandCount List of supported band count. May be nullptr
    1720             :                                    if any source band count is supported.
    1721             :  @param nSupportedInputBandCountSize Size of panSupportedInputBandCount
    1722             :  @param pfnInit Initialization function called when a VRTProcessedDataset
    1723             :                 step uses the register function. This initialization function
    1724             :                 will return the output data type, output band count and
    1725             :                 potentially initialize a working structure, typically parsing
    1726             :                 arguments. May be nullptr.
    1727             :                 If not specified, it will be assumed that the input and output
    1728             :                 data types are the same, and that the input number of bands
    1729             :                 and output number of bands are the same.
    1730             :  @param pfnFree Free function that will free the working structure allocated
    1731             :                 by pfnInit. May be nullptr.
    1732             :  @param pfnProcess Processing function called to compute pixel values. Must
    1733             :                    not be nullptr.
    1734             :  @param papszOptions Unused currently. Must be nullptr.
    1735             :  @return CE_None in case of success, error otherwise.
    1736             :  @since 3.9
    1737             :  */
    1738        6900 : CPLErr GDALVRTRegisterProcessedDatasetFunc(
    1739             :     const char *pszFuncName, void *pUserData, const char *pszXMLMetadata,
    1740             :     GDALDataType eRequestedInputDT, const GDALDataType *paeSupportedInputDT,
    1741             :     size_t nSupportedInputDTSize, const int *panSupportedInputBandCount,
    1742             :     size_t nSupportedInputBandCountSize,
    1743             :     GDALVRTProcessedDatasetFuncInit pfnInit,
    1744             :     GDALVRTProcessedDatasetFuncFree pfnFree,
    1745             :     GDALVRTProcessedDatasetFuncProcess pfnProcess,
    1746             :     CPL_UNUSED CSLConstList papszOptions)
    1747             : {
    1748        6900 :     if (pszFuncName == nullptr || pszFuncName[0] == '\0')
    1749             :     {
    1750           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    1751             :                  "pszFuncName should be non-empty");
    1752           0 :         return CE_Failure;
    1753             :     }
    1754             : 
    1755        6900 :     auto &oMap = GetGlobalMapProcessedDatasetFunc();
    1756        6900 :     if (oMap.find(pszFuncName) != oMap.end())
    1757             :     {
    1758           0 :         CPLError(CE_Failure, CPLE_AppDefined, "%s already registered",
    1759             :                  pszFuncName);
    1760           0 :         return CE_Failure;
    1761             :     }
    1762             : 
    1763        6900 :     if (!pfnProcess)
    1764             :     {
    1765           0 :         CPLError(CE_Failure, CPLE_AppDefined, "pfnProcess should not be null");
    1766           0 :         return CE_Failure;
    1767             :     }
    1768             : 
    1769       13800 :     VRTProcessedDatasetFunc oFunc;
    1770        6900 :     oFunc.osFuncName = pszFuncName;
    1771        6900 :     oFunc.pUserData = pUserData;
    1772        6900 :     if (pszXMLMetadata)
    1773             :     {
    1774        6900 :         oFunc.bMetadataSpecified = true;
    1775        6900 :         auto psTree = CPLXMLTreeCloser(CPLParseXMLString(pszXMLMetadata));
    1776        6900 :         if (!psTree)
    1777             :         {
    1778           0 :             CPLError(CE_Failure, CPLE_AppDefined,
    1779             :                      "Cannot parse pszXMLMetadata=%s for %s", pszXMLMetadata,
    1780             :                      pszFuncName);
    1781           0 :             return CE_Failure;
    1782             :         }
    1783        6900 :         const CPLXMLNode *psRoot = CPLGetXMLNode(
    1784             :             psTree.get(), "=ProcessedDatasetFunctionArgumentsList");
    1785        6900 :         if (!psRoot)
    1786             :         {
    1787           0 :             CPLError(CE_Failure, CPLE_AppDefined,
    1788             :                      "No root ProcessedDatasetFunctionArgumentsList element in "
    1789             :                      "pszXMLMetadata=%s for %s",
    1790             :                      pszXMLMetadata, pszFuncName);
    1791           0 :             return CE_Failure;
    1792             :         }
    1793       52440 :         for (const CPLXMLNode *psIter = psRoot->psChild; psIter;
    1794       45540 :              psIter = psIter->psNext)
    1795             :         {
    1796       45540 :             if (psIter->eType == CXT_Element &&
    1797       45540 :                 strcmp(psIter->pszValue, "Argument") == 0)
    1798             :             {
    1799       45540 :                 const char *pszName = CPLGetXMLValue(psIter, "name", nullptr);
    1800       45540 :                 if (!pszName)
    1801             :                 {
    1802           0 :                     CPLError(CE_Failure, CPLE_AppDefined,
    1803             :                              "Missing Argument.name attribute in "
    1804             :                              "pszXMLMetadata=%s for %s",
    1805             :                              pszXMLMetadata, pszFuncName);
    1806           0 :                     return CE_Failure;
    1807             :                 }
    1808       45540 :                 const char *pszType = CPLGetXMLValue(psIter, "type", nullptr);
    1809       45540 :                 if (!pszType)
    1810             :                 {
    1811           0 :                     CPLError(CE_Failure, CPLE_AppDefined,
    1812             :                              "Missing Argument.type attribute in "
    1813             :                              "pszXMLMetadata=%s for %s",
    1814             :                              pszXMLMetadata, pszFuncName);
    1815           0 :                     return CE_Failure;
    1816             :                 }
    1817       45540 :                 if (strcmp(pszType, "constant") == 0)
    1818             :                 {
    1819             :                     const char *pszValue =
    1820           0 :                         CPLGetXMLValue(psIter, "value", nullptr);
    1821           0 :                     if (!pszValue)
    1822             :                     {
    1823           0 :                         CPLError(CE_Failure, CPLE_AppDefined,
    1824             :                                  "Missing Argument.value attribute in "
    1825             :                                  "pszXMLMetadata=%s for %s",
    1826             :                                  pszXMLMetadata, pszFuncName);
    1827           0 :                         return CE_Failure;
    1828             :                     }
    1829           0 :                     oFunc.oMapConstantArguments[CPLString(pszName).tolower()] =
    1830           0 :                         pszValue;
    1831             :                 }
    1832       45540 :                 else if (strcmp(pszType, "builtin") == 0)
    1833             :                 {
    1834           0 :                     if (EQUAL(pszName, "nodata") ||
    1835           0 :                         EQUAL(pszName, "offset_{band}") ||
    1836           0 :                         EQUAL(pszName, "scale_{band}"))
    1837             :                     {
    1838           0 :                         oFunc.oSetBuiltinArguments.insert(
    1839           0 :                             CPLString(pszName).tolower());
    1840             :                     }
    1841             :                     else
    1842             :                     {
    1843           0 :                         CPLError(CE_Failure, CPLE_NotSupported,
    1844             :                                  "Unsupported builtin parameter name %s in "
    1845             :                                  "pszXMLMetadata=%s for %s. Only nodata, "
    1846             :                                  "offset_{band} and scale_{band} are supported",
    1847             :                                  pszName, pszXMLMetadata, pszFuncName);
    1848           0 :                         return CE_Failure;
    1849             :                     }
    1850             :                 }
    1851       45540 :                 else if (strcmp(pszType, "boolean") == 0 ||
    1852       42780 :                          strcmp(pszType, "string") == 0 ||
    1853       33120 :                          strcmp(pszType, "integer") == 0 ||
    1854       24840 :                          strcmp(pszType, "double") == 0 ||
    1855        1380 :                          strcmp(pszType, "double_list") == 0)
    1856             :                 {
    1857       45540 :                     VRTProcessedDatasetFunc::OtherArgument otherArgument;
    1858       45540 :                     otherArgument.bRequired = CPLTestBool(
    1859             :                         CPLGetXMLValue(psIter, "required", "false"));
    1860       45540 :                     otherArgument.osType = pszType;
    1861       91080 :                     oFunc.oOtherArguments[CPLString(pszName).tolower()] =
    1862      136620 :                         std::move(otherArgument);
    1863             :                 }
    1864             :                 else
    1865             :                 {
    1866           0 :                     CPLError(CE_Failure, CPLE_NotSupported,
    1867             :                              "Unsupported type for parameter %s in "
    1868             :                              "pszXMLMetadata=%s for %s. Only boolean, string, "
    1869             :                              "integer, double and double_list are supported",
    1870             :                              pszName, pszXMLMetadata, pszFuncName);
    1871           0 :                     return CE_Failure;
    1872             :                 }
    1873             :             }
    1874             :         }
    1875             :     }
    1876        6900 :     oFunc.eRequestedInputDT = eRequestedInputDT;
    1877        6900 :     if (nSupportedInputDTSize)
    1878             :     {
    1879             :         oFunc.aeSupportedInputDT.insert(
    1880           0 :             oFunc.aeSupportedInputDT.end(), paeSupportedInputDT,
    1881           0 :             paeSupportedInputDT + nSupportedInputDTSize);
    1882             :     }
    1883        6900 :     if (nSupportedInputBandCountSize)
    1884             :     {
    1885             :         oFunc.anSupportedInputBandCount.insert(
    1886           0 :             oFunc.anSupportedInputBandCount.end(), panSupportedInputBandCount,
    1887           0 :             panSupportedInputBandCount + nSupportedInputBandCountSize);
    1888             :     }
    1889        6900 :     oFunc.pfnInit = pfnInit;
    1890        6900 :     oFunc.pfnFree = pfnFree;
    1891        6900 :     oFunc.pfnProcess = pfnProcess;
    1892             : 
    1893        6900 :     oMap[pszFuncName] = std::move(oFunc);
    1894             : 
    1895        6900 :     return CE_None;
    1896             : }

Generated by: LCOV version 1.14