LCOV - code coverage report
Current view: top level - frmts/vrt - vrtprocesseddatasetfunctions.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 575 636 90.4 %
Date: 2024-05-06 22:33:47 Functions: 19 19 100.0 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  Virtual GDAL Datasets
       4             :  * Purpose:  Implementation of VRTProcessedDataset processing functions
       5             :  * Author:   Even Rouault <even.rouault at spatialys.com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2024, Even Rouault <even.rouault at spatialys.com>
       9             :  *
      10             :  * Permission is hereby granted, free of charge, to any person obtaining a
      11             :  * copy of this software and associated documentation files (the "Software"),
      12             :  * to deal in the Software without restriction, including without limitation
      13             :  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
      14             :  * and/or sell copies of the Software, and to permit persons to whom the
      15             :  * Software is furnished to do so, subject to the following conditions:
      16             :  *
      17             :  * The above copyright notice and this permission notice shall be included
      18             :  * in all copies or substantial portions of the Software.
      19             :  *
      20             :  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
      21             :  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
      22             :  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
      23             :  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
      24             :  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
      25             :  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
      26             :  * DEALINGS IN THE SOFTWARE.
      27             :  ****************************************************************************/
      28             : 
      29             : #include "cpl_minixml.h"
      30             : #include "cpl_string.h"
      31             : #include "vrtdataset.h"
      32             : 
      33             : #include <algorithm>
      34             : #include <limits>
      35             : #include <map>
      36             : #include <set>
      37             : #include <vector>
      38             : 
      39             : /************************************************************************/
      40             : /*                               GetDstValue()                          */
      41             : /************************************************************************/
      42             : 
      43             : /** Return a destination value given an initial value, the destination no data
      44             :  * value and its replacement value
      45             :  */
      46         466 : static inline double GetDstValue(double dfVal, double dfDstNoData,
      47             :                                  double dfReplacementDstNodata,
      48             :                                  GDALDataType eIntendedDstDT,
      49             :                                  bool bDstIntendedDTIsInteger)
      50             : {
      51         466 :     if (bDstIntendedDTIsInteger && std::round(dfVal) == dfDstNoData)
      52             :     {
      53           1 :         return dfReplacementDstNodata;
      54             :     }
      55         465 :     else if (eIntendedDstDT == GDT_Float32 &&
      56           0 :              static_cast<float>(dfVal) == static_cast<float>(dfDstNoData))
      57             :     {
      58           0 :         return dfReplacementDstNodata;
      59             :     }
      60         465 :     else if (eIntendedDstDT == GDT_Float64 && dfVal == dfDstNoData)
      61             :     {
      62           1 :         return dfReplacementDstNodata;
      63             :     }
      64             :     else
      65             :     {
      66         464 :         return dfVal;
      67             :     }
      68             : }
      69             : 
      70             : /************************************************************************/
      71             : /*                    BandAffineCombinationData                         */
      72             : /************************************************************************/
      73             : 
      74             : namespace
      75             : {
      76             : /** Working structure for 'BandAffineCombination' builtin function. */
      77             : struct BandAffineCombinationData
      78             : {
      79             :     static constexpr const char *const EXPECTED_SIGNATURE =
      80             :         "BandAffineCombination";
      81             :     //! Signature (to make sure callback functions are called with the right argument)
      82             :     const std::string m_osSignature = EXPECTED_SIGNATURE;
      83             : 
      84             :     /** Replacement nodata value */
      85             :     std::vector<double> m_adfReplacementDstNodata{};
      86             : 
      87             :     /** Intended destination data type. */
      88             :     GDALDataType m_eIntendedDstDT = GDT_Float64;
      89             : 
      90             :     /** Affine transformation coefficients.
      91             :      * m_aadfCoefficients[i][0] is the constant term for the i(th) dst band
      92             :      * m_aadfCoefficients[i][j] is the weight of the j(th) src band for the
      93             :      * i(th) dst vand.
      94             :      * Said otherwise dst[i] = m_aadfCoefficients[i][0] +
      95             :      *      sum(m_aadfCoefficients[i][j + 1] * src[j] for j in 0...nSrcBands-1)
      96             :      */
      97             :     std::vector<std::vector<double>> m_aadfCoefficients{};
      98             : 
      99             :     //! Minimum clamping value.
     100             :     double m_dfClampMin = std::numeric_limits<double>::quiet_NaN();
     101             : 
     102             :     //! Maximum clamping value.
     103             :     double m_dfClampMax = std::numeric_limits<double>::quiet_NaN();
     104             : };
     105             : }  // namespace
     106             : 
     107             : /************************************************************************/
     108             : /*               SetOutputValuesForInNoDataAndOutNoData()               */
     109             : /************************************************************************/
     110             : 
     111          17 : static std::vector<double> SetOutputValuesForInNoDataAndOutNoData(
     112             :     int nInBands, double *padfInNoData, int *pnOutBands,
     113             :     double **ppadfOutNoData, bool bSrcNodataSpecified, double dfSrcNoData,
     114             :     bool bDstNodataSpecified, double dfDstNoData, bool bIsFinalStep)
     115             : {
     116          17 :     if (bSrcNodataSpecified)
     117             :     {
     118           3 :         std::vector<double> adfNoData(nInBands, dfSrcNoData);
     119           3 :         memcpy(padfInNoData, adfNoData.data(),
     120           3 :                adfNoData.size() * sizeof(double));
     121             :     }
     122             : 
     123          17 :     std::vector<double> adfDstNoData;
     124          17 :     if (bDstNodataSpecified)
     125             :     {
     126           3 :         adfDstNoData.resize(*pnOutBands, dfDstNoData);
     127             :     }
     128          14 :     else if (bIsFinalStep)
     129             :     {
     130             :         adfDstNoData =
     131          12 :             std::vector<double>(*ppadfOutNoData, *ppadfOutNoData + *pnOutBands);
     132             :     }
     133             :     else
     134             :     {
     135             :         adfDstNoData =
     136           2 :             std::vector<double>(padfInNoData, padfInNoData + nInBands);
     137           2 :         adfDstNoData.resize(*pnOutBands, *padfInNoData);
     138             :     }
     139             : 
     140          17 :     if (*ppadfOutNoData == nullptr)
     141             :     {
     142           2 :         *ppadfOutNoData =
     143           2 :             static_cast<double *>(CPLMalloc(*pnOutBands * sizeof(double)));
     144             :     }
     145          17 :     memcpy(*ppadfOutNoData, adfDstNoData.data(), *pnOutBands * sizeof(double));
     146             : 
     147          17 :     return adfDstNoData;
     148             : }
     149             : 
     150             : /************************************************************************/
     151             : /*                    BandAffineCombinationInit()                       */
     152             : /************************************************************************/
     153             : 
     154             : /** Init function for 'BandAffineCombination' builtin function. */
     155          12 : static CPLErr BandAffineCombinationInit(
     156             :     const char * /*pszFuncName*/, void * /*pUserData*/,
     157             :     CSLConstList papszFunctionArgs, int nInBands, GDALDataType eInDT,
     158             :     double *padfInNoData, int *pnOutBands, GDALDataType *peOutDT,
     159             :     double **ppadfOutNoData, const char * /* pszVRTPath */,
     160             :     VRTPDWorkingDataPtr *ppWorkingData)
     161             : {
     162          12 :     CPLAssert(eInDT == GDT_Float64);
     163             : 
     164          12 :     *peOutDT = eInDT;
     165          12 :     *ppWorkingData = nullptr;
     166             : 
     167          24 :     auto data = std::make_unique<BandAffineCombinationData>();
     168             : 
     169          24 :     std::map<int, std::vector<double>> oMapCoefficients{};
     170          12 :     double dfSrcNoData = std::numeric_limits<double>::quiet_NaN();
     171          12 :     bool bSrcNodataSpecified = false;
     172          12 :     double dfDstNoData = std::numeric_limits<double>::quiet_NaN();
     173          12 :     bool bDstNodataSpecified = false;
     174          12 :     double dfReplacementDstNodata = std::numeric_limits<double>::quiet_NaN();
     175          12 :     bool bReplacementDstNodataSpecified = false;
     176             : 
     177          65 :     for (const auto &[pszKey, pszValue] :
     178          76 :          cpl::IterateNameValue(papszFunctionArgs))
     179             :     {
     180          33 :         if (EQUAL(pszKey, "src_nodata"))
     181             :         {
     182           2 :             bSrcNodataSpecified = true;
     183           2 :             dfSrcNoData = CPLAtof(pszValue);
     184             :         }
     185          31 :         else if (EQUAL(pszKey, "dst_nodata"))
     186             :         {
     187           2 :             bDstNodataSpecified = true;
     188           2 :             dfDstNoData = CPLAtof(pszValue);
     189             :         }
     190          29 :         else if (EQUAL(pszKey, "replacement_nodata"))
     191             :         {
     192           1 :             bReplacementDstNodataSpecified = true;
     193           1 :             dfReplacementDstNodata = CPLAtof(pszValue);
     194             :         }
     195          28 :         else if (EQUAL(pszKey, "dst_intended_datatype"))
     196             :         {
     197           1 :             for (GDALDataType eDT = GDT_Byte; eDT < GDT_TypeCount;
     198           0 :                  eDT = static_cast<GDALDataType>(eDT + 1))
     199             :             {
     200           1 :                 if (EQUAL(GDALGetDataTypeName(eDT), pszValue))
     201             :                 {
     202           1 :                     data->m_eIntendedDstDT = eDT;
     203           1 :                     break;
     204             :                 }
     205             :             }
     206             :         }
     207          27 :         else if (STARTS_WITH_CI(pszKey, "coefficients_"))
     208             :         {
     209          25 :             const int nTargetBand = atoi(pszKey + strlen("coefficients_"));
     210          25 :             if (nTargetBand <= 0 || nTargetBand > 65536)
     211             :             {
     212           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
     213             :                          "Invalid band in argument '%s'", pszKey);
     214           1 :                 return CE_Failure;
     215             :             }
     216          25 :             const CPLStringList aosTokens(CSLTokenizeString2(pszValue, ",", 0));
     217          25 :             if (aosTokens.size() != 1 + nInBands)
     218             :             {
     219           1 :                 CPLError(CE_Failure, CPLE_AppDefined,
     220             :                          "Argument %s has %d values, whereas %d are expected",
     221             :                          pszKey, aosTokens.size(), 1 + nInBands);
     222           1 :                 return CE_Failure;
     223             :             }
     224          24 :             std::vector<double> adfValues;
     225         110 :             for (int i = 0; i < aosTokens.size(); ++i)
     226             :             {
     227          86 :                 adfValues.push_back(CPLAtof(aosTokens[i]));
     228             :             }
     229          24 :             oMapCoefficients[nTargetBand - 1] = std::move(adfValues);
     230             :         }
     231           2 :         else if (EQUAL(pszKey, "min"))
     232             :         {
     233           1 :             data->m_dfClampMin = CPLAtof(pszValue);
     234             :         }
     235           1 :         else if (EQUAL(pszKey, "max"))
     236             :         {
     237           1 :             data->m_dfClampMax = CPLAtof(pszValue);
     238             :         }
     239             :         else
     240             :         {
     241           0 :             CPLError(CE_Warning, CPLE_AppDefined,
     242             :                      "Unrecognized argument name %s. Ignored", pszKey);
     243             :         }
     244             :     }
     245             : 
     246          11 :     const bool bIsFinalStep = *pnOutBands != 0;
     247          11 :     if (bIsFinalStep)
     248             :     {
     249           9 :         if (*pnOutBands != static_cast<int>(oMapCoefficients.size()))
     250             :         {
     251           1 :             CPLError(CE_Failure, CPLE_AppDefined,
     252             :                      "Final step expect %d bands, but only %d coefficient_XX "
     253             :                      "are provided",
     254           1 :                      *pnOutBands, static_cast<int>(oMapCoefficients.size()));
     255           1 :             return CE_Failure;
     256             :         }
     257             :     }
     258             :     else
     259             :     {
     260           2 :         *pnOutBands = static_cast<int>(oMapCoefficients.size());
     261             :     }
     262             : 
     263             :     const std::vector<double> adfDstNoData =
     264             :         SetOutputValuesForInNoDataAndOutNoData(
     265             :             nInBands, padfInNoData, pnOutBands, ppadfOutNoData,
     266             :             bSrcNodataSpecified, dfSrcNoData, bDstNodataSpecified, dfDstNoData,
     267          20 :             bIsFinalStep);
     268             : 
     269          10 :     if (bReplacementDstNodataSpecified)
     270             :     {
     271           1 :         data->m_adfReplacementDstNodata.resize(*pnOutBands,
     272             :                                                dfReplacementDstNodata);
     273             :     }
     274             :     else
     275             :     {
     276          30 :         for (double dfVal : adfDstNoData)
     277             :         {
     278          21 :             data->m_adfReplacementDstNodata.emplace_back(
     279          21 :                 GDALGetNoDataReplacementValue(data->m_eIntendedDstDT, dfVal));
     280             :         }
     281             :     }
     282             : 
     283             :     // Check we have a set of coefficient for all output bands and
     284             :     // convert the map to a vector
     285          32 :     for (auto &oIter : oMapCoefficients)
     286             :     {
     287          23 :         const int iExpected = static_cast<int>(data->m_aadfCoefficients.size());
     288          23 :         if (oIter.first != iExpected)
     289             :         {
     290           1 :             CPLError(CE_Failure, CPLE_AppDefined,
     291             :                      "Argument coefficients_%d is missing", iExpected + 1);
     292           1 :             return CE_Failure;
     293             :         }
     294          22 :         data->m_aadfCoefficients.emplace_back(std::move(oIter.second));
     295             :     }
     296           9 :     *ppWorkingData = data.release();
     297           9 :     return CE_None;
     298             : }
     299             : 
     300             : /************************************************************************/
     301             : /*                    BandAffineCombinationFree()                       */
     302             : /************************************************************************/
     303             : 
     304             : /** Free function for 'BandAffineCombination' builtin function. */
     305           9 : static void BandAffineCombinationFree(const char * /*pszFuncName*/,
     306             :                                       void * /*pUserData*/,
     307             :                                       VRTPDWorkingDataPtr pWorkingData)
     308             : {
     309           9 :     BandAffineCombinationData *data =
     310             :         static_cast<BandAffineCombinationData *>(pWorkingData);
     311           9 :     CPLAssert(data->m_osSignature ==
     312             :               BandAffineCombinationData::EXPECTED_SIGNATURE);
     313           9 :     CPL_IGNORE_RET_VAL(data->m_osSignature);
     314           9 :     delete data;
     315           9 : }
     316             : 
     317             : /************************************************************************/
     318             : /*                    BandAffineCombinationProcess()                    */
     319             : /************************************************************************/
     320             : 
     321             : /** Processing function for 'BandAffineCombination' builtin function. */
     322           9 : static CPLErr BandAffineCombinationProcess(
     323             :     const char * /*pszFuncName*/, void * /*pUserData*/,
     324             :     VRTPDWorkingDataPtr pWorkingData, CSLConstList /* papszFunctionArgs*/,
     325             :     int nBufXSize, int nBufYSize, const void *pInBuffer, size_t nInBufferSize,
     326             :     GDALDataType eInDT, int nInBands, const double *CPL_RESTRICT padfInNoData,
     327             :     void *pOutBuffer, size_t nOutBufferSize, GDALDataType eOutDT, int nOutBands,
     328             :     const double *CPL_RESTRICT padfOutNoData, double /*dfSrcXOff*/,
     329             :     double /*dfSrcYOff*/, double /*dfSrcXSize*/, double /*dfSrcYSize*/,
     330             :     const double /*adfSrcGT*/[], const char * /* pszVRTPath */,
     331             :     CSLConstList /*papszExtra*/)
     332             : {
     333           9 :     const size_t nElts = static_cast<size_t>(nBufXSize) * nBufYSize;
     334             : 
     335           9 :     CPL_IGNORE_RET_VAL(eInDT);
     336           9 :     CPLAssert(eInDT == GDT_Float64);
     337           9 :     CPL_IGNORE_RET_VAL(eOutDT);
     338           9 :     CPLAssert(eOutDT == GDT_Float64);
     339           9 :     CPL_IGNORE_RET_VAL(nInBufferSize);
     340           9 :     CPLAssert(nInBufferSize == nElts * nInBands * sizeof(double));
     341           9 :     CPL_IGNORE_RET_VAL(nOutBufferSize);
     342           9 :     CPLAssert(nOutBufferSize == nElts * nOutBands * sizeof(double));
     343             : 
     344           9 :     const BandAffineCombinationData *data =
     345             :         static_cast<BandAffineCombinationData *>(pWorkingData);
     346           9 :     CPLAssert(data->m_osSignature ==
     347             :               BandAffineCombinationData::EXPECTED_SIGNATURE);
     348           9 :     const double *CPL_RESTRICT padfSrc = static_cast<const double *>(pInBuffer);
     349           9 :     double *CPL_RESTRICT padfDst = static_cast<double *>(pOutBuffer);
     350             :     const bool bDstIntendedDTIsInteger =
     351           9 :         GDALDataTypeIsInteger(data->m_eIntendedDstDT);
     352           9 :     const double dfClampMin = data->m_dfClampMin;
     353           9 :     const double dfClampMax = data->m_dfClampMax;
     354         171 :     for (size_t i = 0; i < nElts; ++i)
     355             :     {
     356         634 :         for (int iDst = 0; iDst < nOutBands; ++iDst)
     357             :         {
     358         472 :             const auto &adfCoefficients = data->m_aadfCoefficients[iDst];
     359         472 :             double dfVal = adfCoefficients[0];
     360         472 :             bool bSetNoData = false;
     361        1856 :             for (int iSrc = 0; iSrc < nInBands; ++iSrc)
     362             :             {
     363             :                 // written this way to work with a NaN value
     364        1390 :                 if (!(padfSrc[iSrc] != padfInNoData[iSrc]))
     365             :                 {
     366           6 :                     bSetNoData = true;
     367           6 :                     break;
     368             :                 }
     369        1384 :                 dfVal += adfCoefficients[iSrc + 1] * padfSrc[iSrc];
     370             :             }
     371         472 :             if (bSetNoData)
     372             :             {
     373           6 :                 *padfDst = padfOutNoData[iDst];
     374             :             }
     375             :             else
     376             :             {
     377         932 :                 double dfDstVal = GetDstValue(
     378         466 :                     dfVal, padfOutNoData[iDst],
     379         466 :                     data->m_adfReplacementDstNodata[iDst],
     380         466 :                     data->m_eIntendedDstDT, bDstIntendedDTIsInteger);
     381         466 :                 if (dfDstVal < dfClampMin)
     382           1 :                     dfDstVal = dfClampMin;
     383         466 :                 if (dfDstVal > dfClampMax)
     384           1 :                     dfDstVal = dfClampMax;
     385         466 :                 *padfDst = dfDstVal;
     386             :             }
     387         472 :             ++padfDst;
     388             :         }
     389         162 :         padfSrc += nInBands;
     390             :     }
     391             : 
     392           9 :     return CE_None;
     393             : }
     394             : 
     395             : /************************************************************************/
     396             : /*                                LUTData                               */
     397             : /************************************************************************/
     398             : 
     399             : namespace
     400             : {
     401             : /** Working structure for 'LUT' builtin function. */
     402             : struct LUTData
     403             : {
     404             :     static constexpr const char *const EXPECTED_SIGNATURE = "LUT";
     405             :     //! Signature (to make sure callback functions are called with the right argument)
     406             :     const std::string m_osSignature = EXPECTED_SIGNATURE;
     407             : 
     408             :     //! m_aadfLUTInputs[i][j] is the j(th) input value for that LUT of band i.
     409             :     std::vector<std::vector<double>> m_aadfLUTInputs{};
     410             : 
     411             :     //! m_aadfLUTOutputs[i][j] is the j(th) output value for that LUT of band i.
     412             :     std::vector<std::vector<double>> m_aadfLUTOutputs{};
     413             : 
     414             :     /************************************************************************/
     415             :     /*                              LookupValue()                           */
     416             :     /************************************************************************/
     417             : 
     418          18 :     double LookupValue(int iBand, double dfInput) const
     419             :     {
     420          18 :         const auto &adfInput = m_aadfLUTInputs[iBand];
     421          18 :         const auto &afdOutput = m_aadfLUTOutputs[iBand];
     422             : 
     423             :         // Find the index of the first element in the LUT input array that
     424             :         // is not smaller than the input value.
     425             :         int i = static_cast<int>(
     426          18 :             std::lower_bound(adfInput.data(), adfInput.data() + adfInput.size(),
     427          18 :                              dfInput) -
     428          18 :             adfInput.data());
     429             : 
     430          18 :         if (i == 0)
     431           6 :             return afdOutput[0];
     432             : 
     433             :         // If the index is beyond the end of the LUT input array, the input
     434             :         // value is larger than all the values in the array.
     435          12 :         if (i == static_cast<int>(adfInput.size()))
     436           6 :             return afdOutput.back();
     437             : 
     438           6 :         if (adfInput[i] == dfInput)
     439           0 :             return afdOutput[i];
     440             : 
     441             :         // Otherwise, interpolate.
     442           6 :         return afdOutput[i - 1] + (dfInput - adfInput[i - 1]) *
     443           6 :                                       ((afdOutput[i] - afdOutput[i - 1]) /
     444           6 :                                        (adfInput[i] - adfInput[i - 1]));
     445             :     }
     446             : };
     447             : }  // namespace
     448             : 
     449             : /************************************************************************/
     450             : /*                                LUTInit()                             */
     451             : /************************************************************************/
     452             : 
     453             : /** Init function for 'LUT' builtin function. */
     454           6 : static CPLErr LUTInit(const char * /*pszFuncName*/, void * /*pUserData*/,
     455             :                       CSLConstList papszFunctionArgs, int nInBands,
     456             :                       GDALDataType eInDT, double *padfInNoData, int *pnOutBands,
     457             :                       GDALDataType *peOutDT, double **ppadfOutNoData,
     458             :                       const char * /* pszVRTPath */,
     459             :                       VRTPDWorkingDataPtr *ppWorkingData)
     460             : {
     461           6 :     CPLAssert(eInDT == GDT_Float64);
     462             : 
     463           6 :     const bool bIsFinalStep = *pnOutBands != 0;
     464           6 :     *peOutDT = eInDT;
     465           6 :     *ppWorkingData = nullptr;
     466             : 
     467           6 :     if (!bIsFinalStep)
     468             :     {
     469           0 :         *pnOutBands = nInBands;
     470             :     }
     471             : 
     472          12 :     auto data = std::make_unique<LUTData>();
     473             : 
     474           6 :     double dfSrcNoData = std::numeric_limits<double>::quiet_NaN();
     475           6 :     bool bSrcNodataSpecified = false;
     476           6 :     double dfDstNoData = std::numeric_limits<double>::quiet_NaN();
     477           6 :     bool bDstNodataSpecified = false;
     478             : 
     479          12 :     std::map<int, std::pair<std::vector<double>, std::vector<double>>> oMap{};
     480             : 
     481          22 :     for (const auto &[pszKey, pszValue] :
     482          26 :          cpl::IterateNameValue(papszFunctionArgs))
     483             :     {
     484          12 :         if (EQUAL(pszKey, "src_nodata"))
     485             :         {
     486           1 :             bSrcNodataSpecified = true;
     487           1 :             dfSrcNoData = CPLAtof(pszValue);
     488             :         }
     489          11 :         else if (EQUAL(pszKey, "dst_nodata"))
     490             :         {
     491           1 :             bDstNodataSpecified = true;
     492           1 :             dfDstNoData = CPLAtof(pszValue);
     493             :         }
     494          10 :         else if (STARTS_WITH_CI(pszKey, "lut_"))
     495             :         {
     496          10 :             const int nBand = atoi(pszKey + strlen("lut_"));
     497          10 :             if (nBand <= 0 || nBand > nInBands)
     498             :             {
     499           1 :                 CPLError(CE_Failure, CPLE_AppDefined,
     500             :                          "Invalid band in argument '%s'", pszKey);
     501           2 :                 return CE_Failure;
     502             :             }
     503           9 :             const CPLStringList aosTokens(CSLTokenizeString2(pszValue, ",", 0));
     504           9 :             std::vector<double> adfInputValues;
     505           9 :             std::vector<double> adfOutputValues;
     506          26 :             for (int i = 0; i < aosTokens.size(); ++i)
     507             :             {
     508             :                 const CPLStringList aosTokens2(
     509          18 :                     CSLTokenizeString2(aosTokens[i], ":", 0));
     510          18 :                 if (aosTokens2.size() != 2)
     511             :                 {
     512           1 :                     CPLError(CE_Failure, CPLE_AppDefined,
     513             :                              "Invalid value for argument '%s'", pszKey);
     514           1 :                     return CE_Failure;
     515             :                 }
     516          17 :                 adfInputValues.push_back(CPLAtof(aosTokens2[0]));
     517          17 :                 adfOutputValues.push_back(CPLAtof(aosTokens2[1]));
     518             :             }
     519          16 :             oMap[nBand - 1] = std::pair(std::move(adfInputValues),
     520          16 :                                         std::move(adfOutputValues));
     521             :         }
     522             :         else
     523             :         {
     524           0 :             CPLError(CE_Warning, CPLE_AppDefined,
     525             :                      "Unrecognized argument name %s. Ignored", pszKey);
     526             :         }
     527             :     }
     528             : 
     529           4 :     SetOutputValuesForInNoDataAndOutNoData(
     530             :         nInBands, padfInNoData, pnOutBands, ppadfOutNoData, bSrcNodataSpecified,
     531             :         dfSrcNoData, bDstNodataSpecified, dfDstNoData, bIsFinalStep);
     532             : 
     533           4 :     int iExpected = 0;
     534             :     // Check we have values for all bands and convert to vector
     535          11 :     for (auto &oIter : oMap)
     536             :     {
     537           7 :         if (oIter.first != iExpected)
     538             :         {
     539           0 :             CPLError(CE_Failure, CPLE_AppDefined, "Argument lut_%d is missing",
     540             :                      iExpected + 1);
     541           0 :             return CE_Failure;
     542             :         }
     543           7 :         ++iExpected;
     544           7 :         data->m_aadfLUTInputs.emplace_back(std::move(oIter.second.first));
     545           7 :         data->m_aadfLUTOutputs.emplace_back(std::move(oIter.second.second));
     546             :     }
     547             : 
     548           4 :     if (static_cast<int>(oMap.size()) < *pnOutBands)
     549             :     {
     550           1 :         CPLError(CE_Failure, CPLE_AppDefined, "Missing lut_XX element(s)");
     551           1 :         return CE_Failure;
     552             :     }
     553             : 
     554           3 :     *ppWorkingData = data.release();
     555           3 :     return CE_None;
     556             : }
     557             : 
     558             : /************************************************************************/
     559             : /*                                LUTFree()                             */
     560             : /************************************************************************/
     561             : 
     562             : /** Free function for 'LUT' builtin function. */
     563           3 : static void LUTFree(const char * /*pszFuncName*/, void * /*pUserData*/,
     564             :                     VRTPDWorkingDataPtr pWorkingData)
     565             : {
     566           3 :     LUTData *data = static_cast<LUTData *>(pWorkingData);
     567           3 :     CPLAssert(data->m_osSignature == LUTData::EXPECTED_SIGNATURE);
     568           3 :     CPL_IGNORE_RET_VAL(data->m_osSignature);
     569           3 :     delete data;
     570           3 : }
     571             : 
     572             : /************************************************************************/
     573             : /*                             LUTProcess()                             */
     574             : /************************************************************************/
     575             : 
     576             : /** Processing function for 'LUT' builtin function. */
     577             : static CPLErr
     578           3 : LUTProcess(const char * /*pszFuncName*/, void * /*pUserData*/,
     579             :            VRTPDWorkingDataPtr pWorkingData,
     580             :            CSLConstList /* papszFunctionArgs*/, int nBufXSize, int nBufYSize,
     581             :            const void *pInBuffer, size_t nInBufferSize, GDALDataType eInDT,
     582             :            int nInBands, const double *CPL_RESTRICT padfInNoData,
     583             :            void *pOutBuffer, size_t nOutBufferSize, GDALDataType eOutDT,
     584             :            int nOutBands, const double *CPL_RESTRICT padfOutNoData,
     585             :            double /*dfSrcXOff*/, double /*dfSrcYOff*/, double /*dfSrcXSize*/,
     586             :            double /*dfSrcYSize*/, const double /*adfSrcGT*/[],
     587             :            const char * /* pszVRTPath */, CSLConstList /*papszExtra*/)
     588             : {
     589           3 :     const size_t nElts = static_cast<size_t>(nBufXSize) * nBufYSize;
     590             : 
     591           3 :     CPL_IGNORE_RET_VAL(eInDT);
     592           3 :     CPLAssert(eInDT == GDT_Float64);
     593           3 :     CPL_IGNORE_RET_VAL(eOutDT);
     594           3 :     CPLAssert(eOutDT == GDT_Float64);
     595           3 :     CPL_IGNORE_RET_VAL(nInBufferSize);
     596           3 :     CPLAssert(nInBufferSize == nElts * nInBands * sizeof(double));
     597           3 :     CPL_IGNORE_RET_VAL(nOutBufferSize);
     598           3 :     CPLAssert(nOutBufferSize == nElts * nOutBands * sizeof(double));
     599           3 :     CPLAssert(nInBands == nOutBands);
     600           3 :     CPL_IGNORE_RET_VAL(nOutBands);
     601             : 
     602           3 :     const LUTData *data = static_cast<LUTData *>(pWorkingData);
     603           3 :     CPLAssert(data->m_osSignature == LUTData::EXPECTED_SIGNATURE);
     604           3 :     const double *CPL_RESTRICT padfSrc = static_cast<const double *>(pInBuffer);
     605           3 :     double *CPL_RESTRICT padfDst = static_cast<double *>(pOutBuffer);
     606          14 :     for (size_t i = 0; i < nElts; ++i)
     607             :     {
     608          33 :         for (int iBand = 0; iBand < nInBands; ++iBand)
     609             :         {
     610             :             // written this way to work with a NaN value
     611          22 :             if (!(*padfSrc != padfInNoData[iBand]))
     612           4 :                 *padfDst = padfOutNoData[iBand];
     613             :             else
     614          18 :                 *padfDst = data->LookupValue(iBand, *padfSrc);
     615          22 :             ++padfSrc;
     616          22 :             ++padfDst;
     617             :         }
     618             :     }
     619             : 
     620           3 :     return CE_None;
     621             : }
     622             : 
     623             : /************************************************************************/
     624             : /*                        LocalScaleOffsetData                          */
     625             : /************************************************************************/
     626             : 
     627             : namespace
     628             : {
     629             : /** Working structure for 'LocalScaleOffset' builtin function. */
     630             : struct LocalScaleOffsetData
     631             : {
     632             :     static constexpr const char *const EXPECTED_SIGNATURE = "LocalScaleOffset";
     633             :     //! Signature (to make sure callback functions are called with the right argument)
     634             :     const std::string m_osSignature = EXPECTED_SIGNATURE;
     635             : 
     636             :     //! Nodata value for gain dataset(s)
     637             :     double m_dfGainNodata = std::numeric_limits<double>::quiet_NaN();
     638             : 
     639             :     //! Nodata value for offset dataset(s)
     640             :     double m_dfOffsetNodata = std::numeric_limits<double>::quiet_NaN();
     641             : 
     642             :     //! Minimum clamping value.
     643             :     double m_dfClampMin = std::numeric_limits<double>::quiet_NaN();
     644             : 
     645             :     //! Maximum clamping value.
     646             :     double m_dfClampMax = std::numeric_limits<double>::quiet_NaN();
     647             : 
     648             :     //! Map from gain/offset dataset name to datasets
     649             :     std::map<std::string, std::unique_ptr<GDALDataset>> m_oDatasetMap{};
     650             : 
     651             :     //! Vector of size nInBands that point to the raster band from which to read gains.
     652             :     std::vector<GDALRasterBand *> m_oGainBands{};
     653             : 
     654             :     //! Vector of size nInBands that point to the raster band from which to read offsets.
     655             :     std::vector<GDALRasterBand *> m_oOffsetBands{};
     656             : 
     657             :     //! Working buffer that contain gain values.
     658             :     std::vector<VRTProcessedDataset::NoInitByte> m_abyGainBuffer{};
     659             : 
     660             :     //! Working buffer that contain offset values.
     661             :     std::vector<VRTProcessedDataset::NoInitByte> m_abyOffsetBuffer{};
     662             : };
     663             : }  // namespace
     664             : 
     665             : /************************************************************************/
     666             : /*                           CheckAllBands()                            */
     667             : /************************************************************************/
     668             : 
     669             : /** Return true if the key of oMap is the sequence of all integers between
     670             :  * 0 and nExpectedBandCount-1.
     671             :  */
     672             : template <class T>
     673          20 : static bool CheckAllBands(const std::map<int, T> &oMap, int nExpectedBandCount)
     674             : {
     675          20 :     int iExpected = 0;
     676          44 :     for (const auto &kv : oMap)
     677             :     {
     678          24 :         if (kv.first != iExpected)
     679           0 :             return false;
     680          24 :         ++iExpected;
     681             :     }
     682          20 :     return iExpected == nExpectedBandCount;
     683             : }
     684             : 
     685             : /************************************************************************/
     686             : /*                        LocalScaleOffsetInit()                        */
     687             : /************************************************************************/
     688             : 
     689             : /** Init function for 'LocalScaleOffset' builtin function. */
     690             : static CPLErr
     691           9 : LocalScaleOffsetInit(const char * /*pszFuncName*/, void * /*pUserData*/,
     692             :                      CSLConstList papszFunctionArgs, int nInBands,
     693             :                      GDALDataType eInDT, double *padfInNoData, int *pnOutBands,
     694             :                      GDALDataType *peOutDT, double **ppadfOutNoData,
     695             :                      const char *pszVRTPath, VRTPDWorkingDataPtr *ppWorkingData)
     696             : {
     697           9 :     CPLAssert(eInDT == GDT_Float64);
     698             : 
     699           9 :     const bool bIsFinalStep = *pnOutBands != 0;
     700           9 :     *peOutDT = eInDT;
     701           9 :     *ppWorkingData = nullptr;
     702             : 
     703           9 :     if (!bIsFinalStep)
     704             :     {
     705           0 :         *pnOutBands = nInBands;
     706             :     }
     707             : 
     708          18 :     auto data = std::make_unique<LocalScaleOffsetData>();
     709             : 
     710           9 :     bool bNodataSpecified = false;
     711           9 :     double dfNoData = std::numeric_limits<double>::quiet_NaN();
     712             : 
     713           9 :     bool bGainNodataSpecified = false;
     714           9 :     bool bOffsetNodataSpecified = false;
     715             : 
     716          18 :     std::map<int, std::string> oGainDatasetNameMap;
     717          18 :     std::map<int, int> oGainDatasetBandMap;
     718             : 
     719          18 :     std::map<int, std::string> oOffsetDatasetNameMap;
     720          18 :     std::map<int, int> oOffsetDatasetBandMap;
     721             : 
     722           9 :     bool bRelativeToVRT = false;
     723             : 
     724          68 :     for (const auto &[pszKey, pszValue] :
     725          73 :          cpl::IterateNameValue(papszFunctionArgs))
     726             :     {
     727          36 :         if (EQUAL(pszKey, "relativeToVRT"))
     728             :         {
     729           0 :             bRelativeToVRT = CPLTestBool(pszValue);
     730             :         }
     731          36 :         else if (EQUAL(pszKey, "nodata"))
     732             :         {
     733           0 :             bNodataSpecified = true;
     734           0 :             dfNoData = CPLAtof(pszValue);
     735             :         }
     736          36 :         else if (EQUAL(pszKey, "gain_nodata"))
     737             :         {
     738           0 :             bGainNodataSpecified = true;
     739           0 :             data->m_dfGainNodata = CPLAtof(pszValue);
     740             :         }
     741          36 :         else if (EQUAL(pszKey, "offset_nodata"))
     742             :         {
     743           0 :             bOffsetNodataSpecified = true;
     744           0 :             data->m_dfOffsetNodata = CPLAtof(pszValue);
     745             :         }
     746          36 :         else if (STARTS_WITH_CI(pszKey, "gain_dataset_filename_"))
     747             :         {
     748          10 :             const int nBand = atoi(pszKey + strlen("gain_dataset_filename_"));
     749          10 :             if (nBand <= 0 || nBand > nInBands)
     750             :             {
     751           1 :                 CPLError(CE_Failure, CPLE_AppDefined,
     752             :                          "Invalid band in argument '%s'", pszKey);
     753           4 :                 return CE_Failure;
     754             :             }
     755           9 :             oGainDatasetNameMap[nBand - 1] = pszValue;
     756             :         }
     757          26 :         else if (STARTS_WITH_CI(pszKey, "gain_dataset_band_"))
     758             :         {
     759           9 :             const int nBand = atoi(pszKey + strlen("gain_dataset_band_"));
     760           9 :             if (nBand <= 0 || nBand > nInBands)
     761             :             {
     762           1 :                 CPLError(CE_Failure, CPLE_AppDefined,
     763             :                          "Invalid band in argument '%s'", pszKey);
     764           1 :                 return CE_Failure;
     765             :             }
     766           8 :             oGainDatasetBandMap[nBand - 1] = atoi(pszValue);
     767             :         }
     768          17 :         else if (STARTS_WITH_CI(pszKey, "offset_dataset_filename_"))
     769             :         {
     770           8 :             const int nBand = atoi(pszKey + strlen("offset_dataset_filename_"));
     771           8 :             if (nBand <= 0 || nBand > nInBands)
     772             :             {
     773           1 :                 CPLError(CE_Failure, CPLE_AppDefined,
     774             :                          "Invalid band in argument '%s'", pszKey);
     775           1 :                 return CE_Failure;
     776             :             }
     777           7 :             oOffsetDatasetNameMap[nBand - 1] = pszValue;
     778             :         }
     779           9 :         else if (STARTS_WITH_CI(pszKey, "offset_dataset_band_"))
     780             :         {
     781           7 :             const int nBand = atoi(pszKey + strlen("offset_dataset_band_"));
     782           7 :             if (nBand <= 0 || nBand > nInBands)
     783             :             {
     784           1 :                 CPLError(CE_Failure, CPLE_AppDefined,
     785             :                          "Invalid band in argument '%s'", pszKey);
     786           1 :                 return CE_Failure;
     787             :             }
     788           6 :             oOffsetDatasetBandMap[nBand - 1] = atoi(pszValue);
     789             :         }
     790           2 :         else if (EQUAL(pszKey, "min"))
     791             :         {
     792           1 :             data->m_dfClampMin = CPLAtof(pszValue);
     793             :         }
     794           1 :         else if (EQUAL(pszKey, "max"))
     795             :         {
     796           1 :             data->m_dfClampMax = CPLAtof(pszValue);
     797             :         }
     798             :         else
     799             :         {
     800           0 :             CPLError(CE_Warning, CPLE_AppDefined,
     801             :                      "Unrecognized argument name %s. Ignored", pszKey);
     802             :         }
     803             :     }
     804             : 
     805           5 :     if (!CheckAllBands(oGainDatasetNameMap, nInBands))
     806             :     {
     807           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     808             :                  "Missing gain_dataset_filename_XX element(s)");
     809           0 :         return CE_Failure;
     810             :     }
     811           5 :     if (!CheckAllBands(oGainDatasetBandMap, nInBands))
     812             :     {
     813           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     814             :                  "Missing gain_dataset_band_XX element(s)");
     815           0 :         return CE_Failure;
     816             :     }
     817           5 :     if (!CheckAllBands(oOffsetDatasetNameMap, nInBands))
     818             :     {
     819           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     820             :                  "Missing offset_dataset_filename_XX element(s)");
     821           0 :         return CE_Failure;
     822             :     }
     823           5 :     if (!CheckAllBands(oOffsetDatasetBandMap, nInBands))
     824             :     {
     825           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     826             :                  "Missing offset_dataset_band_XX element(s)");
     827           0 :         return CE_Failure;
     828             :     }
     829             : 
     830           5 :     data->m_oGainBands.resize(nInBands);
     831           5 :     data->m_oOffsetBands.resize(nInBands);
     832             : 
     833           5 :     constexpr int IDX_GAIN = 0;
     834           5 :     constexpr int IDX_OFFSET = 1;
     835           9 :     for (int i : {IDX_GAIN, IDX_OFFSET})
     836             :     {
     837           7 :         const auto &oMapNames =
     838             :             (i == IDX_GAIN) ? oGainDatasetNameMap : oOffsetDatasetNameMap;
     839           7 :         const auto &oMapBands =
     840             :             (i == IDX_GAIN) ? oGainDatasetBandMap : oOffsetDatasetBandMap;
     841          13 :         for (const auto &kv : oMapNames)
     842             :         {
     843           9 :             const int nInBandIdx = kv.first;
     844             :             const auto osFilename = VRTDataset::BuildSourceFilename(
     845           9 :                 kv.second.c_str(), pszVRTPath, bRelativeToVRT);
     846           9 :             auto oIter = data->m_oDatasetMap.find(osFilename);
     847           9 :             if (oIter == data->m_oDatasetMap.end())
     848             :             {
     849             :                 auto poDS = std::unique_ptr<GDALDataset>(GDALDataset::Open(
     850             :                     osFilename.c_str(), GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR,
     851           7 :                     nullptr, nullptr, nullptr));
     852           7 :                 if (!poDS)
     853           1 :                     return CE_Failure;
     854             :                 double adfAuxGT[6];
     855           6 :                 if (poDS->GetGeoTransform(adfAuxGT) != CE_None)
     856             :                 {
     857           1 :                     CPLError(CE_Failure, CPLE_AppDefined,
     858             :                              "%s lacks a geotransform", osFilename.c_str());
     859           1 :                     return CE_Failure;
     860             :                 }
     861           5 :                 oIter = data->m_oDatasetMap
     862           5 :                             .insert(std::pair(osFilename, std::move(poDS)))
     863             :                             .first;
     864             :             }
     865           7 :             auto poDS = oIter->second.get();
     866           7 :             const auto oIterBand = oMapBands.find(nInBandIdx);
     867           7 :             CPLAssert(oIterBand != oMapBands.end());
     868           7 :             const int nAuxBand = oIterBand->second;
     869           7 :             if (nAuxBand <= 0 || nAuxBand > poDS->GetRasterCount())
     870             :             {
     871           1 :                 CPLError(CE_Failure, CPLE_AppDefined,
     872             :                          "Invalid band number (%d) for a %s dataset", nAuxBand,
     873             :                          (i == IDX_GAIN) ? "gain" : "offset");
     874           1 :                 return CE_Failure;
     875             :             }
     876           6 :             auto poAuxBand = poDS->GetRasterBand(nAuxBand);
     877           6 :             int bAuxBandHasNoData = false;
     878             :             const double dfAuxNoData =
     879           6 :                 poAuxBand->GetNoDataValue(&bAuxBandHasNoData);
     880           6 :             if (i == IDX_GAIN)
     881             :             {
     882           3 :                 data->m_oGainBands[nInBandIdx] = poAuxBand;
     883           3 :                 if (!bGainNodataSpecified && bAuxBandHasNoData)
     884           2 :                     data->m_dfGainNodata = dfAuxNoData;
     885             :             }
     886             :             else
     887             :             {
     888           3 :                 data->m_oOffsetBands[nInBandIdx] = poAuxBand;
     889           3 :                 if (!bOffsetNodataSpecified && bAuxBandHasNoData)
     890           2 :                     data->m_dfOffsetNodata = dfAuxNoData;
     891             :             }
     892             :         }
     893             :     }
     894             : 
     895           2 :     SetOutputValuesForInNoDataAndOutNoData(
     896             :         nInBands, padfInNoData, pnOutBands, ppadfOutNoData, bNodataSpecified,
     897             :         dfNoData, bNodataSpecified, dfNoData, bIsFinalStep);
     898             : 
     899           2 :     *ppWorkingData = data.release();
     900           2 :     return CE_None;
     901             : }
     902             : 
     903             : /************************************************************************/
     904             : /*                       LocalScaleOffsetFree()                         */
     905             : /************************************************************************/
     906             : 
     907             : /** Free function for 'LocalScaleOffset' builtin function. */
     908           2 : static void LocalScaleOffsetFree(const char * /*pszFuncName*/,
     909             :                                  void * /*pUserData*/,
     910             :                                  VRTPDWorkingDataPtr pWorkingData)
     911             : {
     912           2 :     LocalScaleOffsetData *data =
     913             :         static_cast<LocalScaleOffsetData *>(pWorkingData);
     914           2 :     CPLAssert(data->m_osSignature == LocalScaleOffsetData::EXPECTED_SIGNATURE);
     915           2 :     CPL_IGNORE_RET_VAL(data->m_osSignature);
     916           2 :     delete data;
     917           2 : }
     918             : 
     919             : /************************************************************************/
     920             : /*                          LoadAuxData()                               */
     921             : /************************************************************************/
     922             : 
     923             : // Load auxiliary corresponding offset, gain or trimming data.
     924           7 : static bool LoadAuxData(double dfULX, double dfULY, double dfLRX, double dfLRY,
     925             :                         size_t nElts, int nBufXSize, int nBufYSize,
     926             :                         const char *pszAuxType, GDALRasterBand *poAuxBand,
     927             :                         std::vector<VRTProcessedDataset::NoInitByte> &abyBuffer)
     928             : {
     929             :     double adfAuxGT[6];
     930             :     double adfAuxInvGT[6];
     931             : 
     932             :     // Compute pixel/line coordinates from the georeferenced extent
     933          14 :     CPL_IGNORE_RET_VAL(poAuxBand->GetDataset()->GetGeoTransform(
     934           7 :         adfAuxGT));  // return code already tested
     935           7 :     CPL_IGNORE_RET_VAL(GDALInvGeoTransform(adfAuxGT, adfAuxInvGT));
     936           7 :     const double dfULPixel =
     937           7 :         adfAuxInvGT[0] + adfAuxInvGT[1] * dfULX + adfAuxInvGT[2] * dfULY;
     938           7 :     const double dfULLine =
     939           7 :         adfAuxInvGT[3] + adfAuxInvGT[4] * dfULX + adfAuxInvGT[5] * dfULY;
     940           7 :     const double dfLRPixel =
     941           7 :         adfAuxInvGT[0] + adfAuxInvGT[1] * dfLRX + adfAuxInvGT[2] * dfLRY;
     942           7 :     const double dfLRLine =
     943           7 :         adfAuxInvGT[3] + adfAuxInvGT[4] * dfLRX + adfAuxInvGT[5] * dfLRY;
     944           7 :     if (dfULPixel >= dfLRPixel || dfULLine >= dfLRLine)
     945             :     {
     946           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     947             :                  "Unexpected computed %s pixel/line", pszAuxType);
     948           0 :         return false;
     949             :     }
     950           7 :     if (dfULPixel < -1 || dfULLine < -1)
     951             :     {
     952           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     953             :                  "Unexpected computed %s upper left (pixel,line)=(%f,%f)",
     954             :                  pszAuxType, dfULPixel, dfULLine);
     955           0 :         return false;
     956             :     }
     957          14 :     if (dfLRPixel > poAuxBand->GetXSize() + 1 ||
     958           7 :         dfLRLine > poAuxBand->GetYSize() + 1)
     959             :     {
     960           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     961             :                  "Unexpected computed %s lower right (pixel,line)=(%f,%f)",
     962             :                  pszAuxType, dfLRPixel, dfLRLine);
     963           0 :         return false;
     964             :     }
     965             : 
     966           7 :     const int nAuxXOff = std::max(0, static_cast<int>(std::round(dfULPixel)));
     967           7 :     const int nAuxYOff = std::max(0, static_cast<int>(std::round(dfULLine)));
     968          21 :     const int nAuxX2Off = std::min(poAuxBand->GetXSize(),
     969           7 :                                    static_cast<int>(std::round(dfLRPixel)));
     970             :     const int nAuxY2Off =
     971           7 :         std::min(poAuxBand->GetYSize(), static_cast<int>(std::round(dfLRLine)));
     972             : 
     973             :     try
     974             :     {
     975           7 :         abyBuffer.resize(nElts * sizeof(float));
     976             :     }
     977           0 :     catch (const std::bad_alloc &)
     978             :     {
     979           0 :         CPLError(CE_Failure, CPLE_OutOfMemory,
     980             :                  "Out of memory allocating working buffer");
     981           0 :         return false;
     982             :     }
     983             :     GDALRasterIOExtraArg sExtraArg;
     984           7 :     INIT_RASTERIO_EXTRA_ARG(sExtraArg);
     985           7 :     sExtraArg.bFloatingPointWindowValidity = true;
     986           7 :     CPL_IGNORE_RET_VAL(sExtraArg.eResampleAlg);
     987           7 :     sExtraArg.eResampleAlg = GRIORA_Bilinear;
     988           7 :     sExtraArg.dfXOff = std::max(0.0, dfULPixel);
     989           7 :     sExtraArg.dfYOff = std::max(0.0, dfULLine);
     990           7 :     sExtraArg.dfXSize = std::min<double>(poAuxBand->GetXSize(), dfLRPixel) -
     991           7 :                         std::max(0.0, dfULPixel);
     992           7 :     sExtraArg.dfYSize = std::min<double>(poAuxBand->GetYSize(), dfLRLine) -
     993           7 :                         std::max(0.0, dfULLine);
     994           7 :     return (poAuxBand->RasterIO(
     995           7 :                 GF_Read, nAuxXOff, nAuxYOff, std::max(1, nAuxX2Off - nAuxXOff),
     996           7 :                 std::max(1, nAuxY2Off - nAuxYOff), abyBuffer.data(), nBufXSize,
     997           7 :                 nBufYSize, GDT_Float32, 0, 0, &sExtraArg) == CE_None);
     998             : }
     999             : 
    1000             : /************************************************************************/
    1001             : /*                      LocalScaleOffsetProcess()                       */
    1002             : /************************************************************************/
    1003             : 
    1004             : /** Processing function for 'LocalScaleOffset' builtin function. */
    1005           2 : static CPLErr LocalScaleOffsetProcess(
    1006             :     const char * /*pszFuncName*/, void * /*pUserData*/,
    1007             :     VRTPDWorkingDataPtr pWorkingData, CSLConstList /* papszFunctionArgs*/,
    1008             :     int nBufXSize, int nBufYSize, const void *pInBuffer, size_t nInBufferSize,
    1009             :     GDALDataType eInDT, int nInBands, const double *CPL_RESTRICT padfInNoData,
    1010             :     void *pOutBuffer, size_t nOutBufferSize, GDALDataType eOutDT, int nOutBands,
    1011             :     const double *CPL_RESTRICT padfOutNoData, double dfSrcXOff,
    1012             :     double dfSrcYOff, double dfSrcXSize, double dfSrcYSize,
    1013             :     const double adfSrcGT[], const char * /* pszVRTPath */,
    1014             :     CSLConstList /*papszExtra*/)
    1015             : {
    1016           2 :     const size_t nElts = static_cast<size_t>(nBufXSize) * nBufYSize;
    1017             : 
    1018           2 :     CPL_IGNORE_RET_VAL(eInDT);
    1019           2 :     CPLAssert(eInDT == GDT_Float64);
    1020           2 :     CPL_IGNORE_RET_VAL(eOutDT);
    1021           2 :     CPLAssert(eOutDT == GDT_Float64);
    1022           2 :     CPL_IGNORE_RET_VAL(nInBufferSize);
    1023           2 :     CPLAssert(nInBufferSize == nElts * nInBands * sizeof(double));
    1024           2 :     CPL_IGNORE_RET_VAL(nOutBufferSize);
    1025           2 :     CPLAssert(nOutBufferSize == nElts * nOutBands * sizeof(double));
    1026           2 :     CPLAssert(nInBands == nOutBands);
    1027           2 :     CPL_IGNORE_RET_VAL(nOutBands);
    1028             : 
    1029           2 :     LocalScaleOffsetData *data =
    1030             :         static_cast<LocalScaleOffsetData *>(pWorkingData);
    1031           2 :     CPLAssert(data->m_osSignature == LocalScaleOffsetData::EXPECTED_SIGNATURE);
    1032           2 :     const double *CPL_RESTRICT padfSrc = static_cast<const double *>(pInBuffer);
    1033           2 :     double *CPL_RESTRICT padfDst = static_cast<double *>(pOutBuffer);
    1034             : 
    1035             :     // Compute georeferenced extent of input region
    1036           2 :     const double dfULX =
    1037           2 :         adfSrcGT[0] + adfSrcGT[1] * dfSrcXOff + adfSrcGT[2] * dfSrcYOff;
    1038           2 :     const double dfULY =
    1039           2 :         adfSrcGT[3] + adfSrcGT[4] * dfSrcXOff + adfSrcGT[5] * dfSrcYOff;
    1040           2 :     const double dfLRX = adfSrcGT[0] + adfSrcGT[1] * (dfSrcXOff + dfSrcXSize) +
    1041           2 :                          adfSrcGT[2] * (dfSrcYOff + dfSrcYSize);
    1042           2 :     const double dfLRY = adfSrcGT[3] + adfSrcGT[4] * (dfSrcXOff + dfSrcXSize) +
    1043           2 :                          adfSrcGT[5] * (dfSrcYOff + dfSrcYSize);
    1044             : 
    1045           2 :     auto &abyOffsetBuffer = data->m_abyGainBuffer;
    1046           2 :     auto &abyGainBuffer = data->m_abyOffsetBuffer;
    1047             : 
    1048           5 :     for (int iBand = 0; iBand < nInBands; ++iBand)
    1049             :     {
    1050           3 :         if (!LoadAuxData(dfULX, dfULY, dfLRX, dfLRY, nElts, nBufXSize,
    1051           3 :                          nBufYSize, "gain", data->m_oGainBands[iBand],
    1052           6 :                          abyGainBuffer) ||
    1053           3 :             !LoadAuxData(dfULX, dfULY, dfLRX, dfLRY, nElts, nBufXSize,
    1054           3 :                          nBufYSize, "offset", data->m_oOffsetBands[iBand],
    1055             :                          abyOffsetBuffer))
    1056             :         {
    1057           0 :             return CE_Failure;
    1058             :         }
    1059             : 
    1060           3 :         const double *CPL_RESTRICT padfSrcThisBand = padfSrc + iBand;
    1061           3 :         double *CPL_RESTRICT padfDstThisBand = padfDst + iBand;
    1062             :         const float *pafGain =
    1063           3 :             reinterpret_cast<const float *>(abyGainBuffer.data());
    1064             :         const float *pafOffset =
    1065           3 :             reinterpret_cast<const float *>(abyOffsetBuffer.data());
    1066           3 :         const double dfSrcNodata = padfInNoData[iBand];
    1067           3 :         const double dfDstNodata = padfOutNoData[iBand];
    1068           3 :         const double dfGainNodata = data->m_dfGainNodata;
    1069           3 :         const double dfOffsetNodata = data->m_dfOffsetNodata;
    1070           3 :         const double dfClampMin = data->m_dfClampMin;
    1071           3 :         const double dfClampMax = data->m_dfClampMax;
    1072          27 :         for (size_t i = 0; i < nElts; ++i)
    1073             :         {
    1074          24 :             const double dfSrcVal = *padfSrcThisBand;
    1075             :             // written this way to work with a NaN value
    1076          24 :             if (!(dfSrcVal != dfSrcNodata))
    1077             :             {
    1078           2 :                 *padfDstThisBand = dfDstNodata;
    1079             :             }
    1080             :             else
    1081             :             {
    1082          22 :                 const double dfGain = pafGain[i];
    1083          22 :                 const double dfOffset = pafOffset[i];
    1084          22 :                 if (!(dfGain != dfGainNodata) || !(dfOffset != dfOffsetNodata))
    1085             :                 {
    1086           4 :                     *padfDstThisBand = dfDstNodata;
    1087             :                 }
    1088             :                 else
    1089             :                 {
    1090          18 :                     double dfUnscaled = dfSrcVal * dfGain - dfOffset;
    1091          18 :                     if (dfUnscaled < dfClampMin)
    1092           2 :                         dfUnscaled = dfClampMin;
    1093          18 :                     if (dfUnscaled > dfClampMax)
    1094           1 :                         dfUnscaled = dfClampMax;
    1095             : 
    1096          18 :                     *padfDstThisBand = dfUnscaled;
    1097             :                 }
    1098             :             }
    1099          24 :             padfSrcThisBand += nInBands;
    1100          24 :             padfDstThisBand += nInBands;
    1101             :         }
    1102             :     }
    1103             : 
    1104           2 :     return CE_None;
    1105             : }
    1106             : 
    1107             : /************************************************************************/
    1108             : /*                           TrimmingData                               */
    1109             : /************************************************************************/
    1110             : 
    1111             : namespace
    1112             : {
    1113             : /** Working structure for 'Trimming' builtin function. */
    1114             : struct TrimmingData
    1115             : {
    1116             :     static constexpr const char *const EXPECTED_SIGNATURE = "Trimming";
    1117             :     //! Signature (to make sure callback functions are called with the right argument)
    1118             :     const std::string m_osSignature = EXPECTED_SIGNATURE;
    1119             : 
    1120             :     //! Nodata value for trimming dataset
    1121             :     double m_dfTrimmingNodata = std::numeric_limits<double>::quiet_NaN();
    1122             : 
    1123             :     //! Maximum saturating RGB output value.
    1124             :     double m_dfTopRGB = 0;
    1125             : 
    1126             :     //! Maximum threshold beyond which we give up saturation
    1127             :     double m_dfToneCeil = 0;
    1128             : 
    1129             :     //! Margin to allow for dynamics in brighest areas (in [0,1] range)
    1130             :     double m_dfTopMargin = 0;
    1131             : 
    1132             :     //! Index (zero-based) of input/output red band.
    1133             :     int m_nRedBand = 1 - 1;
    1134             : 
    1135             :     //! Index (zero-based) of input/output green band.
    1136             :     int m_nGreenBand = 2 - 1;
    1137             : 
    1138             :     //! Index (zero-based) of input/output blue band.
    1139             :     int m_nBlueBand = 3 - 1;
    1140             : 
    1141             :     //! Trimming dataset
    1142             :     std::unique_ptr<GDALDataset> m_poTrimmingDS{};
    1143             : 
    1144             :     //! Trimming raster band.
    1145             :     GDALRasterBand *m_poTrimmingBand = nullptr;
    1146             : 
    1147             :     //! Working buffer that contain trimming values.
    1148             :     std::vector<VRTProcessedDataset::NoInitByte> m_abyTrimmingBuffer{};
    1149             : };
    1150             : }  // namespace
    1151             : 
    1152             : /************************************************************************/
    1153             : /*                           TrimmingInit()                             */
    1154             : /************************************************************************/
    1155             : 
    1156             : /** Init function for 'Trimming' builtin function. */
    1157          12 : static CPLErr TrimmingInit(const char * /*pszFuncName*/, void * /*pUserData*/,
    1158             :                            CSLConstList papszFunctionArgs, int nInBands,
    1159             :                            GDALDataType eInDT, double *padfInNoData,
    1160             :                            int *pnOutBands, GDALDataType *peOutDT,
    1161             :                            double **ppadfOutNoData, const char *pszVRTPath,
    1162             :                            VRTPDWorkingDataPtr *ppWorkingData)
    1163             : {
    1164          12 :     CPLAssert(eInDT == GDT_Float64);
    1165             : 
    1166          12 :     const bool bIsFinalStep = *pnOutBands != 0;
    1167          12 :     *peOutDT = eInDT;
    1168          12 :     *ppWorkingData = nullptr;
    1169             : 
    1170          12 :     if (!bIsFinalStep)
    1171             :     {
    1172           0 :         *pnOutBands = nInBands;
    1173             :     }
    1174             : 
    1175          24 :     auto data = std::make_unique<TrimmingData>();
    1176             : 
    1177          12 :     bool bNodataSpecified = false;
    1178          12 :     double dfNoData = std::numeric_limits<double>::quiet_NaN();
    1179          24 :     std::string osTrimmingFilename;
    1180          12 :     bool bTrimmingNodataSpecified = false;
    1181          12 :     bool bRelativeToVRT = false;
    1182             : 
    1183          84 :     for (const auto &[pszKey, pszValue] :
    1184          90 :          cpl::IterateNameValue(papszFunctionArgs))
    1185             :     {
    1186          45 :         if (EQUAL(pszKey, "relativeToVRT"))
    1187             :         {
    1188           0 :             bRelativeToVRT = CPLTestBool(pszValue);
    1189             :         }
    1190          45 :         else if (EQUAL(pszKey, "nodata"))
    1191             :         {
    1192           0 :             bNodataSpecified = true;
    1193           0 :             dfNoData = CPLAtof(pszValue);
    1194             :         }
    1195          45 :         else if (EQUAL(pszKey, "trimming_nodata"))
    1196             :         {
    1197           0 :             bTrimmingNodataSpecified = true;
    1198           0 :             data->m_dfTrimmingNodata = CPLAtof(pszValue);
    1199             :         }
    1200          45 :         else if (EQUAL(pszKey, "trimming_dataset_filename"))
    1201             :         {
    1202          12 :             osTrimmingFilename = pszValue;
    1203             :         }
    1204          33 :         else if (EQUAL(pszKey, "red_band"))
    1205             :         {
    1206           5 :             const int nBand = atoi(pszValue) - 1;
    1207           5 :             if (nBand < 0 || nBand >= nInBands)
    1208             :             {
    1209           2 :                 CPLError(CE_Failure, CPLE_AppDefined,
    1210             :                          "Invalid band in argument '%s'", pszKey);
    1211           6 :                 return CE_Failure;
    1212             :             }
    1213           3 :             data->m_nRedBand = nBand;
    1214             :         }
    1215          28 :         else if (EQUAL(pszKey, "green_band"))
    1216             :         {
    1217           5 :             const int nBand = atoi(pszValue) - 1;
    1218           5 :             if (nBand < 0 || nBand >= nInBands)
    1219             :             {
    1220           2 :                 CPLError(CE_Failure, CPLE_AppDefined,
    1221             :                          "Invalid band in argument '%s'", pszKey);
    1222           2 :                 return CE_Failure;
    1223             :             }
    1224           3 :             data->m_nGreenBand = nBand;
    1225             :         }
    1226          23 :         else if (EQUAL(pszKey, "blue_band"))
    1227             :         {
    1228           5 :             const int nBand = atoi(pszValue) - 1;
    1229           5 :             if (nBand < 0 || nBand >= nInBands)
    1230             :             {
    1231           2 :                 CPLError(CE_Failure, CPLE_AppDefined,
    1232             :                          "Invalid band in argument '%s'", pszKey);
    1233           2 :                 return CE_Failure;
    1234             :             }
    1235           3 :             data->m_nBlueBand = nBand;
    1236             :         }
    1237          18 :         else if (EQUAL(pszKey, "top_rgb"))
    1238             :         {
    1239           6 :             data->m_dfTopRGB = CPLAtof(pszValue);
    1240             :         }
    1241          12 :         else if (EQUAL(pszKey, "tone_ceil"))
    1242             :         {
    1243           6 :             data->m_dfToneCeil = CPLAtof(pszValue);
    1244             :         }
    1245           6 :         else if (EQUAL(pszKey, "top_margin"))
    1246             :         {
    1247           6 :             data->m_dfTopMargin = CPLAtof(pszValue);
    1248             :         }
    1249             :         else
    1250             :         {
    1251           0 :             CPLError(CE_Warning, CPLE_AppDefined,
    1252             :                      "Unrecognized argument name %s. Ignored", pszKey);
    1253             :         }
    1254             :     }
    1255             : 
    1256           6 :     if (data->m_nRedBand == data->m_nGreenBand ||
    1257          10 :         data->m_nRedBand == data->m_nBlueBand ||
    1258           4 :         data->m_nGreenBand == data->m_nBlueBand)
    1259             :     {
    1260           3 :         CPLError(
    1261             :             CE_Failure, CPLE_NotSupported,
    1262             :             "red_band, green_band and blue_band must have distinct values");
    1263           3 :         return CE_Failure;
    1264             :     }
    1265             : 
    1266             :     const auto osFilename = VRTDataset::BuildSourceFilename(
    1267           6 :         osTrimmingFilename.c_str(), pszVRTPath, bRelativeToVRT);
    1268           3 :     data->m_poTrimmingDS.reset(GDALDataset::Open(
    1269             :         osFilename.c_str(), GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, nullptr,
    1270             :         nullptr, nullptr));
    1271           3 :     if (!data->m_poTrimmingDS)
    1272           1 :         return CE_Failure;
    1273           2 :     if (data->m_poTrimmingDS->GetRasterCount() != 1)
    1274             :     {
    1275           1 :         CPLError(CE_Failure, CPLE_NotSupported,
    1276             :                  "Trimming dataset should have a single band");
    1277           1 :         return CE_Failure;
    1278             :     }
    1279           1 :     data->m_poTrimmingBand = data->m_poTrimmingDS->GetRasterBand(1);
    1280             : 
    1281             :     double adfAuxGT[6];
    1282           1 :     if (data->m_poTrimmingDS->GetGeoTransform(adfAuxGT) != CE_None)
    1283             :     {
    1284           0 :         CPLError(CE_Failure, CPLE_AppDefined, "%s lacks a geotransform",
    1285             :                  osFilename.c_str());
    1286           0 :         return CE_Failure;
    1287             :     }
    1288           1 :     int bAuxBandHasNoData = false;
    1289             :     const double dfAuxNoData =
    1290           1 :         data->m_poTrimmingBand->GetNoDataValue(&bAuxBandHasNoData);
    1291           1 :     if (!bTrimmingNodataSpecified && bAuxBandHasNoData)
    1292           0 :         data->m_dfTrimmingNodata = dfAuxNoData;
    1293             : 
    1294           1 :     SetOutputValuesForInNoDataAndOutNoData(
    1295             :         nInBands, padfInNoData, pnOutBands, ppadfOutNoData, bNodataSpecified,
    1296             :         dfNoData, bNodataSpecified, dfNoData, bIsFinalStep);
    1297             : 
    1298           1 :     *ppWorkingData = data.release();
    1299           1 :     return CE_None;
    1300             : }
    1301             : 
    1302             : /************************************************************************/
    1303             : /*                           TrimmingFree()                             */
    1304             : /************************************************************************/
    1305             : 
    1306             : /** Free function for 'Trimming' builtin function. */
    1307           1 : static void TrimmingFree(const char * /*pszFuncName*/, void * /*pUserData*/,
    1308             :                          VRTPDWorkingDataPtr pWorkingData)
    1309             : {
    1310           1 :     TrimmingData *data = static_cast<TrimmingData *>(pWorkingData);
    1311           1 :     CPLAssert(data->m_osSignature == TrimmingData::EXPECTED_SIGNATURE);
    1312           1 :     CPL_IGNORE_RET_VAL(data->m_osSignature);
    1313           1 :     delete data;
    1314           1 : }
    1315             : 
    1316             : /************************************************************************/
    1317             : /*                         TrimmingProcess()                            */
    1318             : /************************************************************************/
    1319             : 
    1320             : /** Processing function for 'Trimming' builtin function. */
    1321           1 : static CPLErr TrimmingProcess(
    1322             :     const char * /*pszFuncName*/, void * /*pUserData*/,
    1323             :     VRTPDWorkingDataPtr pWorkingData, CSLConstList /* papszFunctionArgs*/,
    1324             :     int nBufXSize, int nBufYSize, const void *pInBuffer, size_t nInBufferSize,
    1325             :     GDALDataType eInDT, int nInBands, const double *CPL_RESTRICT padfInNoData,
    1326             :     void *pOutBuffer, size_t nOutBufferSize, GDALDataType eOutDT, int nOutBands,
    1327             :     const double *CPL_RESTRICT padfOutNoData, double dfSrcXOff,
    1328             :     double dfSrcYOff, double dfSrcXSize, double dfSrcYSize,
    1329             :     const double adfSrcGT[], const char * /* pszVRTPath */,
    1330             :     CSLConstList /*papszExtra*/)
    1331             : {
    1332           1 :     const size_t nElts = static_cast<size_t>(nBufXSize) * nBufYSize;
    1333             : 
    1334           1 :     CPL_IGNORE_RET_VAL(eInDT);
    1335           1 :     CPLAssert(eInDT == GDT_Float64);
    1336           1 :     CPL_IGNORE_RET_VAL(eOutDT);
    1337           1 :     CPLAssert(eOutDT == GDT_Float64);
    1338           1 :     CPL_IGNORE_RET_VAL(nInBufferSize);
    1339           1 :     CPLAssert(nInBufferSize == nElts * nInBands * sizeof(double));
    1340           1 :     CPL_IGNORE_RET_VAL(nOutBufferSize);
    1341           1 :     CPLAssert(nOutBufferSize == nElts * nOutBands * sizeof(double));
    1342           1 :     CPLAssert(nInBands == nOutBands);
    1343           1 :     CPL_IGNORE_RET_VAL(nOutBands);
    1344             : 
    1345           1 :     TrimmingData *data = static_cast<TrimmingData *>(pWorkingData);
    1346           1 :     CPLAssert(data->m_osSignature == TrimmingData::EXPECTED_SIGNATURE);
    1347           1 :     const double *CPL_RESTRICT padfSrc = static_cast<const double *>(pInBuffer);
    1348           1 :     double *CPL_RESTRICT padfDst = static_cast<double *>(pOutBuffer);
    1349             : 
    1350             :     // Compute georeferenced extent of input region
    1351           1 :     const double dfULX =
    1352           1 :         adfSrcGT[0] + adfSrcGT[1] * dfSrcXOff + adfSrcGT[2] * dfSrcYOff;
    1353           1 :     const double dfULY =
    1354           1 :         adfSrcGT[3] + adfSrcGT[4] * dfSrcXOff + adfSrcGT[5] * dfSrcYOff;
    1355           1 :     const double dfLRX = adfSrcGT[0] + adfSrcGT[1] * (dfSrcXOff + dfSrcXSize) +
    1356           1 :                          adfSrcGT[2] * (dfSrcYOff + dfSrcYSize);
    1357           1 :     const double dfLRY = adfSrcGT[3] + adfSrcGT[4] * (dfSrcXOff + dfSrcXSize) +
    1358           1 :                          adfSrcGT[5] * (dfSrcYOff + dfSrcYSize);
    1359             : 
    1360           1 :     if (!LoadAuxData(dfULX, dfULY, dfLRX, dfLRY, nElts, nBufXSize, nBufYSize,
    1361             :                      "trimming", data->m_poTrimmingBand,
    1362           1 :                      data->m_abyTrimmingBuffer))
    1363             :     {
    1364           0 :         return CE_Failure;
    1365             :     }
    1366             : 
    1367             :     const float *pafTrimming =
    1368           1 :         reinterpret_cast<const float *>(data->m_abyTrimmingBuffer.data());
    1369           1 :     const int nRedBand = data->m_nRedBand;
    1370           1 :     const int nGreenBand = data->m_nGreenBand;
    1371           1 :     const int nBlueBand = data->m_nBlueBand;
    1372           1 :     const double dfTopMargin = data->m_dfTopMargin;
    1373           1 :     const double dfTopRGB = data->m_dfTopRGB;
    1374           1 :     const double dfToneCeil = data->m_dfToneCeil;
    1375             : #if !defined(trimming_non_optimized_version)
    1376           1 :     const double dfInvToneCeil = 1.0 / dfToneCeil;
    1377             : #endif
    1378             :     const bool bRGBBandsAreFirst =
    1379           1 :         std::max(std::max(nRedBand, nGreenBand), nBlueBand) <= 2;
    1380           1 :     const double dfNoDataTrimming = data->m_dfTrimmingNodata;
    1381           1 :     const double dfNoDataRed = padfInNoData[nRedBand];
    1382           1 :     const double dfNoDataGreen = padfInNoData[nGreenBand];
    1383           1 :     const double dfNoDataBlue = padfInNoData[nBlueBand];
    1384           7 :     for (size_t i = 0; i < nElts; ++i)
    1385             :     {
    1386             :         // Extract local saturation value from trimming image
    1387           6 :         const double dfLocalMaxRGB = pafTrimming[i];
    1388             :         const double dfReducedRGB =
    1389           6 :             std::min((1.0 - dfTopMargin) * dfTopRGB / dfLocalMaxRGB, 1.0);
    1390             : 
    1391           6 :         const double dfRed = padfSrc[nRedBand];
    1392           6 :         const double dfGreen = padfSrc[nGreenBand];
    1393           6 :         const double dfBlue = padfSrc[nBlueBand];
    1394           6 :         bool bNoDataPixel = false;
    1395           6 :         if ((dfLocalMaxRGB != dfNoDataTrimming) && (dfRed != dfNoDataRed) &&
    1396           6 :             (dfGreen != dfNoDataGreen) && (dfBlue != dfNoDataBlue))
    1397             :         {
    1398             :             // RGB bands specific process
    1399           6 :             const double dfMaxRGB = std::max(std::max(dfRed, dfGreen), dfBlue);
    1400             : #if !defined(trimming_non_optimized_version)
    1401           6 :             const double dfRedTimesToneRed = std::min(dfRed, dfToneCeil);
    1402           6 :             const double dfGreenTimesToneGreen = std::min(dfGreen, dfToneCeil);
    1403           6 :             const double dfBlueTimesToneBlue = std::min(dfBlue, dfToneCeil);
    1404             :             const double dfInvToneMaxRGB =
    1405           6 :                 std::max(dfMaxRGB * dfInvToneCeil, 1.0);
    1406           6 :             const double dfReducedRGBTimesInvToneMaxRGB =
    1407             :                 dfReducedRGB * dfInvToneMaxRGB;
    1408           6 :             padfDst[nRedBand] = std::min(
    1409           6 :                 dfRedTimesToneRed * dfReducedRGBTimesInvToneMaxRGB, dfTopRGB);
    1410           6 :             padfDst[nGreenBand] =
    1411          12 :                 std::min(dfGreenTimesToneGreen * dfReducedRGBTimesInvToneMaxRGB,
    1412           6 :                          dfTopRGB);
    1413           6 :             padfDst[nBlueBand] = std::min(
    1414           6 :                 dfBlueTimesToneBlue * dfReducedRGBTimesInvToneMaxRGB, dfTopRGB);
    1415             : #else
    1416             :             // Original formulas. Slightly less optimized than the above ones.
    1417             :             const double dfToneMaxRGB = std::min(dfToneCeil / dfMaxRGB, 1.0);
    1418             :             const double dfToneRed = std::min(dfToneCeil / dfRed, 1.0);
    1419             :             const double dfToneGreen = std::min(dfToneCeil / dfGreen, 1.0);
    1420             :             const double dfToneBlue = std::min(dfToneCeil / dfBlue, 1.0);
    1421             :             padfDst[nRedBand] = std::min(
    1422             :                 dfReducedRGB * dfRed * dfToneRed / dfToneMaxRGB, dfTopRGB);
    1423             :             padfDst[nGreenBand] = std::min(
    1424             :                 dfReducedRGB * dfGreen * dfToneGreen / dfToneMaxRGB, dfTopRGB);
    1425             :             padfDst[nBlueBand] = std::min(
    1426             :                 dfReducedRGB * dfBlue * dfToneBlue / dfToneMaxRGB, dfTopRGB);
    1427             : #endif
    1428             : 
    1429             :             // Other bands processing (NIR, ...): only apply RGB reduction factor
    1430           6 :             if (bRGBBandsAreFirst)
    1431             :             {
    1432             :                 // optimization
    1433          12 :                 for (int iBand = 3; iBand < nInBands; ++iBand)
    1434             :                 {
    1435           6 :                     if (padfSrc[iBand] != padfInNoData[iBand])
    1436             :                     {
    1437           6 :                         padfDst[iBand] = dfReducedRGB * padfSrc[iBand];
    1438             :                     }
    1439             :                     else
    1440             :                     {
    1441           0 :                         bNoDataPixel = true;
    1442           0 :                         break;
    1443             :                     }
    1444             :                 }
    1445             :             }
    1446             :             else
    1447             :             {
    1448           0 :                 for (int iBand = 0; iBand < nInBands; ++iBand)
    1449             :                 {
    1450           0 :                     if (iBand != nRedBand && iBand != nGreenBand &&
    1451           0 :                         iBand != nBlueBand)
    1452             :                     {
    1453           0 :                         if (padfSrc[iBand] != padfInNoData[iBand])
    1454             :                         {
    1455           0 :                             padfDst[iBand] = dfReducedRGB * padfSrc[iBand];
    1456             :                         }
    1457             :                         else
    1458             :                         {
    1459           0 :                             bNoDataPixel = true;
    1460           0 :                             break;
    1461             :                         }
    1462             :                     }
    1463             :                 }
    1464           6 :             }
    1465             :         }
    1466             :         else
    1467             :         {
    1468           0 :             bNoDataPixel = true;
    1469             :         }
    1470           6 :         if (bNoDataPixel)
    1471             :         {
    1472           0 :             for (int iBand = 0; iBand < nInBands; ++iBand)
    1473             :             {
    1474           0 :                 padfDst[iBand] = padfOutNoData[iBand];
    1475             :             }
    1476             :         }
    1477             : 
    1478           6 :         padfSrc += nInBands;
    1479           6 :         padfDst += nInBands;
    1480             :     }
    1481             : 
    1482           1 :     return CE_None;
    1483             : }
    1484             : 
    1485             : /************************************************************************/
    1486             : /*              GDALVRTRegisterDefaultProcessedDatasetFuncs()           */
    1487             : /************************************************************************/
    1488             : 
    1489             : /** Register builtin functions that can be used in a VRTProcessedDataset.
    1490             :  */
    1491        1218 : void GDALVRTRegisterDefaultProcessedDatasetFuncs()
    1492             : {
    1493        1218 :     GDALVRTRegisterProcessedDatasetFunc(
    1494             :         "BandAffineCombination", nullptr,
    1495             :         "<ProcessedDatasetFunctionArgumentsList>"
    1496             :         "   <Argument name='src_nodata' type='double' "
    1497             :         "description='Override input nodata value'/>"
    1498             :         "   <Argument name='dst_nodata' type='double' "
    1499             :         "description='Override output nodata value'/>"
    1500             :         "   <Argument name='replacement_nodata' "
    1501             :         "description='value to substitute to a valid computed value that "
    1502             :         "would be nodata' type='double'/>"
    1503             :         "   <Argument name='dst_intended_datatype' type='string' "
    1504             :         "description='Intented datatype of output (which might be "
    1505             :         "different than the working data type)'/>"
    1506             :         "   <Argument name='coefficients_{band}' "
    1507             :         "description='Comma-separated coefficients for combining bands. "
    1508             :         "First one is constant term' "
    1509             :         "type='double_list' required='true'/>"
    1510             :         "   <Argument name='min' description='clamp min value' type='double'/>"
    1511             :         "   <Argument name='max' description='clamp max value' type='double'/>"
    1512             :         "</ProcessedDatasetFunctionArgumentsList>",
    1513             :         GDT_Float64, nullptr, 0, nullptr, 0, BandAffineCombinationInit,
    1514             :         BandAffineCombinationFree, BandAffineCombinationProcess, nullptr);
    1515             : 
    1516        1218 :     GDALVRTRegisterProcessedDatasetFunc(
    1517             :         "LUT", nullptr,
    1518             :         "<ProcessedDatasetFunctionArgumentsList>"
    1519             :         "   <Argument name='src_nodata' type='double' "
    1520             :         "description='Override input nodata value'/>"
    1521             :         "   <Argument name='dst_nodata' type='double' "
    1522             :         "description='Override output nodata value'/>"
    1523             :         "   <Argument name='lut_{band}' "
    1524             :         "description='List of the form [src value 1]:[dest value 1],"
    1525             :         "[src value 2]:[dest value 2],...' "
    1526             :         "type='string' required='true'/>"
    1527             :         "</ProcessedDatasetFunctionArgumentsList>",
    1528             :         GDT_Float64, nullptr, 0, nullptr, 0, LUTInit, LUTFree, LUTProcess,
    1529             :         nullptr);
    1530             : 
    1531        1218 :     GDALVRTRegisterProcessedDatasetFunc(
    1532             :         "LocalScaleOffset", nullptr,
    1533             :         "<ProcessedDatasetFunctionArgumentsList>"
    1534             :         "   <Argument name='relativeToVRT' "
    1535             :         "description='Whether gain and offset filenames are relative to "
    1536             :         "the VRT' type='boolean' default='false'/>"
    1537             :         "   <Argument name='gain_dataset_filename_{band}' "
    1538             :         "description='Filename to the gain dataset' "
    1539             :         "type='string' required='true'/>"
    1540             :         "   <Argument name='gain_dataset_band_{band}' "
    1541             :         "description='Band of the gain dataset' "
    1542             :         "type='integer' required='true'/>"
    1543             :         "   <Argument name='offset_dataset_filename_{band}' "
    1544             :         "description='Filename to the offset dataset' "
    1545             :         "type='string' required='true'/>"
    1546             :         "   <Argument name='offset_dataset_band_{band}' "
    1547             :         "description='Band of the offset dataset' "
    1548             :         "type='integer' required='true'/>"
    1549             :         "   <Argument name='min' description='clamp min value' type='double'/>"
    1550             :         "   <Argument name='max' description='clamp max value' type='double'/>"
    1551             :         "   <Argument name='nodata' type='double' "
    1552             :         "description='Override dataset nodata value'/>"
    1553             :         "   <Argument name='gain_nodata' type='double' "
    1554             :         "description='Override gain dataset nodata value'/>"
    1555             :         "   <Argument name='offset_nodata' type='double' "
    1556             :         "description='Override offset dataset nodata value'/>"
    1557             :         "</ProcessedDatasetFunctionArgumentsList>",
    1558             :         GDT_Float64, nullptr, 0, nullptr, 0, LocalScaleOffsetInit,
    1559             :         LocalScaleOffsetFree, LocalScaleOffsetProcess, nullptr);
    1560             : 
    1561        1218 :     GDALVRTRegisterProcessedDatasetFunc(
    1562             :         "Trimming", nullptr,
    1563             :         "<ProcessedDatasetFunctionArgumentsList>"
    1564             :         "   <Argument name='relativeToVRT' "
    1565             :         "description='Whether trimming_dataset_filename is relative to the VRT'"
    1566             :         " type='boolean' default='false'/>"
    1567             :         "   <Argument name='trimming_dataset_filename' "
    1568             :         "description='Filename to the trimming dataset' "
    1569             :         "type='string' required='true'/>"
    1570             :         "   <Argument name='red_band' type='integer' default='1'/>"
    1571             :         "   <Argument name='green_band' type='integer' default='2'/>"
    1572             :         "   <Argument name='blue_band' type='integer' default='3'/>"
    1573             :         "   <Argument name='top_rgb' "
    1574             :         "description='Maximum saturating RGB output value' "
    1575             :         "type='double' required='true'/>"
    1576             :         "   <Argument name='tone_ceil' "
    1577             :         "description='Maximum threshold beyond which we give up saturation' "
    1578             :         "type='double' required='true'/>"
    1579             :         "   <Argument name='top_margin' "
    1580             :         "description='Margin to allow for dynamics in brighest areas "
    1581             :         "(between 0 and 1, should be close to 0)' "
    1582             :         "type='double' required='true'/>"
    1583             :         "   <Argument name='nodata' type='double' "
    1584             :         "description='Override dataset nodata value'/>"
    1585             :         "   <Argument name='trimming_nodata' type='double' "
    1586             :         "description='Override trimming dataset nodata value'/>"
    1587             :         "</ProcessedDatasetFunctionArgumentsList>",
    1588             :         GDT_Float64, nullptr, 0, nullptr, 0, TrimmingInit, TrimmingFree,
    1589             :         TrimmingProcess, nullptr);
    1590        1218 : }

Generated by: LCOV version 1.14