LCOV - code coverage report
Current view: top level - frmts/vrt - pixelfunctions.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 629 695 90.5 %
Date: 2024-05-06 22:33:47 Functions: 40 40 100.0 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL
       4             :  * Purpose:  Implementation of a set of GDALDerivedPixelFunc(s) to be used
       5             :  *           with source raster band of virtual GDAL datasets.
       6             :  * Author:   Antonio Valentino <antonio.valentino@tiscali.it>
       7             :  *
       8             :  ******************************************************************************
       9             :  * Copyright (c) 2008-2014,2022 Antonio Valentino <antonio.valentino@tiscali.it>
      10             :  *
      11             :  * Permission is hereby granted, free of charge, to any person obtaining a
      12             :  * copy of this software and associated documentation files (the "Software"),
      13             :  * to deal in the Software without restriction, including without limitation
      14             :  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
      15             :  * and/or sell copies of the Software, and to permit persons to whom the
      16             :  * Software is furnished to do so, subject to the following conditions:
      17             :  *
      18             :  * The above copyright notice and this permission notice shall be included
      19             :  * in all copies or substantial portions of the Software.
      20             :  *
      21             :  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
      22             :  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
      23             :  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
      24             :  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
      25             :  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
      26             :  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
      27             :  * DEALINGS IN THE SOFTWARE.
      28             :  *****************************************************************************/
      29             : 
      30             : #include <cmath>
      31             : #include "gdal.h"
      32             : #include "vrtdataset.h"
      33             : 
      34             : #include <limits>
      35             : 
      36             : template <typename T>
      37       55719 : inline double GetSrcVal(const void *pSource, GDALDataType eSrcType, T ii)
      38             : {
      39       55719 :     switch (eSrcType)
      40             :     {
      41           0 :         case GDT_Unknown:
      42           0 :             return 0;
      43       32950 :         case GDT_Byte:
      44       32950 :             return static_cast<const GByte *>(pSource)[ii];
      45           0 :         case GDT_Int8:
      46           0 :             return static_cast<const GInt8 *>(pSource)[ii];
      47           0 :         case GDT_UInt16:
      48           0 :             return static_cast<const GUInt16 *>(pSource)[ii];
      49           0 :         case GDT_Int16:
      50           0 :             return static_cast<const GInt16 *>(pSource)[ii];
      51           0 :         case GDT_UInt32:
      52           0 :             return static_cast<const GUInt32 *>(pSource)[ii];
      53         400 :         case GDT_Int32:
      54         400 :             return static_cast<const GInt32 *>(pSource)[ii];
      55             :         // Precision loss currently for int64/uint64
      56           0 :         case GDT_UInt64:
      57             :             return static_cast<double>(
      58           0 :                 static_cast<const uint64_t *>(pSource)[ii]);
      59           0 :         case GDT_Int64:
      60             :             return static_cast<double>(
      61           0 :                 static_cast<const int64_t *>(pSource)[ii]);
      62        8636 :         case GDT_Float32:
      63        8636 :             return static_cast<const float *>(pSource)[ii];
      64        8231 :         case GDT_Float64:
      65        8231 :             return static_cast<const double *>(pSource)[ii];
      66         360 :         case GDT_CInt16:
      67         360 :             return static_cast<const GInt16 *>(pSource)[2 * ii];
      68           0 :         case GDT_CInt32:
      69           0 :             return static_cast<const GInt32 *>(pSource)[2 * ii];
      70        4660 :         case GDT_CFloat32:
      71        4660 :             return static_cast<const float *>(pSource)[2 * ii];
      72         482 :         case GDT_CFloat64:
      73         482 :             return static_cast<const double *>(pSource)[2 * ii];
      74           0 :         case GDT_TypeCount:
      75           0 :             break;
      76             :     }
      77           0 :     return 0;
      78             : }
      79             : 
      80          57 : static CPLErr FetchDoubleArg(CSLConstList papszArgs, const char *pszName,
      81             :                              double *pdfX, double *pdfDefault = nullptr)
      82             : {
      83          57 :     const char *pszVal = CSLFetchNameValue(papszArgs, pszName);
      84             : 
      85          57 :     if (pszVal == nullptr)
      86             :     {
      87          15 :         if (pdfDefault == nullptr)
      88             :         {
      89           0 :             CPLError(CE_Failure, CPLE_AppDefined,
      90             :                      "Missing pixel function argument: %s", pszName);
      91           0 :             return CE_Failure;
      92             :         }
      93             :         else
      94             :         {
      95          15 :             *pdfX = *pdfDefault;
      96          15 :             return CE_None;
      97             :         }
      98             :     }
      99             : 
     100          42 :     char *pszEnd = nullptr;
     101          42 :     *pdfX = std::strtod(pszVal, &pszEnd);
     102          42 :     if (pszEnd == pszVal)
     103             :     {
     104           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     105             :                  "Failed to parse pixel function argument: %s", pszName);
     106           0 :         return CE_Failure;
     107             :     }
     108             : 
     109          42 :     return CE_None;
     110             : }
     111             : 
     112           5 : static CPLErr RealPixelFunc(void **papoSources, int nSources, void *pData,
     113             :                             int nXSize, int nYSize, GDALDataType eSrcType,
     114             :                             GDALDataType eBufType, int nPixelSpace,
     115             :                             int nLineSpace)
     116             : {
     117             :     /* ---- Init ---- */
     118           5 :     if (nSources != 1)
     119           1 :         return CE_Failure;
     120             : 
     121           4 :     const int nPixelSpaceSrc = GDALGetDataTypeSizeBytes(eSrcType);
     122           4 :     const size_t nLineSpaceSrc = static_cast<size_t>(nPixelSpaceSrc) * nXSize;
     123             : 
     124             :     /* ---- Set pixels ---- */
     125          56 :     for (int iLine = 0; iLine < nYSize; ++iLine)
     126             :     {
     127          52 :         GDALCopyWords(static_cast<GByte *>(papoSources[0]) +
     128          52 :                           nLineSpaceSrc * iLine,
     129             :                       eSrcType, nPixelSpaceSrc,
     130             :                       static_cast<GByte *>(pData) +
     131          52 :                           static_cast<GSpacing>(nLineSpace) * iLine,
     132             :                       eBufType, nPixelSpace, nXSize);
     133             :     }
     134             : 
     135             :     /* ---- Return success ---- */
     136           4 :     return CE_None;
     137             : }  // RealPixelFunc
     138             : 
     139           7 : static CPLErr ImagPixelFunc(void **papoSources, int nSources, void *pData,
     140             :                             int nXSize, int nYSize, GDALDataType eSrcType,
     141             :                             GDALDataType eBufType, int nPixelSpace,
     142             :                             int nLineSpace)
     143             : {
     144             :     /* ---- Init ---- */
     145           7 :     if (nSources != 1)
     146           1 :         return CE_Failure;
     147             : 
     148           6 :     if (GDALDataTypeIsComplex(eSrcType))
     149             :     {
     150           5 :         const GDALDataType eSrcBaseType = GDALGetNonComplexDataType(eSrcType);
     151           5 :         const int nPixelSpaceSrc = GDALGetDataTypeSizeBytes(eSrcType);
     152           5 :         const size_t nLineSpaceSrc =
     153           5 :             static_cast<size_t>(nPixelSpaceSrc) * nXSize;
     154             : 
     155           5 :         const void *const pImag = static_cast<GByte *>(papoSources[0]) +
     156           5 :                                   GDALGetDataTypeSizeBytes(eSrcType) / 2;
     157             : 
     158             :         /* ---- Set pixels ---- */
     159          35 :         for (int iLine = 0; iLine < nYSize; ++iLine)
     160             :         {
     161          30 :             GDALCopyWords(static_cast<const GByte *>(pImag) +
     162          30 :                               nLineSpaceSrc * iLine,
     163             :                           eSrcBaseType, nPixelSpaceSrc,
     164             :                           static_cast<GByte *>(pData) +
     165          30 :                               static_cast<GSpacing>(nLineSpace) * iLine,
     166             :                           eBufType, nPixelSpace, nXSize);
     167             :         }
     168             :     }
     169             :     else
     170             :     {
     171           1 :         const double dfImag = 0;
     172             : 
     173             :         /* ---- Set pixels ---- */
     174          21 :         for (int iLine = 0; iLine < nYSize; ++iLine)
     175             :         {
     176             :             // Always copy from the same location.
     177          20 :             GDALCopyWords(&dfImag, eSrcType, 0,
     178             :                           static_cast<GByte *>(pData) +
     179          20 :                               static_cast<GSpacing>(nLineSpace) * iLine,
     180             :                           eBufType, nPixelSpace, nXSize);
     181             :         }
     182             :     }
     183             : 
     184             :     /* ---- Return success ---- */
     185           6 :     return CE_None;
     186             : }  // ImagPixelFunc
     187             : 
     188           2 : static CPLErr ComplexPixelFunc(void **papoSources, int nSources, void *pData,
     189             :                                int nXSize, int nYSize, GDALDataType eSrcType,
     190             :                                GDALDataType eBufType, int nPixelSpace,
     191             :                                int nLineSpace)
     192             : {
     193             :     /* ---- Init ---- */
     194           2 :     if (nSources != 2)
     195           1 :         return CE_Failure;
     196             : 
     197           1 :     const void *const pReal = papoSources[0];
     198           1 :     const void *const pImag = papoSources[1];
     199             : 
     200             :     /* ---- Set pixels ---- */
     201           1 :     size_t ii = 0;
     202          21 :     for (int iLine = 0; iLine < nYSize; ++iLine)
     203             :     {
     204         420 :         for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
     205             :         {
     206             :             // Source raster pixels may be obtained with GetSrcVal macro.
     207             :             const double adfPixVal[2] = {
     208         400 :                 GetSrcVal(pReal, eSrcType, ii),  // re
     209         400 :                 GetSrcVal(pImag, eSrcType, ii)   // im
     210         400 :             };
     211             : 
     212         400 :             GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
     213             :                           static_cast<GByte *>(pData) +
     214         400 :                               static_cast<GSpacing>(nLineSpace) * iLine +
     215         400 :                               iCol * nPixelSpace,
     216             :                           eBufType, nPixelSpace, 1);
     217             :         }
     218             :     }
     219             : 
     220             :     /* ---- Return success ---- */
     221           1 :     return CE_None;
     222             : }  // ComplexPixelFunc
     223             : 
     224             : typedef enum
     225             : {
     226             :     GAT_amplitude,
     227             :     GAT_intensity,
     228             :     GAT_dB
     229             : } PolarAmplitudeType;
     230             : 
     231             : static const char pszPolarPixelFuncMetadata[] =
     232             :     "<PixelFunctionArgumentsList>"
     233             :     "   <Argument name='amplitude_type' description='Amplitude Type' "
     234             :     "type='string-select' default='AMPLITUDE'>"
     235             :     "       <Value>INTENSITY</Value>"
     236             :     "       <Value>dB</Value>"
     237             :     "       <Value>AMPLITUDE</Value>"
     238             :     "   </Argument>"
     239             :     "</PixelFunctionArgumentsList>";
     240             : 
     241           4 : static CPLErr PolarPixelFunc(void **papoSources, int nSources, void *pData,
     242             :                              int nXSize, int nYSize, GDALDataType eSrcType,
     243             :                              GDALDataType eBufType, int nPixelSpace,
     244             :                              int nLineSpace, CSLConstList papszArgs)
     245             : {
     246             :     /* ---- Init ---- */
     247           4 :     if (nSources != 2)
     248           0 :         return CE_Failure;
     249             : 
     250           4 :     const char pszName[] = "amplitude_type";
     251           4 :     const char *pszVal = CSLFetchNameValue(papszArgs, pszName);
     252           4 :     PolarAmplitudeType amplitudeType = GAT_amplitude;
     253           4 :     if (pszVal != nullptr)
     254             :     {
     255           3 :         if (strcmp(pszVal, "INTENSITY") == 0)
     256           1 :             amplitudeType = GAT_intensity;
     257           2 :         else if (strcmp(pszVal, "dB") == 0)
     258           1 :             amplitudeType = GAT_dB;
     259           1 :         else if (strcmp(pszVal, "AMPLITUDE") != 0)
     260             :         {
     261           0 :             CPLError(CE_Failure, CPLE_AppDefined,
     262             :                      "Invalid value for pixel function argument '%s': %s",
     263             :                      pszName, pszVal);
     264           0 :             return CE_Failure;
     265             :         }
     266             :     }
     267             : 
     268           4 :     const void *const pAmp = papoSources[0];
     269           4 :     const void *const pPhase = papoSources[1];
     270             : 
     271             :     /* ---- Set pixels ---- */
     272           4 :     size_t ii = 0;
     273          84 :     for (int iLine = 0; iLine < nYSize; ++iLine)
     274             :     {
     275        1680 :         for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
     276             :         {
     277             :             // Source raster pixels may be obtained with GetSrcVal macro.
     278        1600 :             double dfAmp = GetSrcVal(pAmp, eSrcType, ii);
     279        1600 :             switch (amplitudeType)
     280             :             {
     281         400 :                 case GAT_intensity:
     282             :                     // clip to zero
     283         400 :                     dfAmp = dfAmp <= 0 ? 0 : std::sqrt(dfAmp);
     284         400 :                     break;
     285         400 :                 case GAT_dB:
     286         400 :                     dfAmp = dfAmp <= 0
     287         400 :                                 ? -std::numeric_limits<double>::infinity()
     288         400 :                                 : pow(10, dfAmp / 20.);
     289         400 :                     break;
     290         800 :                 case GAT_amplitude:
     291         800 :                     break;
     292             :             }
     293        1600 :             const double dfPhase = GetSrcVal(pPhase, eSrcType, ii);
     294             :             const double adfPixVal[2] = {
     295        1600 :                 dfAmp * std::cos(dfPhase),  // re
     296        1600 :                 dfAmp * std::sin(dfPhase)   // im
     297        1600 :             };
     298             : 
     299        1600 :             GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
     300             :                           static_cast<GByte *>(pData) +
     301        1600 :                               static_cast<GSpacing>(nLineSpace) * iLine +
     302        1600 :                               iCol * nPixelSpace,
     303             :                           eBufType, nPixelSpace, 1);
     304             :         }
     305             :     }
     306             : 
     307             :     /* ---- Return success ---- */
     308           4 :     return CE_None;
     309             : }  // PolarPixelFunc
     310             : 
     311           4 : static CPLErr ModulePixelFunc(void **papoSources, int nSources, void *pData,
     312             :                               int nXSize, int nYSize, GDALDataType eSrcType,
     313             :                               GDALDataType eBufType, int nPixelSpace,
     314             :                               int nLineSpace)
     315             : {
     316             :     /* ---- Init ---- */
     317           4 :     if (nSources != 1)
     318           1 :         return CE_Failure;
     319             : 
     320           3 :     if (GDALDataTypeIsComplex(eSrcType))
     321             :     {
     322           2 :         const void *pReal = papoSources[0];
     323           2 :         const void *pImag = static_cast<GByte *>(papoSources[0]) +
     324           2 :                             GDALGetDataTypeSizeBytes(eSrcType) / 2;
     325             : 
     326             :         /* ---- Set pixels ---- */
     327           2 :         size_t ii = 0;
     328          14 :         for (int iLine = 0; iLine < nYSize; ++iLine)
     329             :         {
     330          72 :             for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
     331             :             {
     332             :                 // Source raster pixels may be obtained with GetSrcVal macro.
     333          60 :                 const double dfReal = GetSrcVal(pReal, eSrcType, ii);
     334          60 :                 const double dfImag = GetSrcVal(pImag, eSrcType, ii);
     335             : 
     336          60 :                 const double dfPixVal = sqrt(dfReal * dfReal + dfImag * dfImag);
     337             : 
     338          60 :                 GDALCopyWords(&dfPixVal, GDT_Float64, 0,
     339             :                               static_cast<GByte *>(pData) +
     340          60 :                                   static_cast<GSpacing>(nLineSpace) * iLine +
     341          60 :                                   iCol * nPixelSpace,
     342             :                               eBufType, nPixelSpace, 1);
     343             :             }
     344             :         }
     345             :     }
     346             :     else
     347             :     {
     348             :         /* ---- Set pixels ---- */
     349           1 :         size_t ii = 0;
     350          21 :         for (int iLine = 0; iLine < nYSize; ++iLine)
     351             :         {
     352         420 :             for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
     353             :             {
     354             :                 // Source raster pixels may be obtained with GetSrcVal macro.
     355             :                 const double dfPixVal =
     356         400 :                     fabs(GetSrcVal(papoSources[0], eSrcType, ii));
     357             : 
     358         400 :                 GDALCopyWords(&dfPixVal, GDT_Float64, 0,
     359             :                               static_cast<GByte *>(pData) +
     360         400 :                                   static_cast<GSpacing>(nLineSpace) * iLine +
     361         400 :                                   iCol * nPixelSpace,
     362             :                               eBufType, nPixelSpace, 1);
     363             :             }
     364             :         }
     365             :     }
     366             : 
     367             :     /* ---- Return success ---- */
     368           3 :     return CE_None;
     369             : }  // ModulePixelFunc
     370             : 
     371           5 : static CPLErr PhasePixelFunc(void **papoSources, int nSources, void *pData,
     372             :                              int nXSize, int nYSize, GDALDataType eSrcType,
     373             :                              GDALDataType eBufType, int nPixelSpace,
     374             :                              int nLineSpace)
     375             : {
     376             :     /* ---- Init ---- */
     377           5 :     if (nSources != 1)
     378           1 :         return CE_Failure;
     379             : 
     380           4 :     if (GDALDataTypeIsComplex(eSrcType))
     381             :     {
     382           2 :         const void *const pReal = papoSources[0];
     383           2 :         const void *const pImag = static_cast<GByte *>(papoSources[0]) +
     384           2 :                                   GDALGetDataTypeSizeBytes(eSrcType) / 2;
     385             : 
     386             :         /* ---- Set pixels ---- */
     387           2 :         size_t ii = 0;
     388          14 :         for (int iLine = 0; iLine < nYSize; ++iLine)
     389             :         {
     390          72 :             for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
     391             :             {
     392             :                 // Source raster pixels may be obtained with GetSrcVal macro.
     393          60 :                 const double dfReal = GetSrcVal(pReal, eSrcType, ii);
     394          60 :                 const double dfImag = GetSrcVal(pImag, eSrcType, ii);
     395             : 
     396          60 :                 const double dfPixVal = atan2(dfImag, dfReal);
     397             : 
     398          60 :                 GDALCopyWords(&dfPixVal, GDT_Float64, 0,
     399             :                               static_cast<GByte *>(pData) +
     400          60 :                                   static_cast<GSpacing>(nLineSpace) * iLine +
     401          60 :                                   iCol * nPixelSpace,
     402             :                               eBufType, nPixelSpace, 1);
     403             :             }
     404             :         }
     405             :     }
     406           2 :     else if (GDALDataTypeIsInteger(eSrcType) && !GDALDataTypeIsSigned(eSrcType))
     407             :     {
     408           1 :         constexpr double dfZero = 0;
     409           7 :         for (int iLine = 0; iLine < nYSize; ++iLine)
     410             :         {
     411           6 :             GDALCopyWords(&dfZero, GDT_Float64, 0,
     412             :                           static_cast<GByte *>(pData) +
     413           6 :                               static_cast<GSpacing>(nLineSpace) * iLine,
     414             :                           eBufType, nPixelSpace, nXSize);
     415             :         }
     416             :     }
     417             :     else
     418             :     {
     419             :         /* ---- Set pixels ---- */
     420           1 :         size_t ii = 0;
     421           7 :         for (int iLine = 0; iLine < nYSize; ++iLine)
     422             :         {
     423          36 :             for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
     424             :             {
     425          30 :                 const void *const pReal = papoSources[0];
     426             : 
     427             :                 // Source raster pixels may be obtained with GetSrcVal macro.
     428          30 :                 const double dfReal = GetSrcVal(pReal, eSrcType, ii);
     429          30 :                 const double dfPixVal = (dfReal < 0) ? M_PI : 0.0;
     430             : 
     431          30 :                 GDALCopyWords(&dfPixVal, GDT_Float64, 0,
     432             :                               static_cast<GByte *>(pData) +
     433          30 :                                   static_cast<GSpacing>(nLineSpace) * iLine +
     434          30 :                                   iCol * nPixelSpace,
     435             :                               eBufType, nPixelSpace, 1);
     436             :             }
     437             :         }
     438             :     }
     439             : 
     440             :     /* ---- Return success ---- */
     441           4 :     return CE_None;
     442             : }  // PhasePixelFunc
     443             : 
     444           4 : static CPLErr ConjPixelFunc(void **papoSources, int nSources, void *pData,
     445             :                             int nXSize, int nYSize, GDALDataType eSrcType,
     446             :                             GDALDataType eBufType, int nPixelSpace,
     447             :                             int nLineSpace)
     448             : {
     449             :     /* ---- Init ---- */
     450           4 :     if (nSources != 1)
     451           1 :         return CE_Failure;
     452             : 
     453           3 :     if (GDALDataTypeIsComplex(eSrcType) && GDALDataTypeIsComplex(eBufType))
     454             :     {
     455           2 :         const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
     456           2 :         const void *const pReal = papoSources[0];
     457           2 :         const void *const pImag =
     458           2 :             static_cast<GByte *>(papoSources[0]) + nOffset;
     459             : 
     460             :         /* ---- Set pixels ---- */
     461           2 :         size_t ii = 0;
     462          14 :         for (int iLine = 0; iLine < nYSize; ++iLine)
     463             :         {
     464          72 :             for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
     465             :             {
     466             :                 // Source raster pixels may be obtained with GetSrcVal macro.
     467             :                 const double adfPixVal[2] = {
     468          60 :                     +GetSrcVal(pReal, eSrcType, ii),  // re
     469         120 :                     -GetSrcVal(pImag, eSrcType, ii)   // im
     470          60 :                 };
     471             : 
     472          60 :                 GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
     473             :                               static_cast<GByte *>(pData) +
     474          60 :                                   static_cast<GSpacing>(nLineSpace) * iLine +
     475          60 :                                   iCol * nPixelSpace,
     476             :                               eBufType, nPixelSpace, 1);
     477             :             }
     478             :         }
     479             :     }
     480             :     else
     481             :     {
     482             :         // No complex data type.
     483           1 :         return RealPixelFunc(papoSources, nSources, pData, nXSize, nYSize,
     484           1 :                              eSrcType, eBufType, nPixelSpace, nLineSpace);
     485             :     }
     486             : 
     487             :     /* ---- Return success ---- */
     488           2 :     return CE_None;
     489             : }  // ConjPixelFunc
     490             : 
     491             : static const char pszSumPixelFuncMetadata[] =
     492             :     "<PixelFunctionArgumentsList>"
     493             :     "   <Argument name='k' description='Optional constant term' type='double' "
     494             :     "default='0.0' />"
     495             :     "</PixelFunctionArgumentsList>";
     496             : 
     497           4 : static CPLErr SumPixelFunc(void **papoSources, int nSources, void *pData,
     498             :                            int nXSize, int nYSize, GDALDataType eSrcType,
     499             :                            GDALDataType eBufType, int nPixelSpace,
     500             :                            int nLineSpace, CSLConstList papszArgs)
     501             : {
     502             :     /* ---- Init ---- */
     503           4 :     if (nSources < 2)
     504           1 :         return CE_Failure;
     505             : 
     506           3 :     double dfK = 0.0;
     507           3 :     if (FetchDoubleArg(papszArgs, "k", &dfK, &dfK) != CE_None)
     508           0 :         return CE_Failure;
     509             : 
     510             :     /* ---- Set pixels ---- */
     511           3 :     if (GDALDataTypeIsComplex(eSrcType))
     512             :     {
     513           1 :         const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
     514             : 
     515             :         /* ---- Set pixels ---- */
     516           1 :         size_t ii = 0;
     517           7 :         for (int iLine = 0; iLine < nYSize; ++iLine)
     518             :         {
     519          36 :             for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
     520             :             {
     521          30 :                 double adfSum[2] = {dfK, 0.0};
     522             : 
     523         120 :                 for (int iSrc = 0; iSrc < nSources; ++iSrc)
     524             :                 {
     525          90 :                     const void *const pReal = papoSources[iSrc];
     526          90 :                     const void *const pImag =
     527          90 :                         static_cast<const GByte *>(pReal) + nOffset;
     528             : 
     529             :                     // Source raster pixels may be obtained with GetSrcVal
     530             :                     // macro.
     531          90 :                     adfSum[0] += GetSrcVal(pReal, eSrcType, ii);
     532          90 :                     adfSum[1] += GetSrcVal(pImag, eSrcType, ii);
     533             :                 }
     534             : 
     535          30 :                 GDALCopyWords(adfSum, GDT_CFloat64, 0,
     536             :                               static_cast<GByte *>(pData) +
     537          30 :                                   static_cast<GSpacing>(nLineSpace) * iLine +
     538          30 :                                   iCol * nPixelSpace,
     539             :                               eBufType, nPixelSpace, 1);
     540             :             }
     541             :         }
     542             :     }
     543             :     else
     544             :     {
     545             :         /* ---- Set pixels ---- */
     546           2 :         size_t ii = 0;
     547          42 :         for (int iLine = 0; iLine < nYSize; ++iLine)
     548             :         {
     549         840 :             for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
     550             :             {
     551         800 :                 double dfSum = dfK;  // Not complex.
     552             : 
     553        3200 :                 for (int iSrc = 0; iSrc < nSources; ++iSrc)
     554             :                 {
     555             :                     // Source raster pixels may be obtained with GetSrcVal
     556             :                     // macro.
     557        2400 :                     dfSum += GetSrcVal(papoSources[iSrc], eSrcType, ii);
     558             :                 }
     559             : 
     560         800 :                 GDALCopyWords(&dfSum, GDT_Float64, 0,
     561             :                               static_cast<GByte *>(pData) +
     562         800 :                                   static_cast<GSpacing>(nLineSpace) * iLine +
     563         800 :                                   iCol * nPixelSpace,
     564             :                               eBufType, nPixelSpace, 1);
     565             :             }
     566             :         }
     567             :     }
     568             : 
     569             :     /* ---- Return success ---- */
     570           3 :     return CE_None;
     571             : } /* SumPixelFunc */
     572             : 
     573           3 : static CPLErr DiffPixelFunc(void **papoSources, int nSources, void *pData,
     574             :                             int nXSize, int nYSize, GDALDataType eSrcType,
     575             :                             GDALDataType eBufType, int nPixelSpace,
     576             :                             int nLineSpace)
     577             : {
     578             :     /* ---- Init ---- */
     579           3 :     if (nSources != 2)
     580           1 :         return CE_Failure;
     581             : 
     582           2 :     if (GDALDataTypeIsComplex(eSrcType))
     583             :     {
     584           1 :         const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
     585           1 :         const void *const pReal0 = papoSources[0];
     586           1 :         const void *const pImag0 =
     587           1 :             static_cast<GByte *>(papoSources[0]) + nOffset;
     588           1 :         const void *const pReal1 = papoSources[1];
     589           1 :         const void *const pImag1 =
     590           1 :             static_cast<GByte *>(papoSources[1]) + nOffset;
     591             : 
     592             :         /* ---- Set pixels ---- */
     593           1 :         size_t ii = 0;
     594           7 :         for (int iLine = 0; iLine < nYSize; ++iLine)
     595             :         {
     596          36 :             for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
     597             :             {
     598             :                 // Source raster pixels may be obtained with GetSrcVal macro.
     599          30 :                 double adfPixVal[2] = {GetSrcVal(pReal0, eSrcType, ii) -
     600          30 :                                            GetSrcVal(pReal1, eSrcType, ii),
     601          60 :                                        GetSrcVal(pImag0, eSrcType, ii) -
     602          30 :                                            GetSrcVal(pImag1, eSrcType, ii)};
     603             : 
     604          30 :                 GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
     605             :                               static_cast<GByte *>(pData) +
     606          30 :                                   static_cast<GSpacing>(nLineSpace) * iLine +
     607          30 :                                   iCol * nPixelSpace,
     608             :                               eBufType, nPixelSpace, 1);
     609             :             }
     610             :         }
     611             :     }
     612             :     else
     613             :     {
     614             :         /* ---- Set pixels ---- */
     615           1 :         size_t ii = 0;
     616           7 :         for (int iLine = 0; iLine < nYSize; ++iLine)
     617             :         {
     618          36 :             for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
     619             :             {
     620             :                 // Source raster pixels may be obtained with GetSrcVal macro.
     621             :                 // Not complex.
     622             :                 const double dfPixVal =
     623          30 :                     GetSrcVal(papoSources[0], eSrcType, ii) -
     624          30 :                     GetSrcVal(papoSources[1], eSrcType, ii);
     625             : 
     626          30 :                 GDALCopyWords(&dfPixVal, GDT_Float64, 0,
     627             :                               static_cast<GByte *>(pData) +
     628          30 :                                   static_cast<GSpacing>(nLineSpace) * iLine +
     629          30 :                                   iCol * nPixelSpace,
     630             :                               eBufType, nPixelSpace, 1);
     631             :             }
     632             :         }
     633             :     }
     634             : 
     635             :     /* ---- Return success ---- */
     636           2 :     return CE_None;
     637             : }  // DiffPixelFunc
     638             : 
     639             : static const char pszMulPixelFuncMetadata[] =
     640             :     "<PixelFunctionArgumentsList>"
     641             :     "   <Argument name='k' description='Optional constant factor' "
     642             :     "type='double' default='1.0' />"
     643             :     "</PixelFunctionArgumentsList>";
     644             : 
     645           4 : static CPLErr MulPixelFunc(void **papoSources, int nSources, void *pData,
     646             :                            int nXSize, int nYSize, GDALDataType eSrcType,
     647             :                            GDALDataType eBufType, int nPixelSpace,
     648             :                            int nLineSpace, CSLConstList papszArgs)
     649             : {
     650             :     /* ---- Init ---- */
     651           4 :     if (nSources < 2)
     652           1 :         return CE_Failure;
     653             : 
     654           3 :     double dfK = 1.0;
     655           3 :     if (FetchDoubleArg(papszArgs, "k", &dfK, &dfK) != CE_None)
     656           0 :         return CE_Failure;
     657             : 
     658             :     /* ---- Set pixels ---- */
     659           3 :     if (GDALDataTypeIsComplex(eSrcType))
     660             :     {
     661           1 :         const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
     662             : 
     663             :         /* ---- Set pixels ---- */
     664           1 :         size_t ii = 0;
     665           7 :         for (int iLine = 0; iLine < nYSize; ++iLine)
     666             :         {
     667          36 :             for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
     668             :             {
     669          30 :                 double adfPixVal[2] = {dfK, 0.0};
     670             : 
     671          90 :                 for (int iSrc = 0; iSrc < nSources; ++iSrc)
     672             :                 {
     673          60 :                     const void *const pReal = papoSources[iSrc];
     674          60 :                     const void *const pImag =
     675          60 :                         static_cast<const GByte *>(pReal) + nOffset;
     676             : 
     677          60 :                     const double dfOldR = adfPixVal[0];
     678          60 :                     const double dfOldI = adfPixVal[1];
     679             : 
     680             :                     // Source raster pixels may be obtained with GetSrcVal
     681             :                     // macro.
     682          60 :                     const double dfNewR = GetSrcVal(pReal, eSrcType, ii);
     683          60 :                     const double dfNewI = GetSrcVal(pImag, eSrcType, ii);
     684             : 
     685          60 :                     adfPixVal[0] = dfOldR * dfNewR - dfOldI * dfNewI;
     686          60 :                     adfPixVal[1] = dfOldR * dfNewI + dfOldI * dfNewR;
     687             :                 }
     688             : 
     689          30 :                 GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
     690             :                               static_cast<GByte *>(pData) +
     691          30 :                                   static_cast<GSpacing>(nLineSpace) * iLine +
     692          30 :                                   iCol * nPixelSpace,
     693             :                               eBufType, nPixelSpace, 1);
     694             :             }
     695             :         }
     696             :     }
     697             :     else
     698             :     {
     699             :         /* ---- Set pixels ---- */
     700           2 :         size_t ii = 0;
     701          42 :         for (int iLine = 0; iLine < nYSize; ++iLine)
     702             :         {
     703         840 :             for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
     704             :             {
     705         800 :                 double dfPixVal = dfK;  // Not complex.
     706             : 
     707        3200 :                 for (int iSrc = 0; iSrc < nSources; ++iSrc)
     708             :                 {
     709             :                     // Source raster pixels may be obtained with GetSrcVal
     710             :                     // macro.
     711        2400 :                     dfPixVal *= GetSrcVal(papoSources[iSrc], eSrcType, ii);
     712             :                 }
     713             : 
     714         800 :                 GDALCopyWords(&dfPixVal, GDT_Float64, 0,
     715             :                               static_cast<GByte *>(pData) +
     716         800 :                                   static_cast<GSpacing>(nLineSpace) * iLine +
     717         800 :                                   iCol * nPixelSpace,
     718             :                               eBufType, nPixelSpace, 1);
     719             :             }
     720             :         }
     721             :     }
     722             : 
     723             :     /* ---- Return success ---- */
     724           3 :     return CE_None;
     725             : }  // MulPixelFunc
     726             : 
     727           2 : static CPLErr DivPixelFunc(void **papoSources, int nSources, void *pData,
     728             :                            int nXSize, int nYSize, GDALDataType eSrcType,
     729             :                            GDALDataType eBufType, int nPixelSpace,
     730             :                            int nLineSpace)
     731             : {
     732             :     /* ---- Init ---- */
     733           2 :     if (nSources != 2)
     734           0 :         return CE_Failure;
     735             : 
     736             :     /* ---- Set pixels ---- */
     737           2 :     if (GDALDataTypeIsComplex(eSrcType))
     738             :     {
     739           1 :         const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
     740           1 :         const void *const pReal0 = papoSources[0];
     741           1 :         const void *const pImag0 =
     742           1 :             static_cast<GByte *>(papoSources[0]) + nOffset;
     743           1 :         const void *const pReal1 = papoSources[1];
     744           1 :         const void *const pImag1 =
     745           1 :             static_cast<GByte *>(papoSources[1]) + nOffset;
     746             : 
     747             :         /* ---- Set pixels ---- */
     748           1 :         size_t ii = 0;
     749           7 :         for (int iLine = 0; iLine < nYSize; ++iLine)
     750             :         {
     751          36 :             for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
     752             :             {
     753             :                 // Source raster pixels may be obtained with GetSrcVal macro.
     754          30 :                 const double dfReal0 = GetSrcVal(pReal0, eSrcType, ii);
     755          30 :                 const double dfReal1 = GetSrcVal(pReal1, eSrcType, ii);
     756          30 :                 const double dfImag0 = GetSrcVal(pImag0, eSrcType, ii);
     757          30 :                 const double dfImag1 = GetSrcVal(pImag1, eSrcType, ii);
     758          30 :                 const double dfAux = dfReal1 * dfReal1 + dfImag1 * dfImag1;
     759             : 
     760             :                 const double adfPixVal[2] = {
     761             :                     dfAux == 0
     762          30 :                         ? std::numeric_limits<double>::infinity()
     763          30 :                         : dfReal0 * dfReal1 / dfAux + dfImag0 * dfImag1 / dfAux,
     764           0 :                     dfAux == 0 ? std::numeric_limits<double>::infinity()
     765          30 :                                : dfReal1 / dfAux * dfImag0 -
     766          30 :                                      dfReal0 * dfImag1 / dfAux};
     767             : 
     768          30 :                 GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
     769             :                               static_cast<GByte *>(pData) +
     770          30 :                                   static_cast<GSpacing>(nLineSpace) * iLine +
     771          30 :                                   iCol * nPixelSpace,
     772             :                               eBufType, nPixelSpace, 1);
     773             :             }
     774             :         }
     775             :     }
     776             :     else
     777             :     {
     778             :         /* ---- Set pixels ---- */
     779           1 :         size_t ii = 0;
     780           7 :         for (int iLine = 0; iLine < nYSize; ++iLine)
     781             :         {
     782          36 :             for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
     783             :             {
     784          30 :                 const double dfVal = GetSrcVal(papoSources[1], eSrcType, ii);
     785             :                 double dfPixVal =
     786             :                     dfVal == 0
     787          60 :                         ? std::numeric_limits<double>::infinity()
     788          30 :                         : GetSrcVal(papoSources[0], eSrcType, ii) / dfVal;
     789             : 
     790          30 :                 GDALCopyWords(&dfPixVal, GDT_Float64, 0,
     791             :                               static_cast<GByte *>(pData) +
     792          30 :                                   static_cast<GSpacing>(nLineSpace) * iLine +
     793          30 :                                   iCol * nPixelSpace,
     794             :                               eBufType, nPixelSpace, 1);
     795             :             }
     796             :         }
     797             :     }
     798             : 
     799             :     /* ---- Return success ---- */
     800           2 :     return CE_None;
     801             : }  // DivPixelFunc
     802             : 
     803           3 : static CPLErr CMulPixelFunc(void **papoSources, int nSources, void *pData,
     804             :                             int nXSize, int nYSize, GDALDataType eSrcType,
     805             :                             GDALDataType eBufType, int nPixelSpace,
     806             :                             int nLineSpace)
     807             : {
     808             :     /* ---- Init ---- */
     809           3 :     if (nSources != 2)
     810           1 :         return CE_Failure;
     811             : 
     812             :     /* ---- Set pixels ---- */
     813           2 :     if (GDALDataTypeIsComplex(eSrcType))
     814             :     {
     815           1 :         const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
     816           1 :         const void *const pReal0 = papoSources[0];
     817           1 :         const void *const pImag0 =
     818           1 :             static_cast<GByte *>(papoSources[0]) + nOffset;
     819           1 :         const void *const pReal1 = papoSources[1];
     820           1 :         const void *const pImag1 =
     821           1 :             static_cast<GByte *>(papoSources[1]) + nOffset;
     822             : 
     823           1 :         size_t ii = 0;
     824           7 :         for (int iLine = 0; iLine < nYSize; ++iLine)
     825             :         {
     826          36 :             for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
     827             :             {
     828             :                 // Source raster pixels may be obtained with GetSrcVal macro.
     829          30 :                 const double dfReal0 = GetSrcVal(pReal0, eSrcType, ii);
     830          30 :                 const double dfReal1 = GetSrcVal(pReal1, eSrcType, ii);
     831          30 :                 const double dfImag0 = GetSrcVal(pImag0, eSrcType, ii);
     832          30 :                 const double dfImag1 = GetSrcVal(pImag1, eSrcType, ii);
     833             :                 const double adfPixVal[2] = {
     834          30 :                     dfReal0 * dfReal1 + dfImag0 * dfImag1,
     835          30 :                     dfReal1 * dfImag0 - dfReal0 * dfImag1};
     836             : 
     837          30 :                 GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
     838             :                               static_cast<GByte *>(pData) +
     839          30 :                                   static_cast<GSpacing>(nLineSpace) * iLine +
     840          30 :                                   iCol * nPixelSpace,
     841             :                               eBufType, nPixelSpace, 1);
     842             :             }
     843             :         }
     844             :     }
     845             :     else
     846             :     {
     847           1 :         size_t ii = 0;
     848          21 :         for (int iLine = 0; iLine < nYSize; ++iLine)
     849             :         {
     850         420 :             for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
     851             :             {
     852             :                 // Source raster pixels may be obtained with GetSrcVal macro.
     853             :                 // Not complex.
     854         400 :                 const double adfPixVal[2] = {
     855         400 :                     GetSrcVal(papoSources[0], eSrcType, ii) *
     856         400 :                         GetSrcVal(papoSources[1], eSrcType, ii),
     857         400 :                     0.0};
     858             : 
     859         400 :                 GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
     860             :                               static_cast<GByte *>(pData) +
     861         400 :                                   static_cast<GSpacing>(nLineSpace) * iLine +
     862         400 :                                   iCol * nPixelSpace,
     863             :                               eBufType, nPixelSpace, 1);
     864             :             }
     865             :         }
     866             :     }
     867             : 
     868             :     /* ---- Return success ---- */
     869           2 :     return CE_None;
     870             : }  // CMulPixelFunc
     871             : 
     872             : static const char pszInvPixelFuncMetadata[] =
     873             :     "<PixelFunctionArgumentsList>"
     874             :     "   <Argument name='k' description='Optional constant factor' "
     875             :     "type='double' default='1.0' />"
     876             :     "</PixelFunctionArgumentsList>";
     877             : 
     878           6 : static CPLErr InvPixelFunc(void **papoSources, int nSources, void *pData,
     879             :                            int nXSize, int nYSize, GDALDataType eSrcType,
     880             :                            GDALDataType eBufType, int nPixelSpace,
     881             :                            int nLineSpace, CSLConstList papszArgs)
     882             : {
     883             :     /* ---- Init ---- */
     884           6 :     if (nSources != 1)
     885           1 :         return CE_Failure;
     886             : 
     887           5 :     double dfK = 1.0;
     888           5 :     if (FetchDoubleArg(papszArgs, "k", &dfK, &dfK) != CE_None)
     889           0 :         return CE_Failure;
     890             : 
     891             :     /* ---- Set pixels ---- */
     892           5 :     if (GDALDataTypeIsComplex(eSrcType))
     893             :     {
     894           2 :         const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
     895           2 :         const void *const pReal = papoSources[0];
     896           2 :         const void *const pImag =
     897           2 :             static_cast<GByte *>(papoSources[0]) + nOffset;
     898             : 
     899           2 :         size_t ii = 0;
     900           9 :         for (int iLine = 0; iLine < nYSize; ++iLine)
     901             :         {
     902          38 :             for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
     903             :             {
     904             :                 // Source raster pixels may be obtained with GetSrcVal macro.
     905          31 :                 const double dfReal = GetSrcVal(pReal, eSrcType, ii);
     906          31 :                 const double dfImag = GetSrcVal(pImag, eSrcType, ii);
     907          31 :                 const double dfAux = dfReal * dfReal + dfImag * dfImag;
     908             :                 const double adfPixVal[2] = {
     909          31 :                     dfAux == 0 ? std::numeric_limits<double>::infinity()
     910          30 :                                : dfK * dfReal / dfAux,
     911           1 :                     dfAux == 0 ? std::numeric_limits<double>::infinity()
     912          31 :                                : -dfK * dfImag / dfAux};
     913             : 
     914          31 :                 GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
     915             :                               static_cast<GByte *>(pData) +
     916          31 :                                   static_cast<GSpacing>(nLineSpace) * iLine +
     917          31 :                                   iCol * nPixelSpace,
     918             :                               eBufType, nPixelSpace, 1);
     919             :             }
     920             :         }
     921             :     }
     922             :     else
     923             :     {
     924             :         /* ---- Set pixels ---- */
     925           3 :         size_t ii = 0;
     926          44 :         for (int iLine = 0; iLine < nYSize; ++iLine)
     927             :         {
     928         842 :             for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
     929             :             {
     930             :                 // Source raster pixels may be obtained with GetSrcVal macro.
     931             :                 // Not complex.
     932         801 :                 const double dfVal = GetSrcVal(papoSources[0], eSrcType, ii);
     933             :                 const double dfPixVal =
     934         801 :                     dfVal == 0 ? std::numeric_limits<double>::infinity()
     935         800 :                                : dfK / dfVal;
     936             : 
     937         801 :                 GDALCopyWords(&dfPixVal, GDT_Float64, 0,
     938             :                               static_cast<GByte *>(pData) +
     939         801 :                                   static_cast<GSpacing>(nLineSpace) * iLine +
     940         801 :                                   iCol * nPixelSpace,
     941             :                               eBufType, nPixelSpace, 1);
     942             :             }
     943             :         }
     944             :     }
     945             : 
     946             :     /* ---- Return success ---- */
     947           5 :     return CE_None;
     948             : }  // InvPixelFunc
     949             : 
     950           4 : static CPLErr IntensityPixelFunc(void **papoSources, int nSources, void *pData,
     951             :                                  int nXSize, int nYSize, GDALDataType eSrcType,
     952             :                                  GDALDataType eBufType, int nPixelSpace,
     953             :                                  int nLineSpace)
     954             : {
     955             :     /* ---- Init ---- */
     956           4 :     if (nSources != 1)
     957           1 :         return CE_Failure;
     958             : 
     959           3 :     if (GDALDataTypeIsComplex(eSrcType))
     960             :     {
     961           2 :         const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
     962           2 :         const void *const pReal = papoSources[0];
     963           2 :         const void *const pImag =
     964           2 :             static_cast<GByte *>(papoSources[0]) + nOffset;
     965             : 
     966             :         /* ---- Set pixels ---- */
     967           2 :         size_t ii = 0;
     968          14 :         for (int iLine = 0; iLine < nYSize; ++iLine)
     969             :         {
     970          72 :             for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
     971             :             {
     972             :                 // Source raster pixels may be obtained with GetSrcVal macro.
     973          60 :                 const double dfReal = GetSrcVal(pReal, eSrcType, ii);
     974          60 :                 const double dfImag = GetSrcVal(pImag, eSrcType, ii);
     975             : 
     976          60 :                 const double dfPixVal = dfReal * dfReal + dfImag * dfImag;
     977             : 
     978          60 :                 GDALCopyWords(&dfPixVal, GDT_Float64, 0,
     979             :                               static_cast<GByte *>(pData) +
     980          60 :                                   static_cast<GSpacing>(nLineSpace) * iLine +
     981          60 :                                   iCol * nPixelSpace,
     982             :                               eBufType, nPixelSpace, 1);
     983             :             }
     984             :         }
     985             :     }
     986             :     else
     987             :     {
     988             :         /* ---- Set pixels ---- */
     989           1 :         size_t ii = 0;
     990          21 :         for (int iLine = 0; iLine < nYSize; ++iLine)
     991             :         {
     992         420 :             for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
     993             :             {
     994             :                 // Source raster pixels may be obtained with GetSrcVal macro.
     995         400 :                 double dfPixVal = GetSrcVal(papoSources[0], eSrcType, ii);
     996         400 :                 dfPixVal *= dfPixVal;
     997             : 
     998         400 :                 GDALCopyWords(&dfPixVal, GDT_Float64, 0,
     999             :                               static_cast<GByte *>(pData) +
    1000         400 :                                   static_cast<GSpacing>(nLineSpace) * iLine +
    1001         400 :                                   iCol * nPixelSpace,
    1002             :                               eBufType, nPixelSpace, 1);
    1003             :             }
    1004             :         }
    1005             :     }
    1006             : 
    1007             :     /* ---- Return success ---- */
    1008           3 :     return CE_None;
    1009             : }  // IntensityPixelFunc
    1010             : 
    1011           2 : static CPLErr SqrtPixelFunc(void **papoSources, int nSources, void *pData,
    1012             :                             int nXSize, int nYSize, GDALDataType eSrcType,
    1013             :                             GDALDataType eBufType, int nPixelSpace,
    1014             :                             int nLineSpace)
    1015             : {
    1016             :     /* ---- Init ---- */
    1017           2 :     if (nSources != 1)
    1018           1 :         return CE_Failure;
    1019           1 :     if (GDALDataTypeIsComplex(eSrcType))
    1020           0 :         return CE_Failure;
    1021             : 
    1022             :     /* ---- Set pixels ---- */
    1023           1 :     size_t ii = 0;
    1024          21 :     for (int iLine = 0; iLine < nYSize; ++iLine)
    1025             :     {
    1026         420 :         for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
    1027             :         {
    1028             :             // Source raster pixels may be obtained with GetSrcVal macro.
    1029             :             const double dfPixVal =
    1030         400 :                 sqrt(GetSrcVal(papoSources[0], eSrcType, ii));
    1031             : 
    1032         400 :             GDALCopyWords(&dfPixVal, GDT_Float64, 0,
    1033             :                           static_cast<GByte *>(pData) +
    1034         400 :                               static_cast<GSpacing>(nLineSpace) * iLine +
    1035         400 :                               iCol * nPixelSpace,
    1036             :                           eBufType, nPixelSpace, 1);
    1037             :         }
    1038             :     }
    1039             : 
    1040             :     /* ---- Return success ---- */
    1041           1 :     return CE_None;
    1042             : }  // SqrtPixelFunc
    1043             : 
    1044          11 : static CPLErr Log10PixelFuncHelper(void **papoSources, int nSources,
    1045             :                                    void *pData, int nXSize, int nYSize,
    1046             :                                    GDALDataType eSrcType, GDALDataType eBufType,
    1047             :                                    int nPixelSpace, int nLineSpace, double fact)
    1048             : {
    1049             :     /* ---- Init ---- */
    1050          11 :     if (nSources != 1)
    1051           2 :         return CE_Failure;
    1052             : 
    1053           9 :     if (GDALDataTypeIsComplex(eSrcType))
    1054             :     {
    1055             :         // Complex input datatype.
    1056           5 :         const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
    1057           5 :         const void *const pReal = papoSources[0];
    1058           5 :         const void *const pImag =
    1059           5 :             static_cast<GByte *>(papoSources[0]) + nOffset;
    1060             : 
    1061             :         /* We should compute fact * log10( sqrt( dfReal * dfReal + dfImag *
    1062             :          * dfImag ) ) */
    1063             :         /* Given that log10(sqrt(x)) = 0.5 * log10(x) */
    1064             :         /* we can remove the sqrt() by multiplying fact by 0.5 */
    1065           5 :         fact *= 0.5;
    1066             : 
    1067             :         /* ---- Set pixels ---- */
    1068           5 :         size_t ii = 0;
    1069          35 :         for (int iLine = 0; iLine < nYSize; ++iLine)
    1070             :         {
    1071         180 :             for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
    1072             :             {
    1073             :                 // Source raster pixels may be obtained with GetSrcVal macro.
    1074         150 :                 const double dfReal = GetSrcVal(pReal, eSrcType, ii);
    1075         150 :                 const double dfImag = GetSrcVal(pImag, eSrcType, ii);
    1076             : 
    1077             :                 const double dfPixVal =
    1078         150 :                     fact * log10(dfReal * dfReal + dfImag * dfImag);
    1079             : 
    1080         150 :                 GDALCopyWords(&dfPixVal, GDT_Float64, 0,
    1081             :                               static_cast<GByte *>(pData) +
    1082         150 :                                   static_cast<GSpacing>(nLineSpace) * iLine +
    1083         150 :                                   iCol * nPixelSpace,
    1084             :                               eBufType, nPixelSpace, 1);
    1085             :             }
    1086             :         }
    1087             :     }
    1088             :     else
    1089             :     {
    1090             :         /* ---- Set pixels ---- */
    1091           4 :         size_t ii = 0;
    1092          84 :         for (int iLine = 0; iLine < nYSize; ++iLine)
    1093             :         {
    1094        1680 :             for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
    1095             :             {
    1096             :                 // Source raster pixels may be obtained with GetSrcVal macro.
    1097             :                 const double dfPixVal =
    1098        1600 :                     fact * log10(fabs(GetSrcVal(papoSources[0], eSrcType, ii)));
    1099             : 
    1100        1600 :                 GDALCopyWords(&dfPixVal, GDT_Float64, 0,
    1101             :                               static_cast<GByte *>(pData) +
    1102        1600 :                                   static_cast<GSpacing>(nLineSpace) * iLine +
    1103        1600 :                                   iCol * nPixelSpace,
    1104             :                               eBufType, nPixelSpace, 1);
    1105             :             }
    1106             :         }
    1107             :     }
    1108             : 
    1109             :     /* ---- Return success ---- */
    1110           9 :     return CE_None;
    1111             : }  // Log10PixelFuncHelper
    1112             : 
    1113           4 : static CPLErr Log10PixelFunc(void **papoSources, int nSources, void *pData,
    1114             :                              int nXSize, int nYSize, GDALDataType eSrcType,
    1115             :                              GDALDataType eBufType, int nPixelSpace,
    1116             :                              int nLineSpace)
    1117             : {
    1118           4 :     return Log10PixelFuncHelper(papoSources, nSources, pData, nXSize, nYSize,
    1119             :                                 eSrcType, eBufType, nPixelSpace, nLineSpace,
    1120           4 :                                 1.0);
    1121             : }  // Log10PixelFunc
    1122             : 
    1123             : static const char pszDBPixelFuncMetadata[] =
    1124             :     "<PixelFunctionArgumentsList>"
    1125             :     "   <Argument name='fact' description='Factor' type='double' "
    1126             :     "default='20.0' />"
    1127             :     "</PixelFunctionArgumentsList>";
    1128             : 
    1129           7 : static CPLErr DBPixelFunc(void **papoSources, int nSources, void *pData,
    1130             :                           int nXSize, int nYSize, GDALDataType eSrcType,
    1131             :                           GDALDataType eBufType, int nPixelSpace,
    1132             :                           int nLineSpace, CSLConstList papszArgs)
    1133             : {
    1134           7 :     double dfFact = 20.;
    1135           7 :     if (FetchDoubleArg(papszArgs, "fact", &dfFact, &dfFact) != CE_None)
    1136           0 :         return CE_Failure;
    1137             : 
    1138           7 :     return Log10PixelFuncHelper(papoSources, nSources, pData, nXSize, nYSize,
    1139             :                                 eSrcType, eBufType, nPixelSpace, nLineSpace,
    1140           7 :                                 dfFact);
    1141             : }  // DBPixelFunc
    1142             : 
    1143           7 : static CPLErr ExpPixelFuncHelper(void **papoSources, int nSources, void *pData,
    1144             :                                  int nXSize, int nYSize, GDALDataType eSrcType,
    1145             :                                  GDALDataType eBufType, int nPixelSpace,
    1146             :                                  int nLineSpace, double base, double fact)
    1147             : {
    1148             :     /* ---- Init ---- */
    1149           7 :     if (nSources != 1)
    1150           2 :         return CE_Failure;
    1151           5 :     if (GDALDataTypeIsComplex(eSrcType))
    1152           0 :         return CE_Failure;
    1153             : 
    1154             :     /* ---- Set pixels ---- */
    1155           5 :     size_t ii = 0;
    1156         105 :     for (int iLine = 0; iLine < nYSize; ++iLine)
    1157             :     {
    1158        2100 :         for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
    1159             :         {
    1160             :             // Source raster pixels may be obtained with GetSrcVal macro.
    1161             :             const double dfPixVal =
    1162        2000 :                 pow(base, GetSrcVal(papoSources[0], eSrcType, ii) * fact);
    1163             : 
    1164        2000 :             GDALCopyWords(&dfPixVal, GDT_Float64, 0,
    1165             :                           static_cast<GByte *>(pData) +
    1166        2000 :                               static_cast<GSpacing>(nLineSpace) * iLine +
    1167        2000 :                               iCol * nPixelSpace,
    1168             :                           eBufType, nPixelSpace, 1);
    1169             :         }
    1170             :     }
    1171             : 
    1172             :     /* ---- Return success ---- */
    1173           5 :     return CE_None;
    1174             : }  // ExpPixelFuncHelper
    1175             : 
    1176             : static const char pszExpPixelFuncMetadata[] =
    1177             :     "<PixelFunctionArgumentsList>"
    1178             :     "   <Argument name='base' description='Base' type='double' "
    1179             :     "default='2.7182818284590452353602874713526624' />"
    1180             :     "   <Argument name='fact' description='Factor' type='double' default='1' />"
    1181             :     "</PixelFunctionArgumentsList>";
    1182             : 
    1183           3 : static CPLErr ExpPixelFunc(void **papoSources, int nSources, void *pData,
    1184             :                            int nXSize, int nYSize, GDALDataType eSrcType,
    1185             :                            GDALDataType eBufType, int nPixelSpace,
    1186             :                            int nLineSpace, CSLConstList papszArgs)
    1187             : {
    1188           3 :     double dfBase = 2.7182818284590452353602874713526624;
    1189           3 :     double dfFact = 1.;
    1190             : 
    1191           3 :     if (FetchDoubleArg(papszArgs, "base", &dfBase, &dfBase) != CE_None)
    1192           0 :         return CE_Failure;
    1193             : 
    1194           3 :     if (FetchDoubleArg(papszArgs, "fact", &dfFact, &dfFact) != CE_None)
    1195           0 :         return CE_Failure;
    1196             : 
    1197           3 :     return ExpPixelFuncHelper(papoSources, nSources, pData, nXSize, nYSize,
    1198             :                               eSrcType, eBufType, nPixelSpace, nLineSpace,
    1199           3 :                               dfBase, dfFact);
    1200             : }  // ExpPixelFunc
    1201             : 
    1202           2 : static CPLErr dB2AmpPixelFunc(void **papoSources, int nSources, void *pData,
    1203             :                               int nXSize, int nYSize, GDALDataType eSrcType,
    1204             :                               GDALDataType eBufType, int nPixelSpace,
    1205             :                               int nLineSpace)
    1206             : {
    1207           2 :     return ExpPixelFuncHelper(papoSources, nSources, pData, nXSize, nYSize,
    1208             :                               eSrcType, eBufType, nPixelSpace, nLineSpace, 10.0,
    1209           2 :                               1. / 20);
    1210             : }  // dB2AmpPixelFunc
    1211             : 
    1212           2 : static CPLErr dB2PowPixelFunc(void **papoSources, int nSources, void *pData,
    1213             :                               int nXSize, int nYSize, GDALDataType eSrcType,
    1214             :                               GDALDataType eBufType, int nPixelSpace,
    1215             :                               int nLineSpace)
    1216             : {
    1217           2 :     return ExpPixelFuncHelper(papoSources, nSources, pData, nXSize, nYSize,
    1218             :                               eSrcType, eBufType, nPixelSpace, nLineSpace, 10.0,
    1219           2 :                               1. / 10);
    1220             : }  // dB2PowPixelFunc
    1221             : 
    1222             : static const char pszPowPixelFuncMetadata[] =
    1223             :     "<PixelFunctionArgumentsList>"
    1224             :     "   <Argument name='power' description='Exponent' type='double' "
    1225             :     "mandatory='1' />"
    1226             :     "</PixelFunctionArgumentsList>";
    1227             : 
    1228           1 : static CPLErr PowPixelFunc(void **papoSources, int nSources, void *pData,
    1229             :                            int nXSize, int nYSize, GDALDataType eSrcType,
    1230             :                            GDALDataType eBufType, int nPixelSpace,
    1231             :                            int nLineSpace, CSLConstList papszArgs)
    1232             : {
    1233             :     /* ---- Init ---- */
    1234           1 :     if (nSources != 1)
    1235           0 :         return CE_Failure;
    1236           1 :     if (GDALDataTypeIsComplex(eSrcType))
    1237           0 :         return CE_Failure;
    1238             : 
    1239             :     double power;
    1240           1 :     if (FetchDoubleArg(papszArgs, "power", &power) != CE_None)
    1241           0 :         return CE_Failure;
    1242             : 
    1243             :     /* ---- Set pixels ---- */
    1244           1 :     size_t ii = 0;
    1245          21 :     for (int iLine = 0; iLine < nYSize; ++iLine)
    1246             :     {
    1247         420 :         for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
    1248             :         {
    1249             :             const double dfPixVal =
    1250         400 :                 std::pow(GetSrcVal(papoSources[0], eSrcType, ii), power);
    1251             : 
    1252         400 :             GDALCopyWords(&dfPixVal, GDT_Float64, 0,
    1253             :                           static_cast<GByte *>(pData) +
    1254         400 :                               static_cast<GSpacing>(nLineSpace) * iLine +
    1255         400 :                               iCol * nPixelSpace,
    1256             :                           eBufType, nPixelSpace, 1);
    1257             :         }
    1258             :     }
    1259             : 
    1260             :     /* ---- Return success ---- */
    1261           1 :     return CE_None;
    1262             : }
    1263             : 
    1264             : // Given nt intervals spaced by dt and beginning at t0, return the index of
    1265             : // the lower bound of the interval that should be used to
    1266             : // interpolate/extrapolate a value for t.
    1267           7 : static std::size_t intervalLeft(double t0, double dt, std::size_t nt, double t)
    1268             : {
    1269           7 :     if (t < t0)
    1270             :     {
    1271           2 :         return 0;
    1272             :     }
    1273             : 
    1274           5 :     std::size_t n = static_cast<std::size_t>((t - t0) / dt);
    1275             : 
    1276           5 :     if (n >= nt - 1)
    1277             :     {
    1278           3 :         return nt - 2;
    1279             :     }
    1280             : 
    1281           2 :     return n;
    1282             : }
    1283             : 
    1284          16 : static double InterpolateLinear(double dfX0, double dfX1, double dfY0,
    1285             :                                 double dfY1, double dfX)
    1286             : {
    1287          16 :     return dfY0 + (dfX - dfX0) * (dfY1 - dfY0) / (dfX1 - dfX0);
    1288             : }
    1289             : 
    1290          12 : static double InterpolateExponential(double dfX0, double dfX1, double dfY0,
    1291             :                                      double dfY1, double dfX)
    1292             : {
    1293          12 :     const double r = std::log(dfY1 / dfY0) / (dfX1 - dfX0);
    1294          12 :     return dfY0 * std::exp(r * (dfX - dfX0));
    1295             : }
    1296             : 
    1297             : static const char pszInterpolatePixelFuncMetadata[] =
    1298             :     "<PixelFunctionArgumentsList>"
    1299             :     "   <Argument name='t0' description='t0' type='double' mandatory='1' />"
    1300             :     "   <Argument name='dt' description='dt' type='double' mandatory='1' />"
    1301             :     "   <Argument name='t' description='t' type='double' mandatory='1' />"
    1302             :     "</PixelFunctionArgumentsList>";
    1303             : 
    1304             : template <decltype(InterpolateLinear) InterpolationFunction>
    1305           7 : CPLErr InterpolatePixelFunc(void **papoSources, int nSources, void *pData,
    1306             :                             int nXSize, int nYSize, GDALDataType eSrcType,
    1307             :                             GDALDataType eBufType, int nPixelSpace,
    1308             :                             int nLineSpace, CSLConstList papszArgs)
    1309             : {
    1310             :     /* ---- Init ---- */
    1311           7 :     if (GDALDataTypeIsComplex(eSrcType))
    1312           0 :         return CE_Failure;
    1313             : 
    1314             :     double dfT0;
    1315           7 :     if (FetchDoubleArg(papszArgs, "t0", &dfT0) == CE_Failure)
    1316           0 :         return CE_Failure;
    1317             : 
    1318             :     double dfT;
    1319           7 :     if (FetchDoubleArg(papszArgs, "t", &dfT) == CE_Failure)
    1320           0 :         return CE_Failure;
    1321             : 
    1322             :     double dfDt;
    1323           7 :     if (FetchDoubleArg(papszArgs, "dt", &dfDt) == CE_Failure)
    1324           0 :         return CE_Failure;
    1325             : 
    1326           7 :     if (nSources < 2)
    1327             :     {
    1328           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    1329             :                  "At least two sources required for interpolation.");
    1330           0 :         return CE_Failure;
    1331             :     }
    1332             : 
    1333           7 :     if (dfT == 0 || !std::isfinite(dfT))
    1334             :     {
    1335           0 :         CPLError(CE_Failure, CPLE_AppDefined, "dt must be finite and non-zero");
    1336           0 :         return CE_Failure;
    1337             :     }
    1338             : 
    1339           7 :     const auto i0 = intervalLeft(dfT0, dfDt, nSources, dfT);
    1340           7 :     const auto i1 = i0 + 1;
    1341           7 :     dfT0 = dfT0 + static_cast<double>(i0) * dfDt;
    1342           7 :     double dfX1 = dfT0 + dfDt;
    1343             : 
    1344             :     /* ---- Set pixels ---- */
    1345           7 :     size_t ii = 0;
    1346          21 :     for (int iLine = 0; iLine < nYSize; ++iLine)
    1347             :     {
    1348          42 :         for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
    1349             :         {
    1350          28 :             const double dfY0 = GetSrcVal(papoSources[i0], eSrcType, ii);
    1351          28 :             const double dfY1 = GetSrcVal(papoSources[i1], eSrcType, ii);
    1352             : 
    1353          28 :             const double dfPixVal =
    1354          28 :                 InterpolationFunction(dfT0, dfX1, dfY0, dfY1, dfT);
    1355             : 
    1356          28 :             GDALCopyWords(&dfPixVal, GDT_Float64, 0,
    1357             :                           static_cast<GByte *>(pData) +
    1358          28 :                               static_cast<GSpacing>(nLineSpace) * iLine +
    1359          28 :                               iCol * nPixelSpace,
    1360             :                           eBufType, nPixelSpace, 1);
    1361             :         }
    1362             :     }
    1363             : 
    1364             :     /* ---- Return success ---- */
    1365           7 :     return CE_None;
    1366             : }
    1367             : 
    1368             : static const char pszReplaceNoDataPixelFuncMetadata[] =
    1369             :     "<PixelFunctionArgumentsList>"
    1370             :     "   <Argument type='builtin' value='NoData' />"
    1371             :     "   <Argument name='to' type='double' description='New NoData value to be "
    1372             :     "replaced' default='nan' />"
    1373             :     "</PixelFunctionArgumentsList>";
    1374             : 
    1375           2 : static CPLErr ReplaceNoDataPixelFunc(void **papoSources, int nSources,
    1376             :                                      void *pData, int nXSize, int nYSize,
    1377             :                                      GDALDataType eSrcType,
    1378             :                                      GDALDataType eBufType, int nPixelSpace,
    1379             :                                      int nLineSpace, CSLConstList papszArgs)
    1380             : {
    1381             :     /* ---- Init ---- */
    1382           2 :     if (nSources != 1)
    1383           0 :         return CE_Failure;
    1384           2 :     if (GDALDataTypeIsComplex(eSrcType))
    1385             :     {
    1386           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    1387             :                  "replace_nodata cannot convert complex data types");
    1388           0 :         return CE_Failure;
    1389             :     }
    1390             : 
    1391           2 :     double dfOldNoData, dfNewNoData = NAN;
    1392           2 :     if (FetchDoubleArg(papszArgs, "NoData", &dfOldNoData) != CE_None)
    1393           0 :         return CE_Failure;
    1394           2 :     if (FetchDoubleArg(papszArgs, "to", &dfNewNoData, &dfNewNoData) != CE_None)
    1395           0 :         return CE_Failure;
    1396             : 
    1397           2 :     if (!GDALDataTypeIsFloating(eBufType) && std::isnan(dfNewNoData))
    1398             :     {
    1399           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    1400             :                  "Using nan requires a floating point type output buffer");
    1401           0 :         return CE_Failure;
    1402             :     }
    1403             : 
    1404             :     /* ---- Set pixels ---- */
    1405           2 :     size_t ii = 0;
    1406         102 :     for (int iLine = 0; iLine < nYSize; ++iLine)
    1407             :     {
    1408        5100 :         for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
    1409             :         {
    1410        5000 :             double dfPixVal = GetSrcVal(papoSources[0], eSrcType, ii);
    1411        5000 :             if (dfPixVal == dfOldNoData || std::isnan(dfPixVal))
    1412        3200 :                 dfPixVal = dfNewNoData;
    1413             : 
    1414        5000 :             GDALCopyWords(&dfPixVal, GDT_Float64, 0,
    1415             :                           static_cast<GByte *>(pData) +
    1416        5000 :                               static_cast<GSpacing>(nLineSpace) * iLine +
    1417        5000 :                               iCol * nPixelSpace,
    1418             :                           eBufType, nPixelSpace, 1);
    1419             :         }
    1420             :     }
    1421             : 
    1422             :     /* ---- Return success ---- */
    1423           2 :     return CE_None;
    1424             : }
    1425             : 
    1426             : static const char pszScalePixelFuncMetadata[] =
    1427             :     "<PixelFunctionArgumentsList>"
    1428             :     "   <Argument type='builtin' value='offset' />"
    1429             :     "   <Argument type='builtin' value='scale' />"
    1430             :     "</PixelFunctionArgumentsList>";
    1431             : 
    1432           1 : static CPLErr ScalePixelFunc(void **papoSources, int nSources, void *pData,
    1433             :                              int nXSize, int nYSize, GDALDataType eSrcType,
    1434             :                              GDALDataType eBufType, int nPixelSpace,
    1435             :                              int nLineSpace, CSLConstList papszArgs)
    1436             : {
    1437             :     /* ---- Init ---- */
    1438           1 :     if (nSources != 1)
    1439           0 :         return CE_Failure;
    1440           1 :     if (GDALDataTypeIsComplex(eSrcType))
    1441             :     {
    1442           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    1443             :                  "scale cannot by applied to complex data types");
    1444           0 :         return CE_Failure;
    1445             :     }
    1446             : 
    1447             :     double dfScale, dfOffset;
    1448           1 :     if (FetchDoubleArg(papszArgs, "scale", &dfScale) != CE_None)
    1449           0 :         return CE_Failure;
    1450           1 :     if (FetchDoubleArg(papszArgs, "offset", &dfOffset) != CE_None)
    1451           0 :         return CE_Failure;
    1452             : 
    1453             :     /* ---- Set pixels ---- */
    1454           1 :     size_t ii = 0;
    1455          21 :     for (int iLine = 0; iLine < nYSize; ++iLine)
    1456             :     {
    1457         420 :         for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
    1458             :         {
    1459             :             const double dfPixVal =
    1460         400 :                 GetSrcVal(papoSources[0], eSrcType, ii) * dfScale + dfOffset;
    1461             : 
    1462         400 :             GDALCopyWords(&dfPixVal, GDT_Float64, 0,
    1463             :                           static_cast<GByte *>(pData) +
    1464         400 :                               static_cast<GSpacing>(nLineSpace) * iLine +
    1465         400 :                               iCol * nPixelSpace,
    1466             :                           eBufType, nPixelSpace, 1);
    1467             :         }
    1468             :     }
    1469             : 
    1470             :     /* ---- Return success ---- */
    1471           1 :     return CE_None;
    1472             : }
    1473             : 
    1474           1 : static CPLErr NormDiffPixelFunc(void **papoSources, int nSources, void *pData,
    1475             :                                 int nXSize, int nYSize, GDALDataType eSrcType,
    1476             :                                 GDALDataType eBufType, int nPixelSpace,
    1477             :                                 int nLineSpace)
    1478             : {
    1479             :     /* ---- Init ---- */
    1480           1 :     if (nSources != 2)
    1481           0 :         return CE_Failure;
    1482             : 
    1483           1 :     if (GDALDataTypeIsComplex(eSrcType))
    1484             :     {
    1485           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    1486             :                  "norm_diff cannot by applied to complex data types");
    1487           0 :         return CE_Failure;
    1488             :     }
    1489             : 
    1490             :     /* ---- Set pixels ---- */
    1491           1 :     size_t ii = 0;
    1492           7 :     for (int iLine = 0; iLine < nYSize; ++iLine)
    1493             :     {
    1494          36 :         for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
    1495             :         {
    1496          30 :             const double dfLeftVal = GetSrcVal(papoSources[0], eSrcType, ii);
    1497          30 :             const double dfRightVal = GetSrcVal(papoSources[1], eSrcType, ii);
    1498             : 
    1499          30 :             const double dfDenom = (dfLeftVal + dfRightVal);
    1500             : 
    1501             :             const double dfPixVal =
    1502          30 :                 dfDenom == 0 ? std::numeric_limits<double>::infinity()
    1503          30 :                              : (dfLeftVal - dfRightVal) / dfDenom;
    1504             : 
    1505          30 :             GDALCopyWords(&dfPixVal, GDT_Float64, 0,
    1506             :                           static_cast<GByte *>(pData) +
    1507          30 :                               static_cast<GSpacing>(nLineSpace) * iLine +
    1508          30 :                               iCol * nPixelSpace,
    1509             :                           eBufType, nPixelSpace, 1);
    1510             :         }
    1511             :     }
    1512             : 
    1513             :     /* ---- Return success ---- */
    1514           1 :     return CE_None;
    1515             : }  // NormDiffPixelFunc
    1516             : 
    1517             : /************************************************************************/
    1518             : /*                   pszMinMaxFuncMetadataNodata                        */
    1519             : /************************************************************************/
    1520             : 
    1521             : static const char pszMinMaxFuncMetadataNodata[] =
    1522             :     "<PixelFunctionArgumentsList>"
    1523             :     "   <Argument type='builtin' value='NoData' optional='true' />"
    1524             :     "   <Argument name='propagateNoData' description='Whether the output value "
    1525             :     "should be NoData as as soon as one source is NoData' type='boolean' "
    1526             :     "default='false' />"
    1527             :     "</PixelFunctionArgumentsList>";
    1528             : 
    1529             : template <class Comparator>
    1530           5 : static CPLErr MinOrMaxPixelFunc(void **papoSources, int nSources, void *pData,
    1531             :                                 int nXSize, int nYSize, GDALDataType eSrcType,
    1532             :                                 GDALDataType eBufType, int nPixelSpace,
    1533             :                                 int nLineSpace, CSLConstList papszArgs)
    1534             : {
    1535             :     /* ---- Init ---- */
    1536           5 :     if (nSources < 2)
    1537           0 :         return CE_Failure;
    1538             : 
    1539           5 :     if (GDALDataTypeIsComplex(eSrcType))
    1540             :     {
    1541           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    1542             :                  "Complex data type not supported for min/max().");
    1543           0 :         return CE_Failure;
    1544             :     }
    1545             : 
    1546           5 :     double dfNoData = std::numeric_limits<double>::quiet_NaN();
    1547           5 :     if (FetchDoubleArg(papszArgs, "NoData", &dfNoData, &dfNoData) != CE_None)
    1548           0 :         return CE_Failure;
    1549           5 :     const bool bPropagateNoData = CPLTestBool(
    1550             :         CSLFetchNameValueDef(papszArgs, "propagateNoData", "false"));
    1551             : 
    1552             :     /* ---- Set pixels ---- */
    1553           5 :     size_t ii = 0;
    1554         255 :     for (int iLine = 0; iLine < nYSize; ++iLine)
    1555             :     {
    1556       12750 :         for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
    1557             :         {
    1558       12500 :             double dfRes = std::numeric_limits<double>::quiet_NaN();
    1559             : 
    1560       43400 :             for (int iSrc = 0; iSrc < nSources; ++iSrc)
    1561             :             {
    1562       32950 :                 const double dfVal = GetSrcVal(papoSources[iSrc], eSrcType, ii);
    1563             : 
    1564       32950 :                 if (std::isnan(dfVal) || dfVal == dfNoData)
    1565             :                 {
    1566       14350 :                     if (bPropagateNoData)
    1567             :                     {
    1568        2050 :                         dfRes = dfNoData;
    1569        2050 :                         break;
    1570             :                     }
    1571             :                 }
    1572       18600 :                 else if (Comparator::compare(dfVal, dfRes))
    1573             :                 {
    1574        9950 :                     dfRes = dfVal;
    1575             :                 }
    1576             :             }
    1577             : 
    1578       12500 :             if (!bPropagateNoData && std::isnan(dfRes))
    1579             :             {
    1580        3200 :                 dfRes = dfNoData;
    1581             :             }
    1582             : 
    1583       12500 :             GDALCopyWords(&dfRes, GDT_Float64, 0,
    1584             :                           static_cast<GByte *>(pData) +
    1585       12500 :                               static_cast<GSpacing>(nLineSpace) * iLine +
    1586       12500 :                               iCol * nPixelSpace,
    1587             :                           eBufType, nPixelSpace, 1);
    1588             :         }
    1589             :     }
    1590             : 
    1591             :     /* ---- Return success ---- */
    1592           5 :     return CE_None;
    1593             : } /* MinOrMaxPixelFunc */
    1594             : 
    1595           2 : static CPLErr MinPixelFunc(void **papoSources, int nSources, void *pData,
    1596             :                            int nXSize, int nYSize, GDALDataType eSrcType,
    1597             :                            GDALDataType eBufType, int nPixelSpace,
    1598             :                            int nLineSpace, CSLConstList papszArgs)
    1599             : {
    1600             :     struct Comparator
    1601             :     {
    1602        2700 :         static bool compare(double x, double resVal)
    1603             :         {
    1604             :             // Written this way to deal with resVal being NaN
    1605        2700 :             return !(x >= resVal);
    1606             :         }
    1607             :     };
    1608             : 
    1609           2 :     return MinOrMaxPixelFunc<Comparator>(papoSources, nSources, pData, nXSize,
    1610             :                                          nYSize, eSrcType, eBufType,
    1611           2 :                                          nPixelSpace, nLineSpace, papszArgs);
    1612             : }
    1613             : 
    1614           3 : static CPLErr MaxPixelFunc(void **papoSources, int nSources, void *pData,
    1615             :                            int nXSize, int nYSize, GDALDataType eSrcType,
    1616             :                            GDALDataType eBufType, int nPixelSpace,
    1617             :                            int nLineSpace, CSLConstList papszArgs)
    1618             : {
    1619             :     struct Comparator
    1620             :     {
    1621       15900 :         static bool compare(double x, double resVal)
    1622             :         {
    1623             :             // Written this way to deal with resVal being NaN
    1624       15900 :             return !(x <= resVal);
    1625             :         }
    1626             :     };
    1627             : 
    1628           3 :     return MinOrMaxPixelFunc<Comparator>(papoSources, nSources, pData, nXSize,
    1629             :                                          nYSize, eSrcType, eBufType,
    1630           3 :                                          nPixelSpace, nLineSpace, papszArgs);
    1631             : }
    1632             : 
    1633             : /************************************************************************/
    1634             : /*                     GDALRegisterDefaultPixelFunc()                   */
    1635             : /************************************************************************/
    1636             : 
    1637             : /**
    1638             :  * This adds a default set of pixel functions to the global list of
    1639             :  * available pixel functions for derived bands:
    1640             :  *
    1641             :  * - "real": extract real part from a single raster band (just a copy if the
    1642             :  *           input is non-complex)
    1643             :  * - "imag": extract imaginary part from a single raster band (0 for
    1644             :  *           non-complex)
    1645             :  * - "complex": make a complex band merging two bands used as real and
    1646             :  *              imag values
    1647             :  * - "polar": make a complex band using input bands for amplitude and
    1648             :  *            phase values (b1 * exp( j * b2 ))
    1649             :  * - "mod": extract module from a single raster band (real or complex)
    1650             :  * - "phase": extract phase from a single raster band [-PI,PI] (0 or PI for
    1651             :               non-complex)
    1652             :  * - "conj": computes the complex conjugate of a single raster band (just a
    1653             :  *           copy if the input is non-complex)
    1654             :  * - "sum": sum 2 or more raster bands
    1655             :  * - "diff": computes the difference between 2 raster bands (b1 - b2)
    1656             :  * - "mul": multiply 2 or more raster bands
    1657             :  * - "div": divide one raster band by another (b1 / b2).
    1658             :  * - "min": minimum value of 2 or more raster bands
    1659             :  * - "max": maximum value of 2 or more raster bands
    1660             :  * - "norm_diff": computes the normalized difference between two raster bands:
    1661             :  *                ``(b1 - b2)/(b1 + b2)``.
    1662             :  * - "cmul": multiply the first band for the complex conjugate of the second
    1663             :  * - "inv": inverse (1./x).
    1664             :  * - "intensity": computes the intensity Re(x*conj(x)) of a single raster band
    1665             :  *                (real or complex)
    1666             :  * - "sqrt": perform the square root of a single raster band (real only)
    1667             :  * - "log10": compute the logarithm (base 10) of the abs of a single raster
    1668             :  *            band (real or complex): log10( abs( x ) )
    1669             :  * - "dB": perform conversion to dB of the abs of a single raster
    1670             :  *         band (real or complex): 20. * log10( abs( x ) ).
    1671             :  *         Note: the optional fact parameter can be set to 10. to get the
    1672             :  *         alternative formula: 10. * log10( abs( x ) )
    1673             :  * - "exp": computes the exponential of each element in the input band ``x``
    1674             :  *          (of real values): ``e ^ x``.
    1675             :  *          The function also accepts two optional parameters: ``base`` and
    1676             :  ``fact``
    1677             :  *          that allow to compute the generalized formula: ``base ^ ( fact *
    1678             :  x)``.
    1679             :  *          Note: this function is the recommended one to perform conversion
    1680             :  *          form logarithmic scale (dB): `` 10. ^ (x / 20.)``, in this case
    1681             :  *          ``base = 10.`` and ``fact = 1./20``
    1682             :  * - "dB2amp": perform scale conversion from logarithmic to linear
    1683             :  *             (amplitude) (i.e. 10 ^ ( x / 20 ) ) of a single raster
    1684             :  *             band (real only).
    1685             :  *             Deprecated in GDAL v3.5. Please use the ``exp`` pixel function
    1686             :  with
    1687             :  *             ``base = 10.`` and ``fact = 0.05`` i.e. ``1./20``
    1688             :  * - "dB2pow": perform scale conversion from logarithmic to linear
    1689             :  *             (power) (i.e. 10 ^ ( x / 10 ) ) of a single raster
    1690             :  *             band (real only)
    1691             :  *             Deprecated in GDAL v3.5. Please use the ``exp`` pixel function
    1692             :  with
    1693             :  *             ``base = 10.`` and ``fact = 0.1`` i.e. ``1./10``
    1694             :  * - "pow": raise a single raster band to a constant power
    1695             :  * - "interpolate_linear": interpolate values between two raster bands
    1696             :  *                         using linear interpolation
    1697             :  * - "interpolate_exp": interpolate values between two raster bands using
    1698             :  *                      exponential interpolation
    1699             :  * - "scale": Apply the RasterBand metadata values of "offset" and "scale"
    1700             :  * - "nan": Convert incoming NoData values to IEEE 754 nan
    1701             :  *
    1702             :  * @see GDALAddDerivedBandPixelFunc
    1703             :  *
    1704             :  * @return CE_None
    1705             :  */
    1706        1218 : CPLErr GDALRegisterDefaultPixelFunc()
    1707             : {
    1708        1218 :     GDALAddDerivedBandPixelFunc("real", RealPixelFunc);
    1709        1218 :     GDALAddDerivedBandPixelFunc("imag", ImagPixelFunc);
    1710        1218 :     GDALAddDerivedBandPixelFunc("complex", ComplexPixelFunc);
    1711        1218 :     GDALAddDerivedBandPixelFuncWithArgs("polar", PolarPixelFunc,
    1712             :                                         pszPolarPixelFuncMetadata);
    1713        1218 :     GDALAddDerivedBandPixelFunc("mod", ModulePixelFunc);
    1714        1218 :     GDALAddDerivedBandPixelFunc("phase", PhasePixelFunc);
    1715        1218 :     GDALAddDerivedBandPixelFunc("conj", ConjPixelFunc);
    1716        1218 :     GDALAddDerivedBandPixelFuncWithArgs("sum", SumPixelFunc,
    1717             :                                         pszSumPixelFuncMetadata);
    1718        1218 :     GDALAddDerivedBandPixelFunc("diff", DiffPixelFunc);
    1719        1218 :     GDALAddDerivedBandPixelFuncWithArgs("mul", MulPixelFunc,
    1720             :                                         pszMulPixelFuncMetadata);
    1721        1218 :     GDALAddDerivedBandPixelFunc("div", DivPixelFunc);
    1722        1218 :     GDALAddDerivedBandPixelFunc("cmul", CMulPixelFunc);
    1723        1218 :     GDALAddDerivedBandPixelFuncWithArgs("inv", InvPixelFunc,
    1724             :                                         pszInvPixelFuncMetadata);
    1725        1218 :     GDALAddDerivedBandPixelFunc("intensity", IntensityPixelFunc);
    1726        1218 :     GDALAddDerivedBandPixelFunc("sqrt", SqrtPixelFunc);
    1727        1218 :     GDALAddDerivedBandPixelFunc("log10", Log10PixelFunc);
    1728        1218 :     GDALAddDerivedBandPixelFuncWithArgs("dB", DBPixelFunc,
    1729             :                                         pszDBPixelFuncMetadata);
    1730        1218 :     GDALAddDerivedBandPixelFuncWithArgs("exp", ExpPixelFunc,
    1731             :                                         pszExpPixelFuncMetadata);
    1732        1218 :     GDALAddDerivedBandPixelFunc("dB2amp",
    1733             :                                 dB2AmpPixelFunc);  // deprecated in v3.5
    1734        1218 :     GDALAddDerivedBandPixelFunc("dB2pow",
    1735             :                                 dB2PowPixelFunc);  // deprecated in v3.5
    1736        1218 :     GDALAddDerivedBandPixelFuncWithArgs("pow", PowPixelFunc,
    1737             :                                         pszPowPixelFuncMetadata);
    1738        1218 :     GDALAddDerivedBandPixelFuncWithArgs("interpolate_linear",
    1739             :                                         InterpolatePixelFunc<InterpolateLinear>,
    1740             :                                         pszInterpolatePixelFuncMetadata);
    1741        1218 :     GDALAddDerivedBandPixelFuncWithArgs(
    1742             :         "interpolate_exp", InterpolatePixelFunc<InterpolateExponential>,
    1743             :         pszInterpolatePixelFuncMetadata);
    1744        1218 :     GDALAddDerivedBandPixelFuncWithArgs("replace_nodata",
    1745             :                                         ReplaceNoDataPixelFunc,
    1746             :                                         pszReplaceNoDataPixelFuncMetadata);
    1747        1218 :     GDALAddDerivedBandPixelFuncWithArgs("scale", ScalePixelFunc,
    1748             :                                         pszScalePixelFuncMetadata);
    1749        1218 :     GDALAddDerivedBandPixelFunc("norm_diff", NormDiffPixelFunc);
    1750        1218 :     GDALAddDerivedBandPixelFuncWithArgs("min", MinPixelFunc,
    1751             :                                         pszMinMaxFuncMetadataNodata);
    1752        1218 :     GDALAddDerivedBandPixelFuncWithArgs("max", MaxPixelFunc,
    1753             :                                         pszMinMaxFuncMetadataNodata);
    1754        1218 :     return CE_None;
    1755             : }

Generated by: LCOV version 1.14