Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  Virtual GDAL Datasets
       4             :  * Purpose:  Implementation of some filter types.
       5             :  * Author:   Frank Warmerdam <warmerdam@pobox.com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2003, Frank Warmerdam <warmerdam@pobox.com>
       9             :  * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include "cpl_port.h"
      15             : #include "vrtdataset.h"
      16             : 
      17             : #include <algorithm>
      18             : #include <cmath>
      19             : #include <cstddef>
      20             : #include <cstdlib>
      21             : #include <cstring>
      22             : #include <limits>
      23             : #include <map>
      24             : 
      25             : #include "cpl_conv.h"
      26             : #include "cpl_error.h"
      27             : #include "cpl_minixml.h"
      28             : #include "cpl_string.h"
      29             : #include "cpl_vsi.h"
      30             : #include "gdal.h"
      31             : #include "gdal_priv.h"
      32             : 
      33             : /*! @cond Doxygen_Suppress */
      34             : 
      35             : /************************************************************************/
      36             : /* ==================================================================== */
      37             : /*                          VRTFilteredSource                           */
      38             : /* ==================================================================== */
      39             : /************************************************************************/
      40             : 
      41             : /************************************************************************/
      42             : /*                         VRTFilteredSource()                          */
      43             : /************************************************************************/
      44             : 
      45          43 : VRTFilteredSource::VRTFilteredSource()
      46          43 :     : m_nSupportedTypesCount(1), m_nExtraEdgePixels(0)
      47             : {
      48         903 :     for (size_t i = 0; i < CPL_ARRAYSIZE(m_aeSupportedTypes); ++i)
      49         860 :         m_aeSupportedTypes[i] = GDT_Unknown;
      50             : 
      51          43 :     m_aeSupportedTypes[0] = GDT_Float32;
      52          43 : }
      53             : 
      54             : /************************************************************************/
      55             : /*                         ~VRTFilteredSource()                         */
      56             : /************************************************************************/
      57             : 
      58          43 : VRTFilteredSource::~VRTFilteredSource()
      59             : {
      60          43 : }
      61             : 
      62             : /************************************************************************/
      63             : /*                         SetExtraEdgePixels()                         */
      64             : /************************************************************************/
      65             : 
      66          41 : void VRTFilteredSource::SetExtraEdgePixels(int nEdgePixels)
      67             : 
      68             : {
      69          41 :     m_nExtraEdgePixels = nEdgePixels;
      70          41 : }
      71             : 
      72             : /************************************************************************/
      73             : /*                   SetFilteringDataTypesSupported()                   */
      74             : /************************************************************************/
      75             : 
      76          43 : void VRTFilteredSource::SetFilteringDataTypesSupported(int nTypeCount,
      77             :                                                        GDALDataType *paeTypes)
      78             : 
      79             : {
      80          43 :     if (nTypeCount >
      81             :         static_cast<int>(sizeof(m_aeSupportedTypes) / sizeof(GDALDataType)))
      82             :     {
      83           0 :         CPLAssert(false);
      84             :         nTypeCount =
      85             :             static_cast<int>(sizeof(m_aeSupportedTypes) / sizeof(GDALDataType));
      86             :     }
      87             : 
      88          43 :     m_nSupportedTypesCount = nTypeCount;
      89          43 :     memcpy(m_aeSupportedTypes, paeTypes, sizeof(GDALDataType) * nTypeCount);
      90          43 : }
      91             : 
      92             : /************************************************************************/
      93             : /*                          IsTypeSupported()                           */
      94             : /************************************************************************/
      95             : 
      96          90 : int VRTFilteredSource::IsTypeSupported(GDALDataType eTestType) const
      97             : 
      98             : {
      99         180 :     for (int i = 0; i < m_nSupportedTypesCount; i++)
     100             :     {
     101          90 :         if (eTestType == m_aeSupportedTypes[i])
     102           0 :             return TRUE;
     103             :     }
     104             : 
     105          90 :     return FALSE;
     106             : }
     107             : 
     108             : /************************************************************************/
     109             : /*                              RasterIO()                              */
     110             : /************************************************************************/
     111             : 
     112          45 : CPLErr VRTFilteredSource::RasterIO(GDALDataType eVRTBandDataType, int nXOff,
     113             :                                    int nYOff, int nXSize, int nYSize,
     114             :                                    void *pData, int nBufXSize, int nBufYSize,
     115             :                                    GDALDataType eBufType, GSpacing nPixelSpace,
     116             :                                    GSpacing nLineSpace,
     117             :                                    GDALRasterIOExtraArg *psExtraArg,
     118             :                                    WorkingState &oWorkingState)
     119             : 
     120             : {
     121             :     /* -------------------------------------------------------------------- */
     122             :     /*      For now we don't support filtered access to non-full            */
     123             :     /*      resolution requests. Just collect the data directly without     */
     124             :     /*      any operator.                                                   */
     125             :     /* -------------------------------------------------------------------- */
     126          45 :     if (nBufXSize != nXSize || nBufYSize != nYSize)
     127             :     {
     128           0 :         return VRTComplexSource::RasterIO(
     129             :             eVRTBandDataType, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize,
     130             :             nBufYSize, eBufType, nPixelSpace, nLineSpace, psExtraArg,
     131           0 :             oWorkingState);
     132             :     }
     133             : 
     134          45 :     double dfXOff = nXOff;
     135          45 :     double dfYOff = nYOff;
     136          45 :     double dfXSize = nXSize;
     137          45 :     double dfYSize = nYSize;
     138          45 :     if (psExtraArg != nullptr && psExtraArg->bFloatingPointWindowValidity)
     139             :     {
     140           4 :         dfXOff = psExtraArg->dfXOff;
     141           4 :         dfYOff = psExtraArg->dfYOff;
     142           4 :         dfXSize = psExtraArg->dfXSize;
     143           4 :         dfYSize = psExtraArg->dfYSize;
     144             :     }
     145             : 
     146             :     // The window we will actually request from the source raster band.
     147          45 :     double dfReqXOff = 0.0;
     148          45 :     double dfReqYOff = 0.0;
     149          45 :     double dfReqXSize = 0.0;
     150          45 :     double dfReqYSize = 0.0;
     151          45 :     int nReqXOff = 0;
     152          45 :     int nReqYOff = 0;
     153          45 :     int nReqXSize = 0;
     154          45 :     int nReqYSize = 0;
     155             : 
     156             :     // The window we will actual set _within_ the pData buffer.
     157          45 :     int nOutXOff = 0;
     158          45 :     int nOutYOff = 0;
     159          45 :     int nOutXSize = 0;
     160          45 :     int nOutYSize = 0;
     161             : 
     162          45 :     bool bError = false;
     163          45 :     if (!GetSrcDstWindow(dfXOff, dfYOff, dfXSize, dfYSize, nBufXSize, nBufYSize,
     164             :                          psExtraArg ? psExtraArg->eResampleAlg
     165             :                                     : GRIORA_NearestNeighbour,
     166             :                          &dfReqXOff, &dfReqYOff, &dfReqXSize, &dfReqYSize,
     167             :                          &nReqXOff, &nReqYOff, &nReqXSize, &nReqYSize,
     168             :                          &nOutXOff, &nOutYOff, &nOutXSize, &nOutYSize, bError))
     169             :     {
     170           0 :         return bError ? CE_Failure : CE_None;
     171             :     }
     172             : 
     173          45 :     pData = reinterpret_cast<GByte *>(pData) + nPixelSpace * nOutXOff +
     174          45 :             nLineSpace * nOutYOff;
     175             : 
     176             :     /* -------------------------------------------------------------------- */
     177             :     /*      Determine the data type we want to request.  We try to match    */
     178             :     /*      the source or destination request, and if both those fail we    */
     179             :     /*      fallback to the first supported type at least as expressive     */
     180             :     /*      as the request.                                                 */
     181             :     /* -------------------------------------------------------------------- */
     182          45 :     GDALDataType eOperDataType = GDT_Unknown;
     183             : 
     184          45 :     if (IsTypeSupported(eBufType))
     185           0 :         eOperDataType = eBufType;
     186             : 
     187          45 :     auto l_band = GetRasterBand();
     188          45 :     if (!l_band)
     189           0 :         return CE_Failure;
     190             : 
     191          90 :     if (eOperDataType == GDT_Unknown &&
     192          45 :         IsTypeSupported(l_band->GetRasterDataType()))
     193           0 :         eOperDataType = l_band->GetRasterDataType();
     194             : 
     195          45 :     if (eOperDataType == GDT_Unknown)
     196             :     {
     197          90 :         for (int i = 0; i < m_nSupportedTypesCount; i++)
     198             :         {
     199          45 :             if (GDALDataTypeUnion(m_aeSupportedTypes[i], eBufType) ==
     200          45 :                 m_aeSupportedTypes[i])
     201             :             {
     202           7 :                 eOperDataType = m_aeSupportedTypes[i];
     203             :             }
     204             :         }
     205             :     }
     206             : 
     207          45 :     if (eOperDataType == GDT_Unknown)
     208             :     {
     209          38 :         eOperDataType = m_aeSupportedTypes[0];
     210             : 
     211          38 :         for (int i = 1; i < m_nSupportedTypesCount; i++)
     212             :         {
     213           0 :             if (GDALGetDataTypeSizeBytes(m_aeSupportedTypes[i]) >
     214           0 :                 GDALGetDataTypeSizeBytes(eOperDataType))
     215             :             {
     216           0 :                 eOperDataType = m_aeSupportedTypes[i];
     217             :             }
     218             :         }
     219             :     }
     220             : 
     221             :     /* -------------------------------------------------------------------- */
     222             :     /*      Allocate the buffer of data into which our imagery will be      */
     223             :     /*      read, with the extra edge pixels as well. This will be the      */
     224             :     /*      source data fed into the filter.                                */
     225             :     /* -------------------------------------------------------------------- */
     226          45 :     if (nOutXSize > INT_MAX - 2 * m_nExtraEdgePixels ||
     227          45 :         nOutYSize > INT_MAX - 2 * m_nExtraEdgePixels)
     228             :     {
     229           0 :         return CE_Failure;
     230             :     }
     231          45 :     const int nExtraXSize = nOutXSize + 2 * m_nExtraEdgePixels;
     232          45 :     const int nExtraYSize = nOutYSize + 2 * m_nExtraEdgePixels;
     233             : 
     234          45 :     GByte *pabyWorkData = static_cast<GByte *>(VSI_MALLOC3_VERBOSE(
     235             :         nExtraXSize, nExtraYSize, GDALGetDataTypeSizeBytes(eOperDataType)));
     236             : 
     237          45 :     if (pabyWorkData == nullptr)
     238             :     {
     239           0 :         return CE_Failure;
     240             :     }
     241             : 
     242          45 :     const GPtrDiff_t nPixelOffset = GDALGetDataTypeSizeBytes(eOperDataType);
     243          45 :     const GPtrDiff_t nLineOffset = nPixelOffset * nExtraXSize;
     244             : 
     245          45 :     int bHasNoData = false;
     246          45 :     const double dfSrcNoDataValue = l_band->GetNoDataValue(&bHasNoData);
     247          45 :     if (bHasNoData)
     248             :     {
     249           5 :         GDALCopyWords64(&dfSrcNoDataValue, GDT_Float64, 0, pabyWorkData,
     250             :                         eOperDataType, static_cast<int>(nPixelOffset),
     251           5 :                         static_cast<size_t>(nExtraXSize) * nExtraYSize);
     252             :     }
     253             :     else
     254             :     {
     255          40 :         memset(pabyWorkData, 0, nLineOffset * nExtraYSize);
     256             :     }
     257             : 
     258             :     /* -------------------------------------------------------------------- */
     259             :     /*      Allocate the output buffer in the same dimensions as the work   */
     260             :     /*      buffer. This allows the filter process to write edge pixels     */
     261             :     /*      if needed for two-pass (separable) filtering.                   */
     262             :     /* -------------------------------------------------------------------- */
     263             :     GByte *pabyOutData = static_cast<GByte *>(
     264          45 :         VSI_MALLOC3_VERBOSE(nExtraXSize, nExtraYSize, nPixelOffset));
     265          45 :     if (pabyOutData == nullptr)
     266             :     {
     267           0 :         CPLFree(pabyWorkData);
     268           0 :         return CE_Failure;
     269             :     }
     270             : 
     271             :     /* -------------------------------------------------------------------- */
     272             :     /*      Figure out the extended window that we want to load.  Note      */
     273             :     /*      that we keep track of the file window as well as the amount     */
     274             :     /*      we will need to edge fill past the edge of the source dataset.  */
     275             :     /* -------------------------------------------------------------------- */
     276          45 :     int nFileXOff = nReqXOff - m_nExtraEdgePixels;
     277          45 :     int nFileYOff = nReqYOff - m_nExtraEdgePixels;
     278          45 :     int nFileXSize = nExtraXSize;
     279          45 :     int nFileYSize = nExtraYSize;
     280             : 
     281          45 :     int nTopFill = 0;
     282          45 :     int nLeftFill = 0;
     283          45 :     int nRightFill = 0;
     284          45 :     int nBottomFill = 0;
     285             : 
     286          45 :     if (nFileXOff < 0)
     287             :     {
     288          45 :         nLeftFill = -nFileXOff;
     289          45 :         nFileXOff = 0;
     290          45 :         nFileXSize -= nLeftFill;
     291             :     }
     292             : 
     293          45 :     if (nFileYOff < 0)
     294             :     {
     295          36 :         nTopFill = -nFileYOff;
     296          36 :         nFileYOff = 0;
     297          36 :         nFileYSize -= nTopFill;
     298             :     }
     299             : 
     300          45 :     if (nFileXOff + nFileXSize > l_band->GetXSize())
     301             :     {
     302          45 :         nRightFill = nFileXOff + nFileXSize - l_band->GetXSize();
     303          45 :         nFileXSize -= nRightFill;
     304             :     }
     305             : 
     306          45 :     if (nFileYOff + nFileYSize > l_band->GetYSize())
     307             :     {
     308          36 :         nBottomFill = nFileYOff + nFileYSize - l_band->GetYSize();
     309          36 :         nFileYSize -= nBottomFill;
     310             :     }
     311             : 
     312             :     /* -------------------------------------------------------------------- */
     313             :     /*      Load the data.                                                  */
     314             :     /* -------------------------------------------------------------------- */
     315             :     {
     316             :         GDALRasterIOExtraArg sExtraArgs;
     317          45 :         INIT_RASTERIO_EXTRA_ARG(sExtraArgs);
     318             :         const bool bIsComplex =
     319          45 :             CPL_TO_BOOL(GDALDataTypeIsComplex(eOperDataType));
     320          45 :         const CPLErr eErr = VRTComplexSource::RasterIOInternal<float>(
     321             :             l_band, eVRTBandDataType, nFileXOff, nFileYOff, nFileXSize,
     322             :             nFileYSize,
     323          45 :             pabyWorkData + nLineOffset * nTopFill + nPixelOffset * nLeftFill,
     324             :             nFileXSize, nFileYSize, eOperDataType, nPixelOffset, nLineOffset,
     325             :             &sExtraArgs, bIsComplex ? GDT_CFloat32 : GDT_Float32,
     326             :             oWorkingState);
     327             : 
     328          45 :         if (eErr != CE_None)
     329             :         {
     330           0 :             VSIFree(pabyWorkData);
     331           0 :             VSIFree(pabyOutData);
     332             : 
     333           0 :             return eErr;
     334             :         }
     335             :     }
     336             : 
     337             :     /* -------------------------------------------------------------------- */
     338             :     /*      Fill in missing areas.  Note that we replicate the edge         */
     339             :     /*      valid values out.  We don't using "mirroring" which might be    */
     340             :     /*      more suitable for some times of filters.  We also don't mark    */
     341             :     /*      these pixels as "nodata" though perhaps we should.              */
     342             :     /* -------------------------------------------------------------------- */
     343          45 :     if (nLeftFill != 0 || nRightFill != 0)
     344             :     {
     345        2212 :         for (int i = nTopFill; i < nExtraYSize - nBottomFill; i++)
     346             :         {
     347        2167 :             if (nLeftFill != 0)
     348        2167 :                 GDALCopyWords(
     349        2167 :                     pabyWorkData + nPixelOffset * nLeftFill + i * nLineOffset,
     350        2167 :                     eOperDataType, 0, pabyWorkData + i * nLineOffset,
     351             :                     eOperDataType, static_cast<int>(nPixelOffset), nLeftFill);
     352             : 
     353        2167 :             if (nRightFill != 0)
     354        2167 :                 GDALCopyWords(pabyWorkData + i * nLineOffset +
     355        2167 :                                   nPixelOffset * (nExtraXSize - nRightFill - 1),
     356             :                               eOperDataType, 0,
     357        2167 :                               pabyWorkData + i * nLineOffset +
     358        2167 :                                   nPixelOffset * (nExtraXSize - nRightFill),
     359             :                               eOperDataType, static_cast<int>(nPixelOffset),
     360             :                               nRightFill);
     361             :         }
     362             :     }
     363             : 
     364          88 :     for (int i = 0; i < nTopFill; i++)
     365             :     {
     366          43 :         memcpy(pabyWorkData + i * nLineOffset,
     367          43 :                pabyWorkData + nTopFill * nLineOffset, nLineOffset);
     368             :     }
     369             : 
     370          88 :     for (int i = nExtraYSize - nBottomFill; i < nExtraYSize; i++)
     371             :     {
     372          43 :         memcpy(pabyWorkData + i * nLineOffset,
     373          43 :                pabyWorkData + (nExtraYSize - nBottomFill - 1) * nLineOffset,
     374             :                nLineOffset);
     375             :     }
     376             : 
     377             :     /* -------------------------------------------------------------------- */
     378             :     /*      Filter the data.                                                */
     379             :     /* -------------------------------------------------------------------- */
     380          90 :     const CPLErr eErr = FilterData(nExtraXSize, nExtraYSize, eOperDataType,
     381          45 :                                    pabyWorkData, pabyOutData);
     382             : 
     383          45 :     VSIFree(pabyWorkData);
     384          45 :     if (eErr != CE_None)
     385             :     {
     386           0 :         VSIFree(pabyOutData);
     387             : 
     388           0 :         return eErr;
     389             :     }
     390             : 
     391             :     /* -------------------------------------------------------------------- */
     392             :     /*      Copy from work buffer to target buffer.                         */
     393             :     /* -------------------------------------------------------------------- */
     394          45 :     GByte *pabySrcRow =
     395          45 :         pabyOutData + (nLineOffset + nPixelOffset) * m_nExtraEdgePixels;
     396          45 :     GByte *pabyDstRow = reinterpret_cast<GByte *>(pData);
     397             : 
     398        2194 :     for (int i = 0; i < nOutYSize;
     399        2149 :          i++, pabySrcRow += nLineOffset, pabyDstRow += nLineSpace)
     400             :     {
     401        2149 :         GDALCopyWords(pabySrcRow, eOperDataType, static_cast<int>(nPixelOffset),
     402             :                       pabyDstRow, eBufType, static_cast<int>(nPixelSpace),
     403             :                       nOutXSize);
     404             :     }
     405             : 
     406          45 :     VSIFree(pabyOutData);
     407             : 
     408          45 :     return CE_None;
     409             : }
     410             : 
     411             : /************************************************************************/
     412             : /* ==================================================================== */
     413             : /*                       VRTKernelFilteredSource                        */
     414             : /* ==================================================================== */
     415             : /************************************************************************/
     416             : 
     417             : /************************************************************************/
     418             : /*                      VRTKernelFilteredSource()                       */
     419             : /************************************************************************/
     420             : 
     421          43 : VRTKernelFilteredSource::VRTKernelFilteredSource()
     422             : {
     423          43 :     GDALDataType aeSupTypes[] = {GDT_Float32};
     424          43 :     SetFilteringDataTypesSupported(1, aeSupTypes);
     425          43 : }
     426             : 
     427             : /************************************************************************/
     428             : /*                              GetType()                               */
     429             : /************************************************************************/
     430             : 
     431          86 : const char *VRTKernelFilteredSource::GetType() const
     432             : {
     433             :     static const char *TYPE = "KernelFilteredSource";
     434          86 :     return TYPE;
     435             : }
     436             : 
     437             : /************************************************************************/
     438             : /*                           SetNormalized()                            */
     439             : /************************************************************************/
     440             : 
     441          41 : void VRTKernelFilteredSource::SetNormalized(bool bNormalizedIn)
     442             : 
     443             : {
     444          41 :     m_bNormalized = bNormalizedIn;
     445          41 : }
     446             : 
     447             : /************************************************************************/
     448             : /*                             SetKernel()                              */
     449             : /************************************************************************/
     450             : 
     451             : CPLErr
     452          41 : VRTKernelFilteredSource::SetKernel(int nNewKernelSize, bool bSeparable,
     453             :                                    const std::vector<double> &adfNewCoefs)
     454             : 
     455             : {
     456          41 :     if (nNewKernelSize < 1 || (nNewKernelSize % 2) != 1)
     457             :     {
     458           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     459             :                  "Illegal filtering kernel size %d, "
     460             :                  "must be odd positive number.",
     461             :                  nNewKernelSize);
     462           0 :         return CE_Failure;
     463             :     }
     464          41 :     if (adfNewCoefs.size() !=
     465          41 :         static_cast<size_t>(nNewKernelSize) * (bSeparable ? 1 : nNewKernelSize))
     466             :     {
     467           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     468             :                  "adfNewCoefs[] is not of expected size");
     469           0 :         return CE_Failure;
     470             :     }
     471             : 
     472          41 :     m_nKernelSize = nNewKernelSize;
     473          41 :     m_bSeparable = bSeparable;
     474          41 :     m_adfKernelCoefs = adfNewCoefs;
     475             : 
     476          41 :     SetExtraEdgePixels((nNewKernelSize - 1) / 2);
     477             : 
     478          41 :     return CE_None;
     479             : }
     480             : 
     481             : /************************************************************************/
     482             : /*                             FilterData()                             */
     483             : /************************************************************************/
     484             : 
     485          45 : CPLErr VRTKernelFilteredSource::FilterData(int nXSize, int nYSize,
     486             :                                            GDALDataType eType,
     487             :                                            GByte *pabySrcData,
     488             :                                            GByte *pabyDstData)
     489             : 
     490             : {
     491             :     /* -------------------------------------------------------------------- */
     492             :     /*      Validate data type.                                             */
     493             :     /* -------------------------------------------------------------------- */
     494          45 :     if (eType != GDT_Float32)
     495             :     {
     496           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     497             :                  "Unsupported data type (%s) in "
     498             :                  "VRTKernelFilteredSource::FilterData()",
     499             :                  GDALGetDataTypeName(eType));
     500           0 :         return CE_Failure;
     501             :     }
     502             : 
     503          45 :     CPLAssert(m_nExtraEdgePixels * 2 + 1 == m_nKernelSize ||
     504             :               (m_nKernelSize == 0 && m_nExtraEdgePixels == 0));
     505             : 
     506          45 :     const bool bMin = m_function == "min";
     507          45 :     const bool bMax = m_function == "max";
     508          45 :     const bool bStdDev = m_function == "stddev";
     509          45 :     const bool bMedian = m_function == "median";
     510          45 :     const bool bMode = m_function == "mode";
     511             : 
     512             :     /* -------------------------------------------------------------------- */
     513             :     /*      Float32 case.                                                   */
     514             :     /* -------------------------------------------------------------------- */
     515             :     // if( eType == GDT_Float32 )
     516             :     {
     517          45 :         float *pafSrcData = reinterpret_cast<float *>(pabySrcData);
     518          45 :         float *pafDstData = reinterpret_cast<float *>(pabyDstData);
     519             : 
     520          45 :         int bHasNoData = FALSE;
     521          45 :         auto l_poBand = GetRasterBand();
     522          45 :         if (!l_poBand)
     523           0 :             return CE_Failure;
     524             :         const float fNoData =
     525          45 :             static_cast<float>(l_poBand->GetNoDataValue(&bHasNoData));
     526             : 
     527          45 :         const int nAxisCount = m_bSeparable ? 2 : 1;
     528             : 
     529          91 :         for (int nAxis = 0; nAxis < nAxisCount; ++nAxis)
     530             :         {
     531          46 :             const int nISize = nAxis == 0 ? nYSize : nXSize;
     532          46 :             const int nJSize = nAxis == 0 ? nXSize : nYSize;
     533          46 :             const int nIStride = nAxis == 0 ? nXSize : 1;
     534          46 :             const int nJStride = nAxis == 0 ? 1 : nXSize;
     535             : 
     536          46 :             const int nIMin = m_nExtraEdgePixels;
     537          46 :             const int nIMax = nISize - m_nExtraEdgePixels;
     538          46 :             const int nJMin = (m_bSeparable ? 0 : m_nExtraEdgePixels);
     539          46 :             const int nJMax = nJSize - (m_bSeparable ? 0 : m_nExtraEdgePixels);
     540             : 
     541        7069 :             for (GPtrDiff_t iJ = nJMin; iJ < nJMax; ++iJ)
     542             :             {
     543        7023 :                 if (nAxis == 1)
     544          62 :                     memcpy(pafSrcData + iJ * nJStride,
     545          62 :                            pafDstData + iJ * nJStride, sizeof(float) * nXSize);
     546             : 
     547      768540 :                 for (int iI = nIMin; iI < nIMax; ++iI)
     548             :                 {
     549      761517 :                     const GPtrDiff_t iIndex =
     550      761517 :                         static_cast<GPtrDiff_t>(iI) * nIStride + iJ * nJStride;
     551             : 
     552      761517 :                     if (bHasNoData && pafSrcData[iIndex] == fNoData)
     553             :                     {
     554         204 :                         pafDstData[iIndex] = fNoData;
     555         204 :                         continue;
     556             :                     }
     557             : 
     558      761313 :                     double dfSum = 0.0, dfKernSum = 0.0;
     559      761313 :                     size_t nValidCount = 0;
     560             :                     double dfRes =
     561     1522620 :                         bMin   ? std::numeric_limits<double>::infinity()
     562      761304 :                         : bMax ? -std::numeric_limits<double>::infinity()
     563      761313 :                                : 0.0;
     564      761313 :                     double dfMean = 0.0;
     565      761313 :                     double dfM2 = 0.0;
     566     1522630 :                     std::vector<double> adfVals;
     567             : 
     568             :                     struct CompareNaNAware
     569             :                     {
     570         265 :                         bool operator()(double lhs, double rhs) const
     571             :                         {
     572         265 :                             return (std::isnan(lhs) && !std::isnan(rhs)) ||
     573         265 :                                    lhs < rhs;
     574             :                         }
     575             :                     };
     576             : 
     577      761313 :                     std::map<double, size_t, CompareNaNAware> mapValToCount;
     578      761313 :                     size_t maxCount = 0;
     579             : 
     580     3108850 :                     for (GPtrDiff_t iII = -m_nExtraEdgePixels, iK = 0;
     581     3108850 :                          iII <= m_nExtraEdgePixels; ++iII)
     582             :                     {
     583     9236960 :                         for (GPtrDiff_t iJJ =
     584     2347540 :                                  (m_bSeparable ? 0 : -m_nExtraEdgePixels);
     585     9236960 :                              iJJ <= (m_bSeparable ? 0 : m_nExtraEdgePixels);
     586             :                              ++iJJ, ++iK)
     587             :                         {
     588     6889420 :                             const float *pfData = pafSrcData + iIndex +
     589     6889420 :                                                   iII * nIStride +
     590     6889420 :                                                   iJJ * nJStride;
     591     6889420 :                             if (bHasNoData &&
     592     6889420 :                                 (*pfData == fNoData || std::isnan(*pfData)))
     593             :                             {
     594      977782 :                                 continue;
     595             :                             }
     596     6888560 :                             if (m_adfKernelCoefs[iK] == 0.0)
     597             :                             {
     598      976926 :                                 continue;
     599             :                             }
     600    11823300 :                             const double dfVal = static_cast<double>(*pfData) *
     601     5911640 :                                                  m_adfKernelCoefs[iK];
     602     5911640 :                             ++nValidCount;
     603             : 
     604     5911640 :                             if (bMax)
     605             :                             {
     606         282 :                                 if (dfVal > dfRes)
     607         175 :                                     dfRes = dfVal;
     608             :                             }
     609     5911350 :                             else if (bMin)
     610             :                             {
     611          81 :                                 if (dfVal < dfRes)
     612           9 :                                     dfRes = dfVal;
     613             :                             }
     614     5911270 :                             else if (bStdDev)
     615             :                             {
     616          81 :                                 const double dfDelta = dfVal - dfMean;
     617          81 :                                 dfMean += dfDelta / nValidCount;
     618          81 :                                 dfM2 += dfDelta * (dfVal - dfMean);
     619             :                             }
     620     5911190 :                             else if (bMedian)
     621             :                             {
     622         148 :                                 adfVals.push_back(dfVal);
     623             :                             }
     624     5911040 :                             else if (bMode)
     625             :                             {
     626          81 :                                 const size_t nCountVal = ++mapValToCount[dfVal];
     627          81 :                                 if (nCountVal > maxCount)
     628             :                                 {
     629          26 :                                     maxCount = nCountVal;
     630          26 :                                     dfRes = dfVal;
     631             :                                 }
     632             :                             }
     633             :                             else
     634             :                             {
     635     5910960 :                                 dfSum += dfVal;
     636     5910960 :                                 dfKernSum += m_adfKernelCoefs[iK];
     637             :                             }
     638             :                         }
     639             :                     }
     640             : 
     641             :                     float fResult;
     642      761313 :                     if (bMax || bMin || bMode)
     643             :                     {
     644          51 :                         fResult = nValidCount  ? static_cast<float>(dfRes)
     645           0 :                                   : bHasNoData ? fNoData
     646           0 :                                                : 0.0f;
     647             :                     }
     648      761262 :                     else if (bStdDev)
     649             :                     {
     650           9 :                         fResult =
     651             :                             nValidCount
     652           9 :                                 ? static_cast<float>(sqrt(
     653           9 :                                       dfM2 / static_cast<double>(nValidCount)))
     654           0 :                             : bHasNoData ? fNoData
     655           0 :                                          : 0.0f;
     656             :                     }
     657      761253 :                     else if (bMedian)
     658             :                     {
     659          17 :                         if (!adfVals.empty())
     660             :                         {
     661          17 :                             if ((adfVals.size() % 2) == 1)
     662             :                             {
     663          32 :                                 std::nth_element(adfVals.begin(),
     664           0 :                                                  adfVals.begin() +
     665          16 :                                                      adfVals.size() / 2,
     666             :                                                  adfVals.end());
     667          16 :                                 dfRes = adfVals[adfVals.size() / 2];
     668             :                             }
     669             :                             else
     670             :                             {
     671           2 :                                 std::nth_element(adfVals.begin(),
     672           0 :                                                  adfVals.begin() +
     673           1 :                                                      adfVals.size() / 2 - 1,
     674             :                                                  adfVals.end());
     675           1 :                                 dfRes = adfVals[adfVals.size() / 2 - 1];
     676           2 :                                 std::nth_element(adfVals.begin(),
     677           0 :                                                  adfVals.begin() +
     678           1 :                                                      adfVals.size() / 2,
     679             :                                                  adfVals.end());
     680           1 :                                 dfRes =
     681           1 :                                     (dfRes + adfVals[adfVals.size() / 2]) / 2;
     682             :                             }
     683          17 :                             fResult = static_cast<float>(dfRes);
     684             :                         }
     685             :                         else
     686           0 :                             fResult = bHasNoData ? fNoData : 0.0f;
     687             :                     }
     688      761236 :                     else if (!m_bNormalized)
     689      512027 :                         fResult = static_cast<float>(dfSum);
     690      249209 :                     else if (nValidCount == 0 || dfKernSum == 0.0)
     691           0 :                         fResult = bHasNoData ? fNoData : 0.0f;
     692             :                     else
     693      249209 :                         fResult = static_cast<float>(dfSum / dfKernSum);
     694      761313 :                     pafDstData[iIndex] = fResult;
     695             :                 }
     696             :             }
     697             :         }
     698             :     }
     699             : 
     700          45 :     return CE_None;
     701             : }
     702             : 
     703             : /************************************************************************/
     704             : /*                              XMLInit()                               */
     705             : /************************************************************************/
     706             : 
     707             : CPLErr
     708          14 : VRTKernelFilteredSource::XMLInit(const CPLXMLNode *psTree,
     709             :                                  const char *pszVRTPath,
     710             :                                  VRTMapSharedResources &oMapSharedSources)
     711             : 
     712             : {
     713             :     {
     714             :         const CPLErr eErr =
     715          14 :             VRTFilteredSource::XMLInit(psTree, pszVRTPath, oMapSharedSources);
     716          14 :         if (eErr != CE_None)
     717           0 :             return eErr;
     718             :     }
     719             : 
     720          14 :     const int nNewKernelSize = atoi(CPLGetXMLValue(psTree, "Kernel.Size", "0"));
     721             : 
     722          14 :     if (nNewKernelSize == 0)
     723           0 :         return CE_None;
     724             :     // To prevent a integer overflow when computing nNewKernelSize *
     725             :     // nNewKernelSize
     726          27 :     if (nNewKernelSize < 0 ||
     727          13 :         nNewKernelSize > static_cast<int>(std::sqrt(static_cast<double>(
     728          13 :                              std::numeric_limits<int>::max()))))
     729             :     {
     730           2 :         CPLError(CE_Failure, CPLE_IllegalArg,
     731             :                  "Invalid value for kernel size: %d", nNewKernelSize);
     732           2 :         return CE_Failure;
     733             :     }
     734             : 
     735             :     const CPLStringList aosCoefItems(
     736          24 :         CSLTokenizeString(CPLGetXMLValue(psTree, "Kernel.Coefs", "")));
     737             : 
     738          12 :     const int nCoefs = aosCoefItems.size();
     739             : 
     740          12 :     const bool bSquare = nCoefs == nNewKernelSize * nNewKernelSize;
     741          12 :     const bool bSeparable = nCoefs == nNewKernelSize && nCoefs != 1;
     742             : 
     743          12 :     if (!bSquare && !bSeparable)
     744             :     {
     745           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     746             :                  "Got wrong number of filter kernel coefficients (%s).  "
     747             :                  "Expected %d or %d, got %d.",
     748             :                  CPLGetXMLValue(psTree, "Kernel.Coefs", ""),
     749             :                  nNewKernelSize * nNewKernelSize, nNewKernelSize, nCoefs);
     750           0 :         return CE_Failure;
     751             :     }
     752             : 
     753          24 :     std::vector<double> adfNewCoefs;
     754          12 :     adfNewCoefs.reserve(nCoefs);
     755         122 :     for (int i = 0; i < nCoefs; i++)
     756         110 :         adfNewCoefs.push_back(CPLAtof(aosCoefItems[i]));
     757             : 
     758          12 :     const CPLErr eErr = SetKernel(nNewKernelSize, bSeparable, adfNewCoefs);
     759          12 :     if (eErr == CE_None)
     760             :     {
     761          12 :         SetNormalized(atoi(CPLGetXMLValue(psTree, "Kernel.normalized", "0")) !=
     762             :                       0);
     763             :     }
     764             : 
     765          12 :     const char *pszFunction = CPLGetXMLValue(psTree, "Function", nullptr);
     766          12 :     if (pszFunction)
     767             :     {
     768           0 :         if (!EQUAL(pszFunction, "max") && !EQUAL(pszFunction, "min") &&
     769           0 :             !EQUAL(pszFunction, "median") && !EQUAL(pszFunction, "mode") &&
     770           0 :             !EQUAL(pszFunction, "stddev"))
     771             :         {
     772           0 :             CPLError(CE_Failure, CPLE_IllegalArg, "Unsupported function: %s",
     773             :                      pszFunction);
     774           0 :             return CE_Failure;
     775             :         }
     776           0 :         SetFunction(pszFunction);
     777             :     }
     778             : 
     779          12 :     return eErr;
     780             : }
     781             : 
     782             : /************************************************************************/
     783             : /*                           SerializeToXML()                           */
     784             : /************************************************************************/
     785             : 
     786           2 : CPLXMLNode *VRTKernelFilteredSource::SerializeToXML(const char *pszVRTPath)
     787             : 
     788             : {
     789           2 :     CPLXMLNode *psSrc = VRTFilteredSource::SerializeToXML(pszVRTPath);
     790             : 
     791           2 :     if (psSrc == nullptr)
     792           0 :         return nullptr;
     793             : 
     794           2 :     CPLFree(psSrc->pszValue);
     795           2 :     psSrc->pszValue = CPLStrdup("KernelFilteredSource");
     796             : 
     797           2 :     if (m_nKernelSize == 0)
     798           0 :         return psSrc;
     799             : 
     800           2 :     CPLXMLNode *psKernel = CPLCreateXMLNode(psSrc, CXT_Element, "Kernel");
     801             : 
     802           2 :     CPLCreateXMLNode(CPLCreateXMLNode(psKernel, CXT_Attribute, "normalized"),
     803           2 :                      CXT_Text, m_bNormalized ? "1" : "0");
     804             : 
     805           2 :     std::string osCoefs;
     806          14 :     for (auto dfVal : m_adfKernelCoefs)
     807             :     {
     808          12 :         if (!osCoefs.empty())
     809          10 :             osCoefs += ' ';
     810          12 :         osCoefs += CPLSPrintf("%.8g", dfVal);
     811             :     }
     812             : 
     813           2 :     CPLSetXMLValue(psKernel, "Size", CPLSPrintf("%d", m_nKernelSize));
     814           2 :     CPLSetXMLValue(psKernel, "Coefs", osCoefs.c_str());
     815             : 
     816           2 :     if (!m_function.empty())
     817           0 :         CPLCreateXMLElementAndValue(psSrc, "Function", m_function.c_str());
     818             : 
     819           2 :     return psSrc;
     820             : }
     821             : 
     822             : /************************************************************************/
     823             : /*                       VRTParseFilterSources()                        */
     824             : /************************************************************************/
     825             : 
     826          14 : VRTSource *VRTParseFilterSources(const CPLXMLNode *psChild,
     827             :                                  const char *pszVRTPath,
     828             :                                  VRTMapSharedResources &oMapSharedSources)
     829             : 
     830             : {
     831          14 :     if (EQUAL(psChild->pszValue, "KernelFilteredSource"))
     832             :     {
     833          14 :         VRTSource *poSrc = new VRTKernelFilteredSource();
     834          14 :         if (poSrc->XMLInit(psChild, pszVRTPath, oMapSharedSources) == CE_None)
     835          12 :             return poSrc;
     836             : 
     837           2 :         delete poSrc;
     838             :     }
     839             : 
     840           2 :     return nullptr;
     841             : }
     842             : 
     843             : /*! @endcond */