LCOV - code coverage report
Current view: top level - alg - gdalpansharpen.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 609 749 81.3 %
Date: 2025-10-22 13:51:22 Functions: 40 63 63.5 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL Pansharpening module
       4             :  * Purpose:  Implementation of pansharpening.
       5             :  * Author:   Even Rouault <even.rouault at spatialys.com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2015, Even Rouault <even.rouault at spatialys.com>
       9             :  * Copyright (c) 2015, Airbus DS Geo SA (weighted Brovey algorithm)
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include "cpl_port.h"
      15             : #include "cpl_worker_thread_pool.h"
      16             : #include "gdalpansharpen.h"
      17             : 
      18             : #include <algorithm>
      19             : #include <array>
      20             : #include <cstddef>
      21             : #include <cstdio>
      22             : #include <cstdlib>
      23             : #include <cstring>
      24             : #include <limits>
      25             : #include <new>
      26             : 
      27             : #include "cpl_conv.h"
      28             : #include "cpl_error.h"
      29             : #include "cpl_float.h"
      30             : #include "cpl_multiproc.h"
      31             : #include "cpl_vsi.h"
      32             : #include "../frmts/mem/memdataset.h"
      33             : #include "../frmts/vrt/vrtdataset.h"
      34             : #include "gdal_priv.h"
      35             : #include "gdal_priv_templates.hpp"
      36             : // #include "gdalsse_priv.h"
      37             : 
      38             : // Limit types to practical use cases.
      39             : #define LIMIT_TYPES 1
      40             : 
      41             : /************************************************************************/
      42             : /*                     GDALCreatePansharpenOptions()                    */
      43             : /************************************************************************/
      44             : 
      45             : /** Create pansharpening options.
      46             :  *
      47             :  * @return a newly allocated pansharpening option structure that must be freed
      48             :  * with GDALDestroyPansharpenOptions().
      49             :  *
      50             :  */
      51             : 
      52         139 : GDALPansharpenOptions *GDALCreatePansharpenOptions()
      53             : {
      54             :     GDALPansharpenOptions *psOptions = static_cast<GDALPansharpenOptions *>(
      55         139 :         CPLCalloc(1, sizeof(GDALPansharpenOptions)));
      56         139 :     psOptions->ePansharpenAlg = GDAL_PSH_WEIGHTED_BROVEY;
      57         139 :     psOptions->eResampleAlg = GRIORA_Cubic;
      58         139 :     return psOptions;
      59             : }
      60             : 
      61             : /************************************************************************/
      62             : /*                     GDALDestroyPansharpenOptions()                   */
      63             : /************************************************************************/
      64             : 
      65             : /** Destroy pansharpening options.
      66             :  *
      67             :  * @param psOptions a pansharpening option structure allocated with
      68             :  * GDALCreatePansharpenOptions()
      69             :  *
      70             :  */
      71             : 
      72         140 : void GDALDestroyPansharpenOptions(GDALPansharpenOptions *psOptions)
      73             : {
      74         140 :     if (psOptions == nullptr)
      75           1 :         return;
      76         139 :     CPLFree(psOptions->padfWeights);
      77         139 :     CPLFree(psOptions->pahInputSpectralBands);
      78         139 :     CPLFree(psOptions->panOutPansharpenedBands);
      79         139 :     CPLFree(psOptions);
      80             : }
      81             : 
      82             : /************************************************************************/
      83             : /*                      GDALClonePansharpenOptions()                    */
      84             : /************************************************************************/
      85             : 
      86             : /** Clone pansharpening options.
      87             :  *
      88             :  * @param psOptions a pansharpening option structure allocated with
      89             :  * GDALCreatePansharpenOptions()
      90             :  * @return a newly allocated pansharpening option structure that must be freed
      91             :  * with GDALDestroyPansharpenOptions().
      92             :  *
      93             :  */
      94             : 
      95             : GDALPansharpenOptions *
      96          72 : GDALClonePansharpenOptions(const GDALPansharpenOptions *psOptions)
      97             : {
      98          72 :     GDALPansharpenOptions *psNewOptions = GDALCreatePansharpenOptions();
      99          72 :     psNewOptions->ePansharpenAlg = psOptions->ePansharpenAlg;
     100          72 :     psNewOptions->eResampleAlg = psOptions->eResampleAlg;
     101          72 :     psNewOptions->nBitDepth = psOptions->nBitDepth;
     102          72 :     psNewOptions->nWeightCount = psOptions->nWeightCount;
     103          72 :     if (psOptions->padfWeights)
     104             :     {
     105          72 :         psNewOptions->padfWeights = static_cast<double *>(
     106          72 :             CPLMalloc(sizeof(double) * psOptions->nWeightCount));
     107          72 :         memcpy(psNewOptions->padfWeights, psOptions->padfWeights,
     108          72 :                sizeof(double) * psOptions->nWeightCount);
     109             :     }
     110          72 :     psNewOptions->hPanchroBand = psOptions->hPanchroBand;
     111          72 :     psNewOptions->nInputSpectralBands = psOptions->nInputSpectralBands;
     112          72 :     if (psOptions->pahInputSpectralBands)
     113             :     {
     114          72 :         const size_t nSize =
     115          72 :             sizeof(GDALRasterBandH) * psOptions->nInputSpectralBands;
     116          72 :         psNewOptions->pahInputSpectralBands =
     117          72 :             static_cast<GDALRasterBandH *>(CPLMalloc(nSize));
     118          72 :         memcpy(psNewOptions->pahInputSpectralBands,
     119          72 :                psOptions->pahInputSpectralBands, nSize);
     120             :     }
     121          72 :     psNewOptions->nOutPansharpenedBands = psOptions->nOutPansharpenedBands;
     122          72 :     if (psOptions->panOutPansharpenedBands)
     123             :     {
     124          71 :         psNewOptions->panOutPansharpenedBands = static_cast<int *>(
     125          71 :             CPLMalloc(sizeof(int) * psOptions->nOutPansharpenedBands));
     126          71 :         memcpy(psNewOptions->panOutPansharpenedBands,
     127          71 :                psOptions->panOutPansharpenedBands,
     128          71 :                sizeof(int) * psOptions->nOutPansharpenedBands);
     129             :     }
     130          72 :     psNewOptions->bHasNoData = psOptions->bHasNoData;
     131          72 :     psNewOptions->dfNoData = psOptions->dfNoData;
     132          72 :     psNewOptions->nThreads = psOptions->nThreads;
     133          72 :     return psNewOptions;
     134             : }
     135             : 
     136             : /************************************************************************/
     137             : /*                        GDALPansharpenOperation()                     */
     138             : /************************************************************************/
     139             : 
     140             : /** Pansharpening operation constructor.
     141             :  *
     142             :  * The object is ready to be used after Initialize() has been called.
     143             :  */
     144             : GDALPansharpenOperation::GDALPansharpenOperation() = default;
     145             : 
     146             : /************************************************************************/
     147             : /*                       ~GDALPansharpenOperation()                     */
     148             : /************************************************************************/
     149             : 
     150             : /** Pansharpening operation destructor.
     151             :  */
     152             : 
     153          70 : GDALPansharpenOperation::~GDALPansharpenOperation()
     154             : {
     155          70 :     GDALDestroyPansharpenOptions(psOptions);
     156          79 :     for (size_t i = 0; i < aVDS.size(); i++)
     157           9 :         delete aVDS[i];
     158          70 :     delete poThreadPool;
     159          70 : }
     160             : 
     161             : /************************************************************************/
     162             : /*                              Initialize()                            */
     163             : /************************************************************************/
     164             : 
     165             : /** Initialize the pansharpening operation.
     166             :  *
     167             :  * @param psOptionsIn pansharpening options. Must not be NULL.
     168             :  *
     169             :  * @return CE_None in case of success, CE_Failure in case of failure.
     170             :  */
     171             : CPLErr
     172          70 : GDALPansharpenOperation::Initialize(const GDALPansharpenOptions *psOptionsIn)
     173             : {
     174          70 :     if (psOptionsIn->hPanchroBand == nullptr)
     175             :     {
     176           0 :         CPLError(CE_Failure, CPLE_AppDefined, "hPanchroBand not set");
     177           0 :         return CE_Failure;
     178             :     }
     179          70 :     if (psOptionsIn->nInputSpectralBands <= 0)
     180             :     {
     181           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     182             :                  "No input spectral bands defined");
     183           0 :         return CE_Failure;
     184             :     }
     185          70 :     if (psOptionsIn->padfWeights == nullptr ||
     186          70 :         psOptionsIn->nWeightCount != psOptionsIn->nInputSpectralBands)
     187             :     {
     188           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     189             :                  "No weights defined, or not the same number as input "
     190             :                  "spectral bands");
     191           0 :         return CE_Failure;
     192             :     }
     193             : 
     194          70 :     auto poPanchroBand = GDALRasterBand::FromHandle(psOptionsIn->hPanchroBand);
     195          70 :     auto poPanchroDS = poPanchroBand->GetDataset();
     196          70 :     if (poPanchroDS == nullptr)
     197             :     {
     198           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     199             :                  "Cannot retrieve dataset associated with hPanchroBand");
     200           0 :         return CE_Failure;
     201             :     }
     202             :     // Make sure that the band is really a first level child of the owning dataset
     203          70 :     if (poPanchroDS->GetRasterBand(poPanchroBand->GetBand()) != poPanchroBand)
     204             :     {
     205           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     206             :                  "poPanchroDS->GetRasterBand(poPanchroBand->GetBand()) != "
     207             :                  "poPanchroBand");
     208           0 :         return CE_Failure;
     209             :     }
     210          70 :     GDALGeoTransform panchroGT;
     211          70 :     if (poPanchroDS->GetGeoTransform(panchroGT) != CE_None)
     212             :     {
     213           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     214             :                  "Panchromatic band has no associated geotransform");
     215           0 :         return CE_Failure;
     216             :     }
     217             : 
     218          70 :     GDALRasterBandH hRefBand = psOptionsIn->pahInputSpectralBands[0];
     219          70 :     auto poRefBand = GDALRasterBand::FromHandle(hRefBand);
     220          70 :     auto poRefBandDS = poRefBand->GetDataset();
     221          70 :     if (poRefBandDS == nullptr)
     222             :     {
     223           0 :         CPLError(
     224             :             CE_Failure, CPLE_AppDefined,
     225             :             "Cannot retrieve dataset associated with ahInputSpectralBands[0]");
     226           0 :         return CE_Failure;
     227             :     }
     228             :     // Make sure that the band is really a first level child of the owning dataset
     229          70 :     if (poRefBandDS->GetRasterBand(poRefBand->GetBand()) != poRefBand)
     230             :     {
     231           0 :         CPLError(
     232             :             CE_Failure, CPLE_AppDefined,
     233             :             "poRefBandDS->GetRasterBand(poRefBand->GetBand()) != poRefBand");
     234           0 :         return CE_Failure;
     235             :     }
     236             : 
     237          70 :     GDALGeoTransform refMSGT;
     238          70 :     if (poRefBandDS->GetGeoTransform(refMSGT) != CE_None)
     239             :     {
     240           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     241             :                  "ahInputSpectralBands[0] band has no associated geotransform");
     242           0 :         return CE_Failure;
     243             :     }
     244             : 
     245          70 :     GDALGeoTransform invMSGT;
     246          70 :     if (!refMSGT.GetInverse(invMSGT))
     247             :     {
     248           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     249             :                  "ahInputSpectralBands[0] geotransform is not invertible");
     250           0 :         return CE_Failure;
     251             :     }
     252             : 
     253             :     // Do InvMSGT * PanchroGT multiplication
     254          70 :     m_panToMSGT[1] = invMSGT[1] * panchroGT[1] + invMSGT[2] * panchroGT[4];
     255          70 :     m_panToMSGT[2] = invMSGT[1] * panchroGT[2] + invMSGT[2] * panchroGT[5];
     256         140 :     m_panToMSGT[0] =
     257          70 :         invMSGT[1] * panchroGT[0] + invMSGT[2] * panchroGT[3] + invMSGT[0];
     258          70 :     m_panToMSGT[4] = invMSGT[4] * panchroGT[1] + invMSGT[5] * panchroGT[4];
     259          70 :     m_panToMSGT[5] = invMSGT[4] * panchroGT[2] + invMSGT[5] * panchroGT[5];
     260         140 :     m_panToMSGT[3] =
     261          70 :         invMSGT[4] * panchroGT[0] + invMSGT[5] * panchroGT[3] + invMSGT[3];
     262             : #if 0
     263             :     CPLDebug("GDAL", "m_panToMSGT[] = %g %g %g %g %g %g",
     264             :            m_panToMSGT[0], m_panToMSGT[1], m_panToMSGT[2],
     265             :            m_panToMSGT[3], m_panToMSGT[4], m_panToMSGT[5]);
     266             : #endif
     267          70 :     if (std::fabs(m_panToMSGT[2]) > 1e-10 || std::fabs(m_panToMSGT[4]) > 1e-10)
     268             :     {
     269           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     270             :                  "Composition of panchromatic and multispectral geotransform "
     271             :                  "has rotational terms");
     272           0 :         return CE_Failure;
     273             :     }
     274             : 
     275          70 :     bool bSameDataset = psOptionsIn->nInputSpectralBands > 1;
     276          70 :     if (bSameDataset)
     277          55 :         anInputBands.push_back(GDALGetBandNumber(hRefBand));
     278         174 :     for (int i = 1; i < psOptionsIn->nInputSpectralBands; i++)
     279             :     {
     280         105 :         GDALRasterBandH hBand = psOptionsIn->pahInputSpectralBands[i];
     281         209 :         if (GDALGetRasterBandXSize(hBand) != GDALGetRasterBandXSize(hRefBand) ||
     282         104 :             GDALGetRasterBandYSize(hBand) != GDALGetRasterBandYSize(hRefBand))
     283             :         {
     284           1 :             CPLError(CE_Failure, CPLE_AppDefined,
     285             :                      "Dimensions of input spectral band %d different from "
     286             :                      "first spectral band",
     287             :                      i + 1);
     288           1 :             return CE_Failure;
     289             :         }
     290             : 
     291         104 :         auto poBand = GDALRasterBand::FromHandle(hBand);
     292         104 :         auto poBandDS = poBand->GetDataset();
     293         104 :         if (poBandDS == nullptr)
     294             :         {
     295           0 :             CPLError(CE_Failure, CPLE_AppDefined,
     296             :                      "Cannot retrieve dataset associated with "
     297             :                      "input spectral band %d",
     298             :                      i + 1);
     299           0 :             return CE_Failure;
     300             :         }
     301             :         // Make sure that the band is really a first level child of the owning dataset
     302         104 :         if (poBandDS->GetRasterBand(poBand->GetBand()) != poBand)
     303             :         {
     304           0 :             CPLError(CE_Failure, CPLE_AppDefined,
     305             :                      "poBandDS->GetRasterBand(poBand->GetBand()) != poBand");
     306           0 :             return CE_Failure;
     307             :         }
     308             : 
     309         104 :         GDALGeoTransform MSGT;
     310         104 :         if (poBandDS->GetGeoTransform(MSGT) != CE_None)
     311             :         {
     312           0 :             CPLError(
     313             :                 CE_Failure, CPLE_AppDefined,
     314             :                 "input spectral band %d band has no associated geotransform",
     315             :                 i + 1);
     316           0 :             return CE_Failure;
     317             :         }
     318         104 :         if (MSGT != refMSGT)
     319             :         {
     320           0 :             CPLError(CE_Failure, CPLE_AppDefined,
     321             :                      "input spectral band %d has a different "
     322             :                      "geotransform than the first spectral band",
     323             :                      i + 1);
     324           0 :             return CE_Failure;
     325             :         }
     326             : 
     327         104 :         if (bSameDataset)
     328             :         {
     329         101 :             if (GDALGetBandDataset(hBand) != GDALGetBandDataset(hRefBand))
     330             :             {
     331           5 :                 anInputBands.resize(0);
     332           5 :                 bSameDataset = false;
     333             :             }
     334             :             else
     335             :             {
     336          96 :                 anInputBands.push_back(GDALGetBandNumber(hBand));
     337             :             }
     338             :         }
     339             :     }
     340          69 :     if (psOptionsIn->nOutPansharpenedBands == 0)
     341             :     {
     342           1 :         CPLError(CE_Warning, CPLE_AppDefined,
     343             :                  "No output pansharpened band defined");
     344             :     }
     345         236 :     for (int i = 0; i < psOptionsIn->nOutPansharpenedBands; i++)
     346             :     {
     347         167 :         if (psOptionsIn->panOutPansharpenedBands[i] < 0 ||
     348         167 :             psOptionsIn->panOutPansharpenedBands[i] >=
     349         167 :                 psOptionsIn->nInputSpectralBands)
     350             :         {
     351           0 :             CPLError(CE_Failure, CPLE_AppDefined,
     352             :                      "Invalid value panOutPansharpenedBands[%d] = %d", i,
     353           0 :                      psOptionsIn->panOutPansharpenedBands[i]);
     354           0 :             return CE_Failure;
     355             :         }
     356             :     }
     357             : 
     358          69 :     GDALDataType eWorkDataType = poPanchroBand->GetRasterDataType();
     359          69 :     if (psOptionsIn->nBitDepth)
     360             :     {
     361          13 :         if (psOptionsIn->nBitDepth < 0 || psOptionsIn->nBitDepth > 31 ||
     362          13 :             (eWorkDataType == GDT_Byte && psOptionsIn->nBitDepth > 8) ||
     363           3 :             (eWorkDataType == GDT_UInt16 && psOptionsIn->nBitDepth > 16))
     364             :         {
     365           0 :             CPLError(CE_Failure, CPLE_AppDefined,
     366             :                      "Invalid value nBitDepth = %d for type %s",
     367           0 :                      psOptionsIn->nBitDepth,
     368             :                      GDALGetDataTypeName(eWorkDataType));
     369           0 :             return CE_Failure;
     370             :         }
     371             :     }
     372             : 
     373          69 :     psOptions = GDALClonePansharpenOptions(psOptionsIn);
     374          69 :     if (psOptions->nBitDepth == GDALGetDataTypeSizeBits(eWorkDataType))
     375           7 :         psOptions->nBitDepth = 0;
     376          69 :     if (psOptions->nBitDepth &&
     377           3 :         !(eWorkDataType == GDT_Byte || eWorkDataType == GDT_UInt16 ||
     378             :           eWorkDataType == GDT_UInt32 || eWorkDataType == GDT_UInt64))
     379             :     {
     380           0 :         CPLError(CE_Warning, CPLE_AppDefined,
     381           0 :                  "Ignoring nBitDepth = %d for type %s", psOptions->nBitDepth,
     382             :                  GDALGetDataTypeName(eWorkDataType));
     383           0 :         psOptions->nBitDepth = 0;
     384             :     }
     385             : 
     386             :     // Detect negative weights.
     387         242 :     for (int i = 0; i < psOptions->nInputSpectralBands; i++)
     388             :     {
     389         173 :         if (psOptions->padfWeights[i] < 0.0)
     390             :         {
     391           0 :             bPositiveWeights = FALSE;
     392           0 :             break;
     393             :         }
     394             :     }
     395             : 
     396         242 :     for (int i = 0; i < psOptions->nInputSpectralBands; i++)
     397             :     {
     398         346 :         aMSBands.push_back(
     399         173 :             GDALRasterBand::FromHandle(psOptions->pahInputSpectralBands[i]));
     400             :     }
     401             : 
     402          69 :     if (psOptions->bHasNoData)
     403             :     {
     404           9 :         bool bNeedToWrapInVRT = false;
     405          29 :         for (int i = 0; i < psOptions->nInputSpectralBands; i++)
     406             :         {
     407             :             GDALRasterBand *poBand =
     408          20 :                 GDALRasterBand::FromHandle(psOptions->pahInputSpectralBands[i]);
     409          20 :             int bHasNoData = FALSE;
     410          20 :             double dfNoData = poBand->GetNoDataValue(&bHasNoData);
     411          20 :             if (!bHasNoData || dfNoData != psOptions->dfNoData)
     412          17 :                 bNeedToWrapInVRT = true;
     413             :         }
     414             : 
     415           9 :         if (bNeedToWrapInVRT)
     416             :         {
     417             :             // Wrap spectral bands in a VRT if they don't have the nodata value.
     418           8 :             VRTDataset *poVDS = nullptr;
     419          25 :             for (int i = 0; i < psOptions->nInputSpectralBands; i++)
     420             :             {
     421          17 :                 GDALRasterBand *poSrcBand = aMSBands[i];
     422          17 :                 int iVRTBand = 1;
     423          17 :                 if (anInputBands.empty() || i == 0)
     424             :                 {
     425           9 :                     poVDS = new VRTDataset(poSrcBand->GetXSize(),
     426           9 :                                            poSrcBand->GetYSize());
     427           9 :                     aVDS.push_back(poVDS);
     428             :                 }
     429          17 :                 if (!anInputBands.empty())
     430             :                 {
     431          13 :                     anInputBands[i] = i + 1;
     432          13 :                     iVRTBand = i + 1;
     433             :                 }
     434          17 :                 poVDS->AddBand(poSrcBand->GetRasterDataType(), nullptr);
     435             :                 VRTSourcedRasterBand *poVRTBand =
     436           0 :                     dynamic_cast<VRTSourcedRasterBand *>(
     437          17 :                         poVDS->GetRasterBand(iVRTBand));
     438          17 :                 if (poVRTBand == nullptr)
     439           0 :                     return CE_Failure;
     440          17 :                 aMSBands[i] = poVRTBand;
     441          17 :                 poVRTBand->SetNoDataValue(psOptions->dfNoData);
     442             :                 const char *pszNBITS =
     443          17 :                     poSrcBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE");
     444          17 :                 if (pszNBITS)
     445           0 :                     poVRTBand->SetMetadataItem("NBITS", pszNBITS,
     446           0 :                                                "IMAGE_STRUCTURE");
     447             : 
     448          17 :                 VRTSimpleSource *poSimpleSource = new VRTSimpleSource();
     449          68 :                 poVRTBand->ConfigureSource(
     450             :                     poSimpleSource, poSrcBand, FALSE, 0, 0,
     451          17 :                     poSrcBand->GetXSize(), poSrcBand->GetYSize(), 0, 0,
     452          17 :                     poSrcBand->GetXSize(), poSrcBand->GetYSize());
     453          17 :                 poVRTBand->AddSource(poSimpleSource);
     454             :             }
     455             :         }
     456             :     }
     457             : 
     458             :     // Setup thread pool.
     459          69 :     int nThreads = psOptions->nThreads;
     460          69 :     if (nThreads == -1)
     461          23 :         nThreads = CPLGetNumCPUs();
     462          46 :     else if (nThreads == 0)
     463             :     {
     464             :         const char *pszNumThreads =
     465          46 :             CPLGetConfigOption("GDAL_NUM_THREADS", nullptr);
     466          46 :         if (pszNumThreads)
     467             :         {
     468           0 :             if (EQUAL(pszNumThreads, "ALL_CPUS"))
     469           0 :                 nThreads = CPLGetNumCPUs();
     470             :             else
     471           0 :                 nThreads = std::max(0, std::min(128, atoi(pszNumThreads)));
     472             :         }
     473             :     }
     474          69 :     if (nThreads > 1)
     475             :     {
     476          23 :         CPLDebug("PANSHARPEN", "Using %d threads", nThreads);
     477          46 :         poThreadPool = new (std::nothrow) CPLWorkerThreadPool();
     478             :         // coverity[tainted_data]
     479          46 :         if (poThreadPool == nullptr ||
     480          23 :             !poThreadPool->Setup(nThreads, nullptr, nullptr))
     481             :         {
     482           0 :             delete poThreadPool;
     483           0 :             poThreadPool = nullptr;
     484             :         }
     485             :     }
     486             : 
     487          69 :     GDALRIOResampleAlg eResampleAlg = psOptions->eResampleAlg;
     488          69 :     if (eResampleAlg != GRIORA_NearestNeighbour)
     489             :     {
     490          69 :         const char *pszResampling =
     491         138 :             (eResampleAlg == GRIORA_Bilinear)      ? "BILINEAR"
     492          69 :             : (eResampleAlg == GRIORA_Cubic)       ? "CUBIC"
     493           0 :             : (eResampleAlg == GRIORA_CubicSpline) ? "CUBICSPLINE"
     494           0 :             : (eResampleAlg == GRIORA_Lanczos)     ? "LANCZOS"
     495           0 :             : (eResampleAlg == GRIORA_Average)     ? "AVERAGE"
     496           0 :             : (eResampleAlg == GRIORA_RMS)         ? "RMS"
     497           0 :             : (eResampleAlg == GRIORA_Mode)        ? "MODE"
     498           0 :             : (eResampleAlg == GRIORA_Gauss)       ? "GAUSS"
     499             :                                                    : "UNKNOWN";
     500             : 
     501          69 :         GDALGetResampleFunction(pszResampling, &nKernelRadius);
     502             :     }
     503             : 
     504          69 :     return CE_None;
     505             : }
     506             : 
     507             : /************************************************************************/
     508             : /*                    WeightedBroveyWithNoData()                        */
     509             : /************************************************************************/
     510             : 
     511             : template <class WorkDataType, class OutDataType>
     512           8 : void GDALPansharpenOperation::WeightedBroveyWithNoData(
     513             :     const WorkDataType *pPanBuffer,
     514             :     const WorkDataType *pUpsampledSpectralBuffer, OutDataType *pDataBuf,
     515             :     size_t nValues, size_t nBandValues, WorkDataType nMaxValue) const
     516             : {
     517             :     WorkDataType noData, validValue;
     518           8 :     GDALCopyWord(psOptions->dfNoData, noData);
     519             : 
     520             :     if constexpr (!(std::numeric_limits<WorkDataType>::is_integer))
     521           0 :         validValue = static_cast<WorkDataType>(noData + 1e-5);
     522           8 :     else if (noData == std::numeric_limits<WorkDataType>::min())
     523           4 :         validValue = std::numeric_limits<WorkDataType>::min() + 1;
     524             :     else
     525           4 :         validValue = noData - 1;
     526             : 
     527     1578140 :     for (size_t j = 0; j < nValues; j++)
     528             :     {
     529     1578140 :         double dfPseudoPanchro = 0.0;
     530     5691390 :         for (int i = 0; i < psOptions->nInputSpectralBands; i++)
     531             :         {
     532     4113250 :             WorkDataType nSpectralVal =
     533     4113250 :                 pUpsampledSpectralBuffer[i * nBandValues + j];
     534     4113250 :             if (nSpectralVal == noData)
     535             :             {
     536           0 :                 dfPseudoPanchro = 0.0;
     537           0 :                 break;
     538             :             }
     539     4113250 :             dfPseudoPanchro += psOptions->padfWeights[i] * nSpectralVal;
     540             :         }
     541     1578140 :         if (dfPseudoPanchro != 0.0 && pPanBuffer[j] != noData)
     542             :         {
     543     1591670 :             const double dfFactor = pPanBuffer[j] / dfPseudoPanchro;
     544     5648550 :             for (int i = 0; i < psOptions->nOutPansharpenedBands; i++)
     545             :             {
     546     4065740 :                 WorkDataType nRawValue = pUpsampledSpectralBuffer
     547     4065740 :                     [psOptions->panOutPansharpenedBands[i] * nBandValues + j];
     548             :                 WorkDataType nPansharpenedValue;
     549     4065740 :                 GDALCopyWord(nRawValue * dfFactor, nPansharpenedValue);
     550     3989350 :                 if (nMaxValue != 0 && nPansharpenedValue > nMaxValue)
     551           0 :                     nPansharpenedValue = nMaxValue;
     552             :                 // We don't want a valid value to be mapped to NoData.
     553     3989350 :                 if (nPansharpenedValue == noData)
     554        5170 :                     nPansharpenedValue = validValue;
     555     3989350 :                 GDALCopyWord(nPansharpenedValue, pDataBuf[i * nBandValues + j]);
     556     1582820 :             }
     557             :         }
     558             :         else
     559             :         {
     560        1724 :             for (int i = 0; i < psOptions->nOutPansharpenedBands; i++)
     561             :             {
     562        6404 :                 GDALCopyWord(noData, pDataBuf[i * nBandValues + j]);
     563             :             }
     564             :         }
     565             :     }
     566           8 : }
     567             : 
     568             : /************************************************************************/
     569             : /*                           ComputeFactor()                            */
     570             : /************************************************************************/
     571             : 
     572             : template <class T>
     573    23267180 : static inline double ComputeFactor(T panValue, double dfPseudoPanchro)
     574             : {
     575    23267180 :     if (dfPseudoPanchro == 0.0)
     576        4163 :         return 0.0;
     577             : 
     578    23263020 :     return panValue / dfPseudoPanchro;
     579             : }
     580             : 
     581             : /************************************************************************/
     582             : /*                           ClampAndRound()                            */
     583             : /************************************************************************/
     584             : 
     585     1761948 : template <class T> static inline T ClampAndRound(double dfVal, T nMaxValue)
     586             : {
     587     1761948 :     if (dfVal > nMaxValue)
     588          50 :         return nMaxValue;
     589             :     else
     590     1761894 :         return static_cast<T>(dfVal + 0.5);
     591             : }
     592             : 
     593             : /************************************************************************/
     594             : /*                         WeightedBrovey()                             */
     595             : /************************************************************************/
     596             : 
     597             : template <class WorkDataType, class OutDataType, int bHasBitDepth>
     598         158 : void GDALPansharpenOperation::WeightedBrovey3(
     599             :     const WorkDataType *pPanBuffer,
     600             :     const WorkDataType *pUpsampledSpectralBuffer, OutDataType *pDataBuf,
     601             :     size_t nValues, size_t nBandValues, WorkDataType nMaxValue) const
     602             : {
     603         158 :     if (psOptions->bHasNoData)
     604             :     {
     605           8 :         WeightedBroveyWithNoData<WorkDataType, OutDataType>(
     606             :             pPanBuffer, pUpsampledSpectralBuffer, pDataBuf, nValues,
     607             :             nBandValues, nMaxValue);
     608           8 :         return;
     609             :     }
     610             : 
     611    22423312 :     for (size_t j = 0; j < nValues; j++)
     612             :     {
     613    22137190 :         double dfFactor = 0.0;
     614             :         // if( pPanBuffer[j] == 0 )
     615             :         //     dfFactor = 1.0;
     616             :         // else
     617             :         {
     618    22137190 :             double dfPseudoPanchro = 0.0;
     619    95761400 :             for (int i = 0; i < psOptions->nInputSpectralBands; i++)
     620    73624200 :                 dfPseudoPanchro +=
     621    73624200 :                     psOptions->padfWeights[i] *
     622    73624200 :                     pUpsampledSpectralBuffer[i * nBandValues + j];
     623    22137190 :             dfFactor = ComputeFactor(pPanBuffer[j], dfPseudoPanchro);
     624             :         }
     625             : 
     626    90722300 :         for (int i = 0; i < psOptions->nOutPansharpenedBands; i++)
     627             :         {
     628    68299000 :             WorkDataType nRawValue =
     629    68299000 :                 pUpsampledSpectralBuffer[psOptions->panOutPansharpenedBands[i] *
     630    68299000 :                                              nBandValues +
     631             :                                          j];
     632             :             WorkDataType nPansharpenedValue;
     633    68299000 :             GDALCopyWord(nRawValue * dfFactor, nPansharpenedValue);
     634             :             if constexpr (bHasBitDepth)
     635             :             {
     636           0 :                 if (nPansharpenedValue > nMaxValue)
     637           0 :                     nPansharpenedValue = nMaxValue;
     638             :             }
     639    66803800 :             GDALCopyWord(nPansharpenedValue, pDataBuf[i * nBandValues + j]);
     640             :         }
     641             :     }
     642             : }
     643             : 
     644             : /* We restrict to 64bit processors because they are guaranteed to have SSE2 */
     645             : /* Could possibly be used too on 32bit, but we would need to check at runtime */
     646             : #if defined(__x86_64) || defined(_M_X64) || defined(USE_NEON_OPTIMIZATIONS)
     647             : 
     648             : #define USE_SSE2
     649             : #include "gdalsse_priv.h"
     650             : 
     651             : template <class T, int NINPUT, int NOUTPUT>
     652          41 : size_t GDALPansharpenOperation::WeightedBroveyPositiveWeightsInternal(
     653             :     const T *pPanBuffer, const T *pUpsampledSpectralBuffer, T *pDataBuf,
     654             :     size_t nValues, size_t nBandValues, T nMaxValue) const
     655             : {
     656             :     static_assert(NINPUT == 3 || NINPUT == 4);
     657             :     static_assert(NOUTPUT == 3 || NOUTPUT == 4);
     658          41 :     const XMMReg4Double w0 =
     659          41 :         XMMReg4Double::Load1ValHighAndLow(psOptions->padfWeights + 0);
     660          39 :     const XMMReg4Double w1 =
     661          39 :         XMMReg4Double::Load1ValHighAndLow(psOptions->padfWeights + 1);
     662          39 :     const XMMReg4Double w2 =
     663          39 :         XMMReg4Double::Load1ValHighAndLow(psOptions->padfWeights + 2);
     664          40 :     [[maybe_unused]] const XMMReg4Double w3 =
     665             :         (NINPUT == 3)
     666             :             ? XMMReg4Double::Zero()
     667           8 :             : XMMReg4Double::Load1ValHighAndLow(psOptions->padfWeights + 3);
     668             : 
     669          40 :     const XMMReg4Double zero = XMMReg4Double::Zero();
     670          40 :     double dfMaxValue = nMaxValue;
     671          40 :     const XMMReg4Double maxValue =
     672             :         XMMReg4Double::Load1ValHighAndLow(&dfMaxValue);
     673             : 
     674        2890 :     size_t j = 0;  // Used after for.
     675     1337455 :     for (; j + 3 < nValues; j += 4)
     676             :     {
     677     1337414 :         XMMReg4Double pseudoPanchro = zero;
     678             : 
     679     1334867 :         XMMReg4Double val0 = XMMReg4Double::Load4Val(pUpsampledSpectralBuffer +
     680      777189 :                                                      0 * nBandValues + j);
     681     1338861 :         XMMReg4Double val1 = XMMReg4Double::Load4Val(pUpsampledSpectralBuffer +
     682     1338861 :                                                      1 * nBandValues + j);
     683     1332811 :         XMMReg4Double val2 = XMMReg4Double::Load4Val(pUpsampledSpectralBuffer +
     684     1332811 :                                                      2 * nBandValues + j);
     685     1337391 :         XMMReg4Double val3;
     686             :         if constexpr (NINPUT == 4 || NOUTPUT == 4)
     687             :         {
     688      518418 :             val3 = XMMReg4Double::Load4Val(pUpsampledSpectralBuffer +
     689      520440 :                                            3 * nBandValues + j);
     690             :         }
     691             : 
     692     1336471 :         pseudoPanchro += w0 * val0;
     693     1332845 :         pseudoPanchro += w1 * val1;
     694     1336008 :         pseudoPanchro += w2 * val2;
     695             :         if constexpr (NINPUT == 4)
     696      519215 :             pseudoPanchro += w3 * val3;
     697             : 
     698             :         /* Little trick to avoid use of ternary operator due to one of the
     699             :          * branch being zero */
     700     1337898 :         XMMReg4Double factor = XMMReg4Double::And(
     701             :             XMMReg4Double::NotEquals(pseudoPanchro, zero),
     702      779922 :             XMMReg4Double::Load4Val(pPanBuffer + j) / pseudoPanchro);
     703             : 
     704     1331998 :         val0 = XMMReg4Double::Min(val0 * factor, maxValue);
     705     1333733 :         val1 = XMMReg4Double::Min(val1 * factor, maxValue);
     706     1334434 :         val2 = XMMReg4Double::Min(val2 * factor, maxValue);
     707             :         if constexpr (NOUTPUT == 4)
     708             :         {
     709      260064 :             val3 = XMMReg4Double::Min(val3 * factor, maxValue);
     710             :         }
     711     1333226 :         val0.Store4Val(pDataBuf + 0 * nBandValues + j);
     712     1337221 :         val1.Store4Val(pDataBuf + 1 * nBandValues + j);
     713     1338069 :         val2.Store4Val(pDataBuf + 2 * nBandValues + j);
     714             :         if constexpr (NOUTPUT == 4)
     715             :         {
     716      260665 :             val3.Store4Val(pDataBuf + 3 * nBandValues + j);
     717             :         }
     718             :     }
     719          41 :     return j;
     720             : }
     721             : 
     722             : #else
     723             : 
     724             : template <class T, int NINPUT, int NOUTPUT>
     725             : size_t GDALPansharpenOperation::WeightedBroveyPositiveWeightsInternal(
     726             :     const T *pPanBuffer, const T *pUpsampledSpectralBuffer, T *pDataBuf,
     727             :     size_t nValues, size_t nBandValues, T nMaxValue) const
     728             : {
     729             :     static_assert(NINPUT == 3 || NINPUT == 4);
     730             :     const double dfw0 = psOptions->padfWeights[0];
     731             :     const double dfw1 = psOptions->padfWeights[1];
     732             :     const double dfw2 = psOptions->padfWeights[2];
     733             :     // cppcheck-suppress knownConditionTrueFalse
     734             :     [[maybe_unused]] const double dfw3 =
     735             :         (NINPUT == 3) ? 0 : psOptions->padfWeights[3];
     736             :     size_t j = 0;  // Used after for.
     737             :     for (; j + 1 < nValues; j += 2)
     738             :     {
     739             :         double dfFactor = 0.0;
     740             :         double dfFactor2 = 0.0;
     741             :         double dfPseudoPanchro = 0.0;
     742             :         double dfPseudoPanchro2 = 0.0;
     743             : 
     744             :         dfPseudoPanchro += dfw0 * pUpsampledSpectralBuffer[j];
     745             :         dfPseudoPanchro2 += dfw0 * pUpsampledSpectralBuffer[j + 1];
     746             : 
     747             :         dfPseudoPanchro += dfw1 * pUpsampledSpectralBuffer[nBandValues + j];
     748             :         dfPseudoPanchro2 +=
     749             :             dfw1 * pUpsampledSpectralBuffer[nBandValues + j + 1];
     750             : 
     751             :         dfPseudoPanchro += dfw2 * pUpsampledSpectralBuffer[2 * nBandValues + j];
     752             :         dfPseudoPanchro2 +=
     753             :             dfw2 * pUpsampledSpectralBuffer[2 * nBandValues + j + 1];
     754             : 
     755             :         if constexpr (NINPUT == 4)
     756             :         {
     757             :             dfPseudoPanchro +=
     758             :                 dfw3 * pUpsampledSpectralBuffer[3 * nBandValues + j];
     759             :             dfPseudoPanchro2 +=
     760             :                 dfw3 * pUpsampledSpectralBuffer[3 * nBandValues + j + 1];
     761             :         }
     762             : 
     763             :         dfFactor = ComputeFactor(pPanBuffer[j], dfPseudoPanchro);
     764             :         dfFactor2 = ComputeFactor(pPanBuffer[j + 1], dfPseudoPanchro2);
     765             : 
     766             :         for (int i = 0; i < NOUTPUT; i++)
     767             :         {
     768             :             T nRawValue = pUpsampledSpectralBuffer[i * nBandValues + j];
     769             :             double dfTmp = nRawValue * dfFactor;
     770             :             pDataBuf[i * nBandValues + j] = ClampAndRound(dfTmp, nMaxValue);
     771             : 
     772             :             T nRawValue2 = pUpsampledSpectralBuffer[i * nBandValues + j + 1];
     773             :             double dfTmp2 = nRawValue2 * dfFactor2;
     774             :             pDataBuf[i * nBandValues + j + 1] =
     775             :                 ClampAndRound(dfTmp2, nMaxValue);
     776             :         }
     777             :     }
     778             :     return j;
     779             : }
     780             : #endif
     781             : 
     782             : template <class T>
     783          63 : void GDALPansharpenOperation::WeightedBroveyPositiveWeights(
     784             :     const T *pPanBuffer, const T *pUpsampledSpectralBuffer, T *pDataBuf,
     785             :     size_t nValues, size_t nBandValues, T nMaxValue) const
     786             : {
     787          63 :     if (psOptions->bHasNoData)
     788             :     {
     789           0 :         WeightedBroveyWithNoData<T, T>(pPanBuffer, pUpsampledSpectralBuffer,
     790             :                                        pDataBuf, nValues, nBandValues,
     791             :                                        nMaxValue);
     792           0 :         return;
     793             :     }
     794             : 
     795          63 :     if (nMaxValue == 0)
     796          56 :         nMaxValue = cpl::NumericLimits<T>::max();
     797             :     size_t j;
     798          62 :     if (psOptions->nInputSpectralBands == 3 &&
     799          32 :         psOptions->nOutPansharpenedBands == 3 &&
     800          32 :         psOptions->panOutPansharpenedBands[0] == 0 &&
     801          33 :         psOptions->panOutPansharpenedBands[1] == 1 &&
     802          32 :         psOptions->panOutPansharpenedBands[2] == 2)
     803             :     {
     804          32 :         j = WeightedBroveyPositiveWeightsInternal<T, 3, 3>(
     805             :             pPanBuffer, pUpsampledSpectralBuffer, pDataBuf, nValues,
     806             :             nBandValues, nMaxValue);
     807             :     }
     808          30 :     else if (psOptions->nInputSpectralBands == 4 &&
     809           8 :              psOptions->nOutPansharpenedBands == 4 &&
     810           4 :              psOptions->panOutPansharpenedBands[0] == 0 &&
     811           4 :              psOptions->panOutPansharpenedBands[1] == 1 &&
     812           4 :              psOptions->panOutPansharpenedBands[2] == 2 &&
     813           4 :              psOptions->panOutPansharpenedBands[3] == 3)
     814             :     {
     815           4 :         j = WeightedBroveyPositiveWeightsInternal<T, 4, 4>(
     816             :             pPanBuffer, pUpsampledSpectralBuffer, pDataBuf, nValues,
     817             :             nBandValues, nMaxValue);
     818             :     }
     819          26 :     else if (psOptions->nInputSpectralBands == 4 &&
     820           4 :              psOptions->nOutPansharpenedBands == 3 &&
     821           4 :              psOptions->panOutPansharpenedBands[0] == 0 &&
     822           4 :              psOptions->panOutPansharpenedBands[1] == 1 &&
     823           4 :              psOptions->panOutPansharpenedBands[2] == 2)
     824             :     {
     825           4 :         j = WeightedBroveyPositiveWeightsInternal<T, 4, 3>(
     826             :             pPanBuffer, pUpsampledSpectralBuffer, pDataBuf, nValues,
     827             :             nBandValues, nMaxValue);
     828             :     }
     829             :     else
     830             :     {
     831      433569 :         for (j = 0; j + 1 < nValues; j += 2)
     832             :         {
     833      475408 :             double dfFactor = 0.0;
     834      475408 :             double dfFactor2 = 0.0;
     835      475408 :             double dfPseudoPanchro = 0.0;
     836      475408 :             double dfPseudoPanchro2 = 0.0;
     837     1409228 :             for (int i = 0; i < psOptions->nInputSpectralBands; i++)
     838             :             {
     839      933825 :                 dfPseudoPanchro +=
     840      933825 :                     psOptions->padfWeights[i] *
     841      933825 :                     pUpsampledSpectralBuffer[i * nBandValues + j];
     842      933825 :                 dfPseudoPanchro2 +=
     843      933825 :                     psOptions->padfWeights[i] *
     844      933825 :                     pUpsampledSpectralBuffer[i * nBandValues + j + 1];
     845             :             }
     846             : 
     847      475408 :             dfFactor = ComputeFactor(pPanBuffer[j], dfPseudoPanchro);
     848      476646 :             dfFactor2 = ComputeFactor(pPanBuffer[j + 1], dfPseudoPanchro2);
     849             : 
     850     1256318 :             for (int i = 0; i < psOptions->nOutPansharpenedBands; i++)
     851             :             {
     852      822768 :                 const T nRawValue = pUpsampledSpectralBuffer
     853      822768 :                     [psOptions->panOutPansharpenedBands[i] * nBandValues + j];
     854      822768 :                 const double dfTmp = nRawValue * dfFactor;
     855      822768 :                 pDataBuf[i * nBandValues + j] = ClampAndRound(dfTmp, nMaxValue);
     856             : 
     857      767336 :                 const T nRawValue2 = pUpsampledSpectralBuffer
     858      767336 :                     [psOptions->panOutPansharpenedBands[i] * nBandValues + j +
     859             :                      1];
     860      767336 :                 const double dfTmp2 = nRawValue2 * dfFactor2;
     861      781403 :                 pDataBuf[i * nBandValues + j + 1] =
     862      767336 :                     ClampAndRound(dfTmp2, nMaxValue);
     863             :             }
     864             :         }
     865             :     }
     866          18 :     for (; j < nValues; j++)
     867             :     {
     868          13 :         double dfFactor = 0.0;
     869          13 :         double dfPseudoPanchro = 0.0;
     870          54 :         for (int i = 0; i < psOptions->nInputSpectralBands; i++)
     871          41 :             dfPseudoPanchro += psOptions->padfWeights[i] *
     872          41 :                                pUpsampledSpectralBuffer[i * nBandValues + j];
     873          13 :         dfFactor = ComputeFactor(pPanBuffer[j], dfPseudoPanchro);
     874             : 
     875          53 :         for (int i = 0; i < psOptions->nOutPansharpenedBands; i++)
     876             :         {
     877          40 :             T nRawValue =
     878          40 :                 pUpsampledSpectralBuffer[psOptions->panOutPansharpenedBands[i] *
     879          40 :                                              nBandValues +
     880             :                                          j];
     881          40 :             double dfTmp = nRawValue * dfFactor;
     882          40 :             pDataBuf[i * nBandValues + j] = ClampAndRound(dfTmp, nMaxValue);
     883             :         }
     884             :     }
     885             : }
     886             : 
     887             : template <class WorkDataType, class OutDataType>
     888         154 : void GDALPansharpenOperation::WeightedBrovey(
     889             :     const WorkDataType *pPanBuffer,
     890             :     const WorkDataType *pUpsampledSpectralBuffer, OutDataType *pDataBuf,
     891             :     size_t nValues, size_t nBandValues, WorkDataType nMaxValue) const
     892             : {
     893         154 :     if (nMaxValue == 0)
     894         154 :         WeightedBrovey3<WorkDataType, OutDataType, FALSE>(
     895             :             pPanBuffer, pUpsampledSpectralBuffer, pDataBuf, nValues,
     896             :             nBandValues, 0);
     897             :     else
     898             :     {
     899           0 :         WeightedBrovey3<WorkDataType, OutDataType, TRUE>(
     900             :             pPanBuffer, pUpsampledSpectralBuffer, pDataBuf, nValues,
     901             :             nBandValues, nMaxValue);
     902             :     }
     903         154 : }
     904             : 
     905             : template <class T>
     906          63 : void GDALPansharpenOperation::WeightedBroveyGByteOrUInt16(
     907             :     const T *pPanBuffer, const T *pUpsampledSpectralBuffer, T *pDataBuf,
     908             :     size_t nValues, size_t nBandValues, T nMaxValue) const
     909             : {
     910          63 :     if (bPositiveWeights)
     911             :     {
     912          63 :         WeightedBroveyPositiveWeights(pPanBuffer, pUpsampledSpectralBuffer,
     913             :                                       pDataBuf, nValues, nBandValues,
     914             :                                       nMaxValue);
     915             :     }
     916           0 :     else if (nMaxValue == 0)
     917             :     {
     918           0 :         WeightedBrovey3<T, T, FALSE>(pPanBuffer, pUpsampledSpectralBuffer,
     919             :                                      pDataBuf, nValues, nBandValues, 0);
     920             :     }
     921             :     else
     922             :     {
     923           0 :         WeightedBrovey3<T, T, TRUE>(pPanBuffer, pUpsampledSpectralBuffer,
     924             :                                     pDataBuf, nValues, nBandValues, nMaxValue);
     925             :     }
     926          63 : }
     927             : 
     928             : template <>
     929          48 : void GDALPansharpenOperation::WeightedBrovey<GByte, GByte>(
     930             :     const GByte *pPanBuffer, const GByte *pUpsampledSpectralBuffer,
     931             :     GByte *pDataBuf, size_t nValues, size_t nBandValues, GByte nMaxValue) const
     932             : {
     933          48 :     WeightedBroveyGByteOrUInt16(pPanBuffer, pUpsampledSpectralBuffer, pDataBuf,
     934             :                                 nValues, nBandValues, nMaxValue);
     935          48 : }
     936             : 
     937             : template <>
     938          15 : void GDALPansharpenOperation::WeightedBrovey<GUInt16, GUInt16>(
     939             :     const GUInt16 *pPanBuffer, const GUInt16 *pUpsampledSpectralBuffer,
     940             :     GUInt16 *pDataBuf, size_t nValues, size_t nBandValues,
     941             :     GUInt16 nMaxValue) const
     942             : {
     943          15 :     WeightedBroveyGByteOrUInt16(pPanBuffer, pUpsampledSpectralBuffer, pDataBuf,
     944             :                                 nValues, nBandValues, nMaxValue);
     945          15 : }
     946             : 
     947             : template <class WorkDataType>
     948         216 : CPLErr GDALPansharpenOperation::WeightedBrovey(
     949             :     const WorkDataType *pPanBuffer,
     950             :     const WorkDataType *pUpsampledSpectralBuffer, void *pDataBuf,
     951             :     GDALDataType eBufDataType, size_t nValues, size_t nBandValues,
     952             :     WorkDataType nMaxValue) const
     953             : {
     954         216 :     switch (eBufDataType)
     955             :     {
     956          49 :         case GDT_Byte:
     957          49 :             WeightedBrovey(pPanBuffer, pUpsampledSpectralBuffer,
     958             :                            static_cast<GByte *>(pDataBuf), nValues, nBandValues,
     959             :                            nMaxValue);
     960          49 :             break;
     961             : 
     962          16 :         case GDT_UInt16:
     963          16 :             WeightedBrovey(pPanBuffer, pUpsampledSpectralBuffer,
     964             :                            static_cast<GUInt16 *>(pDataBuf), nValues,
     965             :                            nBandValues, nMaxValue);
     966          16 :             break;
     967             : 
     968             : #ifndef LIMIT_TYPES
     969             :         case GDT_Int8:
     970             :             WeightedBrovey(pPanBuffer, pUpsampledSpectralBuffer,
     971             :                            static_cast<GInt8 *>(pDataBuf), nValues, nBandValues,
     972             :                            nMaxValue);
     973             :             break;
     974             : 
     975             :         case GDT_Int16:
     976             :             WeightedBrovey(pPanBuffer, pUpsampledSpectralBuffer,
     977             :                            static_cast<GInt16 *>(pDataBuf), nValues,
     978             :                            nBandValues, nMaxValue);
     979             :             break;
     980             : 
     981             :         case GDT_UInt32:
     982             :             WeightedBrovey(pPanBuffer, pUpsampledSpectralBuffer,
     983             :                            static_cast<GUInt32 *>(pDataBuf), nValues,
     984             :                            nBandValues, nMaxValue);
     985             :             break;
     986             : 
     987             :         case GDT_Int32:
     988             :             WeightedBrovey(pPanBuffer, pUpsampledSpectralBuffer,
     989             :                            static_cast<GInt32 *>(pDataBuf), nValues,
     990             :                            nBandValues, nMaxValue);
     991             :             break;
     992             : 
     993             :         case GDT_UInt64:
     994             :             WeightedBrovey(pPanBuffer, pUpsampledSpectralBuffer,
     995             :                            static_cast<std::uint64_t *>(pDataBuf), nValues,
     996             :                            nBandValues, nMaxValue);
     997             :             break;
     998             : 
     999             :         case GDT_Int64:
    1000             :             WeightedBrovey(pPanBuffer, pUpsampledSpectralBuffer,
    1001             :                            static_cast<std::int64_t *>(pDataBuf), nValues,
    1002             :                            nBandValues, nMaxValue);
    1003             :             break;
    1004             : 
    1005             :         case GDT_Float16:
    1006             :             WeightedBrovey(pPanBuffer, pUpsampledSpectralBuffer,
    1007             :                            static_cast<GFloat16 *>(pDataBuf), nValues,
    1008             :                            nBandValues, nMaxValue);
    1009             :             break;
    1010             : 
    1011             :         case GDT_Float32:
    1012             :             WeightedBrovey(pPanBuffer, pUpsampledSpectralBuffer,
    1013             :                            static_cast<float *>(pDataBuf), nValues, nBandValues,
    1014             :                            nMaxValue);
    1015             :             break;
    1016             : #endif
    1017             : 
    1018         151 :         case GDT_Float64:
    1019         151 :             WeightedBrovey(pPanBuffer, pUpsampledSpectralBuffer,
    1020             :                            static_cast<double *>(pDataBuf), nValues,
    1021             :                            nBandValues, nMaxValue);
    1022         152 :             break;
    1023             : 
    1024           0 :         default:
    1025           0 :             CPLError(CE_Failure, CPLE_NotSupported,
    1026             :                      "eBufDataType not supported");
    1027           0 :             return CE_Failure;
    1028             :             break;
    1029             :     }
    1030             : 
    1031         217 :     return CE_None;
    1032             : }
    1033             : 
    1034             : template <class WorkDataType>
    1035           6 : CPLErr GDALPansharpenOperation::WeightedBrovey(
    1036             :     const WorkDataType *pPanBuffer,
    1037             :     const WorkDataType *pUpsampledSpectralBuffer, void *pDataBuf,
    1038             :     GDALDataType eBufDataType, size_t nValues, size_t nBandValues) const
    1039             : {
    1040           6 :     switch (eBufDataType)
    1041             :     {
    1042           6 :         case GDT_Byte:
    1043           6 :             WeightedBrovey3<WorkDataType, GByte, FALSE>(
    1044             :                 pPanBuffer, pUpsampledSpectralBuffer,
    1045             :                 static_cast<GByte *>(pDataBuf), nValues, nBandValues, 0);
    1046           6 :             break;
    1047             : 
    1048           0 :         case GDT_UInt16:
    1049           0 :             WeightedBrovey3<WorkDataType, GUInt16, FALSE>(
    1050             :                 pPanBuffer, pUpsampledSpectralBuffer,
    1051             :                 static_cast<GUInt16 *>(pDataBuf), nValues, nBandValues, 0);
    1052           0 :             break;
    1053             : 
    1054             : #ifndef LIMIT_TYPES
    1055             :         case GDT_Int8:
    1056             :             WeightedBrovey3<WorkDataType, GInt8, FALSE>(
    1057             :                 pPanBuffer, pUpsampledSpectralBuffer,
    1058             :                 static_cast<GInt8 *>(pDataBuf), nValues, nBandValues, 0);
    1059             :             break;
    1060             : 
    1061             :         case GDT_Int16:
    1062             :             WeightedBrovey3<WorkDataType, GInt16, FALSE>(
    1063             :                 pPanBuffer, pUpsampledSpectralBuffer,
    1064             :                 static_cast<GInt16 *>(pDataBuf), nValues, nBandValues, 0);
    1065             :             break;
    1066             : 
    1067             :         case GDT_UInt32:
    1068             :             WeightedBrovey3<WorkDataType, GUInt32, FALSE>(
    1069             :                 pPanBuffer, pUpsampledSpectralBuffer,
    1070             :                 static_cast<GUInt32 *>(pDataBuf), nValues, nBandValues, 0);
    1071             :             break;
    1072             : 
    1073             :         case GDT_Int32:
    1074             :             WeightedBrovey3<WorkDataType, GInt32, FALSE>(
    1075             :                 pPanBuffer, pUpsampledSpectralBuffer,
    1076             :                 static_cast<GInt32 *>(pDataBuf), nValues, nBandValues, 0);
    1077             :             break;
    1078             : 
    1079             :         case GDT_UInt64:
    1080             :             WeightedBrovey3<WorkDataType, std::uint64_t, FALSE>(
    1081             :                 pPanBuffer, pUpsampledSpectralBuffer,
    1082             :                 static_cast<std::uint64_t *>(pDataBuf), nValues, nBandValues,
    1083             :                 0);
    1084             :             break;
    1085             : 
    1086             :         case GDT_Int64:
    1087             :             WeightedBrovey3<WorkDataType, std::int64_t, FALSE>(
    1088             :                 pPanBuffer, pUpsampledSpectralBuffer,
    1089             :                 static_cast<std::int64_t *>(pDataBuf), nValues, nBandValues, 0);
    1090             :             break;
    1091             : 
    1092             :         case GDT_Float16:
    1093             :             WeightedBrovey3<WorkDataType, GFloat16, FALSE>(
    1094             :                 pPanBuffer, pUpsampledSpectralBuffer,
    1095             :                 static_cast<GFloat16 *>(pDataBuf), nValues, nBandValues, 0);
    1096             :             break;
    1097             : 
    1098             :         case GDT_Float32:
    1099             :             WeightedBrovey3<WorkDataType, float, FALSE>(
    1100             :                 pPanBuffer, pUpsampledSpectralBuffer,
    1101             :                 static_cast<float *>(pDataBuf), nValues, nBandValues, 0);
    1102             :             break;
    1103             : #endif
    1104             : 
    1105           0 :         case GDT_Float64:
    1106           0 :             WeightedBrovey3<WorkDataType, double, FALSE>(
    1107             :                 pPanBuffer, pUpsampledSpectralBuffer,
    1108             :                 static_cast<double *>(pDataBuf), nValues, nBandValues, 0);
    1109           0 :             break;
    1110             : 
    1111           0 :         default:
    1112           0 :             CPLError(CE_Failure, CPLE_NotSupported,
    1113             :                      "eBufDataType not supported");
    1114           0 :             return CE_Failure;
    1115             :             break;
    1116             :     }
    1117             : 
    1118           6 :     return CE_None;
    1119             : }
    1120             : 
    1121             : /************************************************************************/
    1122             : /*                           ClampValues()                              */
    1123             : /************************************************************************/
    1124             : 
    1125             : template <class T>
    1126           2 : static void ClampValues(T *panBuffer, size_t nValues, T nMaxVal)
    1127             : {
    1128          34 :     for (size_t i = 0; i < nValues; i++)
    1129             :     {
    1130          32 :         if (panBuffer[i] > nMaxVal)
    1131           8 :             panBuffer[i] = nMaxVal;
    1132             :     }
    1133           2 : }
    1134             : 
    1135             : /************************************************************************/
    1136             : /*                         ProcessRegion()                              */
    1137             : /************************************************************************/
    1138             : 
    1139             : /** Executes a pansharpening operation on a rectangular region of the
    1140             :  * resulting dataset.
    1141             :  *
    1142             :  * The window is expressed with respect to the dimensions of the panchromatic
    1143             :  * band.
    1144             :  *
    1145             :  * Spectral bands are upsampled and merged with the panchromatic band according
    1146             :  * to the select algorithm and options.
    1147             :  *
    1148             :  * @param nXOff pixel offset.
    1149             :  * @param nYOff pixel offset.
    1150             :  * @param nXSize width of the pansharpened region to compute.
    1151             :  * @param nYSize height of the pansharpened region to compute.
    1152             :  * @param pDataBuf output buffer. Must be nXSize * nYSize *
    1153             :  *                 GDALGetDataTypeSizeBytes(eBufDataType) *
    1154             :  *                 psOptions->nOutPansharpenedBands large.
    1155             :  *                 It begins with all values of the first output band, followed
    1156             :  *                 by values of the second output band, etc...
    1157             :  * @param eBufDataType data type of the output buffer
    1158             :  *
    1159             :  * @return CE_None in case of success, CE_Failure in case of failure.
    1160             :  *
    1161             :  */
    1162          91 : CPLErr GDALPansharpenOperation::ProcessRegion(int nXOff, int nYOff, int nXSize,
    1163             :                                               int nYSize, void *pDataBuf,
    1164             :                                               GDALDataType eBufDataType)
    1165             : {
    1166          91 :     if (psOptions == nullptr)
    1167           0 :         return CE_Failure;
    1168             : 
    1169             :     // TODO: Avoid allocating buffers each time.
    1170             :     GDALRasterBand *poPanchroBand =
    1171          91 :         GDALRasterBand::FromHandle(psOptions->hPanchroBand);
    1172          91 :     GDALDataType eWorkDataType = poPanchroBand->GetRasterDataType();
    1173             : #ifdef LIMIT_TYPES
    1174          91 :     if (eWorkDataType != GDT_Byte && eWorkDataType != GDT_UInt16)
    1175           6 :         eWorkDataType = GDT_Float64;
    1176             : #endif
    1177          91 :     const int nDataTypeSize = GDALGetDataTypeSizeBytes(eWorkDataType);
    1178          91 :     GByte *pUpsampledSpectralBuffer = static_cast<GByte *>(VSI_MALLOC3_VERBOSE(
    1179             :         nXSize, nYSize,
    1180             :         cpl::fits_on<int>(psOptions->nInputSpectralBands * nDataTypeSize)));
    1181             :     GByte *pPanBuffer = static_cast<GByte *>(
    1182          91 :         VSI_MALLOC3_VERBOSE(nXSize, nYSize, nDataTypeSize));
    1183          91 :     if (pUpsampledSpectralBuffer == nullptr || pPanBuffer == nullptr)
    1184             :     {
    1185           0 :         VSIFree(pUpsampledSpectralBuffer);
    1186           0 :         VSIFree(pPanBuffer);
    1187           0 :         return CE_Failure;
    1188             :     }
    1189             : 
    1190          91 :     CPLErr eErr = poPanchroBand->RasterIO(GF_Read, nXOff, nYOff, nXSize, nYSize,
    1191             :                                           pPanBuffer, nXSize, nYSize,
    1192             :                                           eWorkDataType, 0, 0, nullptr);
    1193          91 :     if (eErr != CE_None)
    1194             :     {
    1195           0 :         VSIFree(pUpsampledSpectralBuffer);
    1196           0 :         VSIFree(pPanBuffer);
    1197           0 :         return CE_Failure;
    1198             :     }
    1199             : 
    1200          91 :     int nTasks = 0;
    1201          91 :     if (poThreadPool)
    1202             :     {
    1203          44 :         nTasks = poThreadPool->GetThreadCount();
    1204          44 :         if (nTasks > nYSize)
    1205           0 :             nTasks = nYSize;
    1206             :     }
    1207             : 
    1208             :     GDALRasterIOExtraArg sExtraArg;
    1209          91 :     INIT_RASTERIO_EXTRA_ARG(sExtraArg);
    1210          91 :     const GDALRIOResampleAlg eResampleAlg = psOptions->eResampleAlg;
    1211             :     // cppcheck-suppress redundantAssignment
    1212          91 :     sExtraArg.eResampleAlg = eResampleAlg;
    1213          91 :     sExtraArg.bFloatingPointWindowValidity = TRUE;
    1214          91 :     sExtraArg.dfXOff = m_panToMSGT[0] + nXOff * m_panToMSGT[1];
    1215          91 :     sExtraArg.dfYOff = m_panToMSGT[3] + nYOff * m_panToMSGT[5];
    1216          91 :     sExtraArg.dfXSize = nXSize * m_panToMSGT[1];
    1217          91 :     sExtraArg.dfYSize = nYSize * m_panToMSGT[5];
    1218          91 :     if (sExtraArg.dfXOff + sExtraArg.dfXSize > aMSBands[0]->GetXSize())
    1219           2 :         sExtraArg.dfXSize = aMSBands[0]->GetXSize() - sExtraArg.dfXOff;
    1220          91 :     if (sExtraArg.dfYOff + sExtraArg.dfYSize > aMSBands[0]->GetYSize())
    1221           2 :         sExtraArg.dfYSize = aMSBands[0]->GetYSize() - sExtraArg.dfYOff;
    1222          91 :     int nSpectralXOff = static_cast<int>(sExtraArg.dfXOff);
    1223          91 :     int nSpectralYOff = static_cast<int>(sExtraArg.dfYOff);
    1224          91 :     int nSpectralXSize = static_cast<int>(0.49999 + sExtraArg.dfXSize);
    1225          91 :     int nSpectralYSize = static_cast<int>(0.49999 + sExtraArg.dfYSize);
    1226          91 :     if (nSpectralXOff + nSpectralXSize > aMSBands[0]->GetXSize())
    1227           0 :         nSpectralXSize = aMSBands[0]->GetXSize() - nSpectralXOff;
    1228          91 :     if (nSpectralYOff + nSpectralYSize > aMSBands[0]->GetYSize())
    1229           0 :         nSpectralYSize = aMSBands[0]->GetYSize() - nSpectralYOff;
    1230          91 :     if (nSpectralXSize == 0)
    1231          10 :         nSpectralXSize = 1;
    1232          91 :     if (nSpectralYSize == 0)
    1233          10 :         nSpectralYSize = 1;
    1234             : 
    1235             :     // When upsampling, extract the multispectral data at
    1236             :     // full resolution in a temp buffer, and then do the upsampling.
    1237          91 :     if (nSpectralXSize < nXSize && nSpectralYSize < nYSize &&
    1238          81 :         eResampleAlg != GRIORA_NearestNeighbour && nYSize > 1)
    1239             :     {
    1240             :         // Take some margin to take into account the radius of the
    1241             :         // resampling kernel.
    1242          81 :         int nXOffExtract = nSpectralXOff - nKernelRadius;
    1243          81 :         int nYOffExtract = nSpectralYOff - nKernelRadius;
    1244          81 :         int nXSizeExtract = nSpectralXSize + 1 + 2 * nKernelRadius;
    1245          81 :         int nYSizeExtract = nSpectralYSize + 1 + 2 * nKernelRadius;
    1246          81 :         if (nXOffExtract < 0)
    1247             :         {
    1248          74 :             nXSizeExtract += nXOffExtract;
    1249          74 :             nXOffExtract = 0;
    1250             :         }
    1251          81 :         if (nYOffExtract < 0)
    1252             :         {
    1253          67 :             nYSizeExtract += nYOffExtract;
    1254          67 :             nYOffExtract = 0;
    1255             :         }
    1256          81 :         if (nXOffExtract + nXSizeExtract > aMSBands[0]->GetXSize())
    1257          74 :             nXSizeExtract = aMSBands[0]->GetXSize() - nXOffExtract;
    1258          81 :         if (nYOffExtract + nYSizeExtract > aMSBands[0]->GetYSize())
    1259          67 :             nYSizeExtract = aMSBands[0]->GetYSize() - nYOffExtract;
    1260             : 
    1261          81 :         GByte *pSpectralBuffer = static_cast<GByte *>(VSI_MALLOC3_VERBOSE(
    1262             :             nXSizeExtract, nYSizeExtract,
    1263             :             cpl::fits_on<int>(psOptions->nInputSpectralBands * nDataTypeSize)));
    1264          81 :         if (pSpectralBuffer == nullptr)
    1265             :         {
    1266           0 :             VSIFree(pUpsampledSpectralBuffer);
    1267           0 :             VSIFree(pPanBuffer);
    1268           0 :             return CE_Failure;
    1269             :         }
    1270             : 
    1271          81 :         if (!anInputBands.empty())
    1272             :         {
    1273             :             // Use dataset RasterIO when possible.
    1274         146 :             eErr = aMSBands[0]->GetDataset()->RasterIO(
    1275             :                 GF_Read, nXOffExtract, nYOffExtract, nXSizeExtract,
    1276             :                 nYSizeExtract, pSpectralBuffer, nXSizeExtract, nYSizeExtract,
    1277          73 :                 eWorkDataType, static_cast<int>(anInputBands.size()),
    1278          73 :                 &anInputBands[0], 0, 0, 0, nullptr);
    1279             :         }
    1280             :         else
    1281             :         {
    1282           8 :             for (int i = 0;
    1283          20 :                  eErr == CE_None && i < psOptions->nInputSpectralBands; i++)
    1284             :             {
    1285          24 :                 eErr = aMSBands[i]->RasterIO(
    1286             :                     GF_Read, nXOffExtract, nYOffExtract, nXSizeExtract,
    1287             :                     nYSizeExtract,
    1288          12 :                     pSpectralBuffer + static_cast<size_t>(i) * nXSizeExtract *
    1289          12 :                                           nYSizeExtract * nDataTypeSize,
    1290             :                     nXSizeExtract, nYSizeExtract, eWorkDataType, 0, 0, nullptr);
    1291             :             }
    1292             :         }
    1293          81 :         if (eErr != CE_None)
    1294             :         {
    1295           0 :             VSIFree(pSpectralBuffer);
    1296           0 :             VSIFree(pUpsampledSpectralBuffer);
    1297           0 :             VSIFree(pPanBuffer);
    1298           0 :             return CE_Failure;
    1299             :         }
    1300             : 
    1301             :         // Create a MEM dataset that wraps the input buffer.
    1302          81 :         auto poMEMDS = MEMDataset::Create("", nXSizeExtract, nYSizeExtract, 0,
    1303             :                                           eWorkDataType, nullptr);
    1304             : 
    1305         322 :         for (int i = 0; i < psOptions->nInputSpectralBands; i++)
    1306             :         {
    1307         241 :             GByte *pabyBuffer =
    1308         241 :                 pSpectralBuffer + static_cast<size_t>(i) * nDataTypeSize *
    1309         241 :                                       nXSizeExtract * nYSizeExtract;
    1310         241 :             GDALRasterBandH hMEMBand = MEMCreateRasterBandEx(
    1311             :                 poMEMDS, i + 1, pabyBuffer, eWorkDataType, 0, 0, false);
    1312         241 :             poMEMDS->AddMEMBand(hMEMBand);
    1313             : 
    1314             :             const char *pszNBITS =
    1315         241 :                 aMSBands[i]->GetMetadataItem("NBITS", "IMAGE_STRUCTURE");
    1316         241 :             if (pszNBITS)
    1317           4 :                 poMEMDS->GetRasterBand(i + 1)->SetMetadataItem(
    1318           4 :                     "NBITS", pszNBITS, "IMAGE_STRUCTURE");
    1319             : 
    1320         241 :             if (psOptions->bHasNoData)
    1321          13 :                 poMEMDS->GetRasterBand(i + 1)->SetNoDataValue(
    1322          13 :                     psOptions->dfNoData);
    1323             :         }
    1324             : 
    1325          81 :         if (nTasks <= 1)
    1326             :         {
    1327          37 :             nSpectralXOff -= nXOffExtract;
    1328          37 :             nSpectralYOff -= nYOffExtract;
    1329          37 :             sExtraArg.dfXOff -= nXOffExtract;
    1330          37 :             sExtraArg.dfYOff -= nYOffExtract;
    1331          37 :             CPL_IGNORE_RET_VAL(poMEMDS->RasterIO(
    1332             :                 GF_Read, nSpectralXOff, nSpectralYOff, nSpectralXSize,
    1333             :                 nSpectralYSize, pUpsampledSpectralBuffer, nXSize, nYSize,
    1334          37 :                 eWorkDataType, psOptions->nInputSpectralBands, nullptr, 0, 0, 0,
    1335             :                 &sExtraArg));
    1336             :         }
    1337             :         else
    1338             :         {
    1339             :             // We are abusing the contract of the GDAL API by using the
    1340             :             // MEMDataset from several threads. In this case, this is safe. In
    1341             :             // case, that would no longer be the case we could create as many
    1342             :             // MEMDataset as threads pointing to the same buffer.
    1343             : 
    1344             :             // To avoid races in threads, we query now the mask flags,
    1345             :             // so that implicit mask bands are created now.
    1346         190 :             for (int i = 0; i < poMEMDS->GetRasterCount(); i++)
    1347             :             {
    1348         146 :                 poMEMDS->GetRasterBand(i + 1)->GetMaskFlags();
    1349             :             }
    1350             : 
    1351          44 :             std::vector<GDALPansharpenResampleJob> asJobs;
    1352          44 :             asJobs.resize(nTasks);
    1353          44 :             GDALPansharpenResampleJob *pasJobs = &(asJobs[0]);
    1354             :             {
    1355          44 :                 std::vector<void *> ahJobData;
    1356          44 :                 ahJobData.resize(nTasks);
    1357             : 
    1358             : #ifdef DEBUG_TIMING
    1359             :                 struct timeval tv;
    1360             : #endif
    1361         220 :                 for (int i = 0; i < nTasks; i++)
    1362             :                 {
    1363         176 :                     const size_t iStartLine =
    1364         176 :                         (static_cast<size_t>(i) * nYSize) / nTasks;
    1365         176 :                     const size_t iNextStartLine =
    1366         176 :                         (static_cast<size_t>(i + 1) * nYSize) / nTasks;
    1367         176 :                     pasJobs[i].poMEMDS = poMEMDS;
    1368         176 :                     pasJobs[i].eResampleAlg = eResampleAlg;
    1369         176 :                     pasJobs[i].dfXOff = sExtraArg.dfXOff - nXOffExtract;
    1370         176 :                     pasJobs[i].dfYOff = m_panToMSGT[3] +
    1371         176 :                                         (nYOff + iStartLine) * m_panToMSGT[5] -
    1372             :                                         nYOffExtract;
    1373         176 :                     pasJobs[i].dfXSize = sExtraArg.dfXSize;
    1374         176 :                     pasJobs[i].dfYSize =
    1375         176 :                         (iNextStartLine - iStartLine) * m_panToMSGT[5];
    1376         176 :                     if (pasJobs[i].dfXOff + pasJobs[i].dfXSize > nXSizeExtract)
    1377             :                     {
    1378           0 :                         pasJobs[i].dfXSize = nYSizeExtract - pasJobs[i].dfXOff;
    1379             :                     }
    1380         352 :                     if (pasJobs[i].dfYOff + pasJobs[i].dfYSize >
    1381         176 :                         aMSBands[0]->GetYSize())
    1382             :                     {
    1383           0 :                         pasJobs[i].dfYSize =
    1384           0 :                             aMSBands[0]->GetYSize() - pasJobs[i].dfYOff;
    1385             :                     }
    1386         176 :                     pasJobs[i].nXOff = static_cast<int>(pasJobs[i].dfXOff);
    1387         176 :                     pasJobs[i].nYOff = static_cast<int>(pasJobs[i].dfYOff);
    1388         176 :                     pasJobs[i].nXSize =
    1389         176 :                         static_cast<int>(0.4999 + pasJobs[i].dfXSize);
    1390         176 :                     pasJobs[i].nYSize =
    1391         176 :                         static_cast<int>(0.4999 + pasJobs[i].dfYSize);
    1392         176 :                     if (pasJobs[i].nXOff + pasJobs[i].nXSize > nXSizeExtract)
    1393             :                     {
    1394           0 :                         pasJobs[i].nXSize = nXSizeExtract - pasJobs[i].nXOff;
    1395             :                     }
    1396         176 :                     if (pasJobs[i].nYOff + pasJobs[i].nYSize > nYSizeExtract)
    1397             :                     {
    1398           0 :                         pasJobs[i].nYSize = nYSizeExtract - pasJobs[i].nYOff;
    1399             :                     }
    1400         176 :                     if (pasJobs[i].nXSize == 0)
    1401           0 :                         pasJobs[i].nXSize = 1;
    1402         176 :                     if (pasJobs[i].nYSize == 0)
    1403           0 :                         pasJobs[i].nYSize = 1;
    1404         176 :                     pasJobs[i].pBuffer = pUpsampledSpectralBuffer +
    1405         176 :                                          static_cast<size_t>(iStartLine) *
    1406         176 :                                              nXSize * nDataTypeSize;
    1407         176 :                     pasJobs[i].eDT = eWorkDataType;
    1408         176 :                     pasJobs[i].nBufXSize = nXSize;
    1409         176 :                     pasJobs[i].nBufYSize =
    1410         176 :                         static_cast<int>(iNextStartLine - iStartLine);
    1411         176 :                     pasJobs[i].nBandCount = psOptions->nInputSpectralBands;
    1412         176 :                     pasJobs[i].nBandSpace =
    1413         176 :                         static_cast<GSpacing>(nXSize) * nYSize * nDataTypeSize;
    1414             : #ifdef DEBUG_TIMING
    1415             :                     pasJobs[i].ptv = &tv;
    1416             : #endif
    1417         176 :                     pasJobs[i].eErr = CE_Failure;
    1418             : 
    1419         176 :                     ahJobData[i] = &(pasJobs[i]);
    1420             :                 }
    1421             : #ifdef DEBUG_TIMING
    1422             :                 gettimeofday(&tv, nullptr);
    1423             : #endif
    1424          44 :                 poThreadPool->SubmitJobs(PansharpenResampleJobThreadFunc,
    1425             :                                          ahJobData);
    1426          44 :                 poThreadPool->WaitCompletion();
    1427             : 
    1428         220 :                 for (int i = 0; i < nTasks; i++)
    1429             :                 {
    1430         176 :                     if (pasJobs[i].eErr == CE_Failure)
    1431             :                     {
    1432           0 :                         CPLError(CE_Failure, CPLE_AppDefined, "%s",
    1433           0 :                                  pasJobs[i].osLastErrorMsg.c_str());
    1434           0 :                         GDALClose(poMEMDS);
    1435           0 :                         VSIFree(pSpectralBuffer);
    1436           0 :                         VSIFree(pUpsampledSpectralBuffer);
    1437           0 :                         VSIFree(pPanBuffer);
    1438             : 
    1439           0 :                         return CE_Failure;
    1440             :                     }
    1441             :                 }
    1442             :             }
    1443             :         }
    1444             : 
    1445          81 :         GDALClose(poMEMDS);
    1446             : 
    1447          81 :         VSIFree(pSpectralBuffer);
    1448             :     }
    1449             :     else
    1450             :     {
    1451          10 :         if (!anInputBands.empty())
    1452             :         {
    1453             :             // Use dataset RasterIO when possible.
    1454          20 :             eErr = aMSBands[0]->GetDataset()->RasterIO(
    1455             :                 GF_Read, nSpectralXOff, nSpectralYOff, nSpectralXSize,
    1456             :                 nSpectralYSize, pUpsampledSpectralBuffer, nXSize, nYSize,
    1457          10 :                 eWorkDataType, static_cast<int>(anInputBands.size()),
    1458          10 :                 &anInputBands[0], 0, 0, 0, &sExtraArg);
    1459             :         }
    1460             :         else
    1461             :         {
    1462           0 :             for (int i = 0;
    1463           0 :                  eErr == CE_None && i < psOptions->nInputSpectralBands; i++)
    1464             :             {
    1465           0 :                 eErr = aMSBands[i]->RasterIO(
    1466             :                     GF_Read, nSpectralXOff, nSpectralYOff, nSpectralXSize,
    1467             :                     nSpectralYSize,
    1468           0 :                     pUpsampledSpectralBuffer + static_cast<size_t>(i) * nXSize *
    1469           0 :                                                    nYSize * nDataTypeSize,
    1470             :                     nXSize, nYSize, eWorkDataType, 0, 0, &sExtraArg);
    1471             :             }
    1472             :         }
    1473          10 :         if (eErr != CE_None)
    1474             :         {
    1475           0 :             VSIFree(pUpsampledSpectralBuffer);
    1476           0 :             VSIFree(pPanBuffer);
    1477           0 :             return CE_Failure;
    1478             :         }
    1479             :     }
    1480             : 
    1481             :     // In case NBITS was not set on the spectral bands, clamp the values
    1482             :     // if overshoot might have occurred.
    1483          91 :     int nBitDepth = psOptions->nBitDepth;
    1484          91 :     if (nBitDepth &&
    1485           0 :         (eResampleAlg == GRIORA_Cubic || eResampleAlg == GRIORA_CubicSpline ||
    1486             :          eResampleAlg == GRIORA_Lanczos))
    1487             :     {
    1488          12 :         for (int i = 0; i < psOptions->nInputSpectralBands; i++)
    1489             :         {
    1490           6 :             GDALRasterBand *poBand = aMSBands[i];
    1491           6 :             int nBandBitDepth = 0;
    1492             :             const char *pszNBITS =
    1493           6 :                 poBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE");
    1494           6 :             if (pszNBITS)
    1495           4 :                 nBandBitDepth = atoi(pszNBITS);
    1496           6 :             if (nBandBitDepth < nBitDepth)
    1497             :             {
    1498           2 :                 if (eWorkDataType == GDT_Byte && nBitDepth >= 0 &&
    1499             :                     nBitDepth <= 8)
    1500             :                 {
    1501           1 :                     ClampValues(
    1502             :                         reinterpret_cast<GByte *>(pUpsampledSpectralBuffer) +
    1503           1 :                             static_cast<size_t>(i) * nXSize * nYSize,
    1504           1 :                         static_cast<size_t>(nXSize) * nYSize,
    1505           1 :                         static_cast<GByte>((1 << nBitDepth) - 1));
    1506             :                 }
    1507           1 :                 else if (eWorkDataType == GDT_UInt16 && nBitDepth >= 0 &&
    1508             :                          nBitDepth <= 16)
    1509             :                 {
    1510           1 :                     ClampValues(
    1511           1 :                         reinterpret_cast<GUInt16 *>(pUpsampledSpectralBuffer) +
    1512           1 :                             static_cast<size_t>(i) * nXSize * nYSize,
    1513           1 :                         static_cast<size_t>(nXSize) * nYSize,
    1514           1 :                         static_cast<GUInt16>((1 << nBitDepth) - 1));
    1515             :                 }
    1516             : #ifndef LIMIT_TYPES
    1517             :                 else if (eWorkDataType == GDT_UInt32)
    1518             :                 {
    1519             :                     ClampValues(reinterpret_cast<GUInt32*>(pUpsampledSpectralBuffer) +
    1520             :                                 static_cast<size_t>(i) * nXSize * nYSize,
    1521             :                                 static_cast<size_t>(nXSize)*nYSize,
    1522             :                                 (static_cast<GUInt32>((1 << nBitDepth)-1));
    1523             :                 }
    1524             : #endif
    1525             :             }
    1526             :         }
    1527             :     }
    1528             : 
    1529          91 :     const GUInt32 nMaxValue = (nBitDepth >= 0 && nBitDepth <= 31)
    1530         182 :                                   ? (1U << nBitDepth) - 1
    1531             :                                   : UINT32_MAX;
    1532             : 
    1533          91 :     double *padfTempBuffer = nullptr;
    1534          91 :     GDALDataType eBufDataTypeOri = eBufDataType;
    1535          91 :     void *pDataBufOri = pDataBuf;
    1536             :     // CFloat64 is the query type used by gdallocationinfo...
    1537             : #ifdef LIMIT_TYPES
    1538          91 :     if (eBufDataType != GDT_Byte && eBufDataType != GDT_UInt16)
    1539             : #else
    1540             :     if (eBufDataType == GDT_CFloat64)
    1541             : #endif
    1542             :     {
    1543          50 :         padfTempBuffer = static_cast<double *>(VSI_MALLOC3_VERBOSE(
    1544             :             nXSize, nYSize, psOptions->nOutPansharpenedBands * sizeof(double)));
    1545          50 :         if (padfTempBuffer == nullptr)
    1546             :         {
    1547           0 :             VSIFree(pUpsampledSpectralBuffer);
    1548           0 :             VSIFree(pPanBuffer);
    1549           0 :             return CE_Failure;
    1550             :         }
    1551          50 :         pDataBuf = padfTempBuffer;
    1552          50 :         eBufDataType = GDT_Float64;
    1553             :     }
    1554             : 
    1555          91 :     if (nTasks > 1)
    1556             :     {
    1557          88 :         std::vector<GDALPansharpenJob> asJobs;
    1558          44 :         asJobs.resize(nTasks);
    1559          44 :         GDALPansharpenJob *pasJobs = &(asJobs[0]);
    1560             :         {
    1561          88 :             std::vector<void *> ahJobData;
    1562          44 :             ahJobData.resize(nTasks);
    1563             : #ifdef DEBUG_TIMING
    1564             :             struct timeval tv;
    1565             : #endif
    1566         220 :             for (int i = 0; i < nTasks; i++)
    1567             :             {
    1568         176 :                 const size_t iStartLine =
    1569         176 :                     (static_cast<size_t>(i) * nYSize) / nTasks;
    1570         176 :                 const size_t iNextStartLine =
    1571         176 :                     (static_cast<size_t>(i + 1) * nYSize) / nTasks;
    1572         176 :                 pasJobs[i].poPansharpenOperation = this;
    1573         176 :                 pasJobs[i].eWorkDataType = eWorkDataType;
    1574         176 :                 pasJobs[i].eBufDataType = eBufDataType;
    1575         176 :                 pasJobs[i].pPanBuffer =
    1576         176 :                     pPanBuffer + iStartLine * nXSize * nDataTypeSize;
    1577         176 :                 pasJobs[i].pUpsampledSpectralBuffer =
    1578         176 :                     pUpsampledSpectralBuffer +
    1579         176 :                     iStartLine * nXSize * nDataTypeSize;
    1580         176 :                 pasJobs[i].pDataBuf =
    1581         176 :                     static_cast<GByte *>(pDataBuf) +
    1582         352 :                     iStartLine * nXSize *
    1583         176 :                         GDALGetDataTypeSizeBytes(eBufDataType);
    1584         176 :                 pasJobs[i].nValues = (iNextStartLine - iStartLine) * nXSize;
    1585         176 :                 pasJobs[i].nBandValues = static_cast<size_t>(nXSize) * nYSize;
    1586         176 :                 pasJobs[i].nMaxValue = nMaxValue;
    1587             : #ifdef DEBUG_TIMING
    1588             :                 pasJobs[i].ptv = &tv;
    1589             : #endif
    1590         176 :                 ahJobData[i] = &(pasJobs[i]);
    1591             :             }
    1592             : #ifdef DEBUG_TIMING
    1593             :             gettimeofday(&tv, nullptr);
    1594             : #endif
    1595          44 :             poThreadPool->SubmitJobs(PansharpenJobThreadFunc, ahJobData);
    1596          44 :             poThreadPool->WaitCompletion();
    1597             :         }
    1598             : 
    1599          44 :         eErr = CE_None;
    1600         220 :         for (int i = 0; i < nTasks; i++)
    1601             :         {
    1602         176 :             if (pasJobs[i].eErr != CE_None)
    1603           0 :                 eErr = CE_Failure;
    1604             :         }
    1605             :     }
    1606             :     else
    1607             :     {
    1608          47 :         eErr = PansharpenChunk(eWorkDataType, eBufDataType, pPanBuffer,
    1609             :                                pUpsampledSpectralBuffer, pDataBuf,
    1610          47 :                                static_cast<size_t>(nXSize) * nYSize,
    1611          47 :                                static_cast<size_t>(nXSize) * nYSize, nMaxValue);
    1612             :     }
    1613             : 
    1614          91 :     if (padfTempBuffer)
    1615             :     {
    1616          50 :         GDALCopyWords64(padfTempBuffer, GDT_Float64, sizeof(double),
    1617             :                         pDataBufOri, eBufDataTypeOri,
    1618             :                         GDALGetDataTypeSizeBytes(eBufDataTypeOri),
    1619          50 :                         static_cast<size_t>(nXSize) * nYSize *
    1620          50 :                             psOptions->nOutPansharpenedBands);
    1621          50 :         VSIFree(padfTempBuffer);
    1622             :     }
    1623             : 
    1624          91 :     VSIFree(pUpsampledSpectralBuffer);
    1625          91 :     VSIFree(pPanBuffer);
    1626             : 
    1627          91 :     return eErr;
    1628             : }
    1629             : 
    1630             : /************************************************************************/
    1631             : /*                   PansharpenResampleJobThreadFunc()                  */
    1632             : /************************************************************************/
    1633             : 
    1634         174 : void GDALPansharpenOperation::PansharpenResampleJobThreadFunc(void *pUserData)
    1635             : {
    1636         174 :     GDALPansharpenResampleJob *psJob =
    1637             :         static_cast<GDALPansharpenResampleJob *>(pUserData);
    1638             : 
    1639             : #ifdef DEBUG_TIMING
    1640             :     struct timeval tv;
    1641             :     gettimeofday(&tv, nullptr);
    1642             :     const GIntBig launch_time =
    1643             :         static_cast<GIntBig>(psJob->ptv->tv_sec) * 1000000 +
    1644             :         static_cast<GIntBig>(psJob->ptv->tv_usec);
    1645             :     const GIntBig start_job = static_cast<GIntBig>(tv.tv_sec) * 1000000 +
    1646             :                               static_cast<GIntBig>(tv.tv_usec);
    1647             : #endif
    1648             : 
    1649             :     GDALRasterIOExtraArg sExtraArg;
    1650         174 :     INIT_RASTERIO_EXTRA_ARG(sExtraArg);
    1651             :     // cppcheck-suppress redundantAssignment
    1652         174 :     sExtraArg.eResampleAlg = psJob->eResampleAlg;
    1653         174 :     sExtraArg.bFloatingPointWindowValidity = TRUE;
    1654         174 :     sExtraArg.dfXOff = psJob->dfXOff;
    1655         174 :     sExtraArg.dfYOff = psJob->dfYOff;
    1656         174 :     sExtraArg.dfXSize = psJob->dfXSize;
    1657         174 :     sExtraArg.dfYSize = psJob->dfYSize;
    1658             : 
    1659         350 :     std::vector<int> anBands;
    1660         759 :     for (int i = 0; i < psJob->nBandCount; ++i)
    1661         583 :         anBands.push_back(i + 1);
    1662             :     // This call to RasterIO() in a thread to poMEMDS shared between several
    1663             :     // threads is really risky, but works given the implementation details...
    1664             :     // Do not do this at home!
    1665         173 :     psJob->eErr = psJob->poMEMDS->RasterIO(
    1666             :         GF_Read, psJob->nXOff, psJob->nYOff, psJob->nXSize, psJob->nYSize,
    1667             :         psJob->pBuffer, psJob->nBufXSize, psJob->nBufYSize, psJob->eDT,
    1668         176 :         psJob->nBandCount, anBands.data(), 0, 0, psJob->nBandSpace, &sExtraArg);
    1669         176 :     if (CPLGetLastErrorType() == CE_Failure)
    1670           0 :         psJob->osLastErrorMsg = CPLGetLastErrorMsg();
    1671             : 
    1672             : #ifdef DEBUG_TIMING
    1673             :     struct timeval tv_end;
    1674             :     gettimeofday(&tv_end, nullptr);
    1675             :     const GIntBig end = static_cast<GIntBig>(tv_end.tv_sec) * 1000000 +
    1676             :                         static_cast<GIntBig>(tv_end.tv_usec);
    1677             :     if (start_job - launch_time > 500)
    1678             :         /*ok*/ printf("Resample: Delay before start=" CPL_FRMT_GIB
    1679             :                       ", completion time=" CPL_FRMT_GIB "\n",
    1680             :                       start_job - launch_time, end - start_job);
    1681             : #endif
    1682         176 : }
    1683             : 
    1684             : /************************************************************************/
    1685             : /*                      PansharpenJobThreadFunc()                       */
    1686             : /************************************************************************/
    1687             : 
    1688         174 : void GDALPansharpenOperation::PansharpenJobThreadFunc(void *pUserData)
    1689             : {
    1690         174 :     GDALPansharpenJob *psJob = static_cast<GDALPansharpenJob *>(pUserData);
    1691             : 
    1692             : #ifdef DEBUG_TIMING
    1693             :     struct timeval tv;
    1694             :     gettimeofday(&tv, nullptr);
    1695             :     const GIntBig launch_time =
    1696             :         static_cast<GIntBig>(psJob->ptv->tv_sec) * 1000000 +
    1697             :         static_cast<GIntBig>(psJob->ptv->tv_usec);
    1698             :     const GIntBig start_job = static_cast<GIntBig>(tv.tv_sec) * 1000000 +
    1699             :                               static_cast<GIntBig>(tv.tv_usec);
    1700             : #endif
    1701             : 
    1702             : #if 0
    1703             :     for( int i = 0; i < 1000000; i++ )
    1704             :         acc += i * i;
    1705             :     psJob->eErr = CE_None;
    1706             : #else
    1707         174 :     psJob->eErr = psJob->poPansharpenOperation->PansharpenChunk(
    1708             :         psJob->eWorkDataType, psJob->eBufDataType, psJob->pPanBuffer,
    1709             :         psJob->pUpsampledSpectralBuffer, psJob->pDataBuf, psJob->nValues,
    1710             :         psJob->nBandValues, psJob->nMaxValue);
    1711             : #endif
    1712             : 
    1713             : #ifdef DEBUG_TIMING
    1714             :     struct timeval tv_end;
    1715             :     gettimeofday(&tv_end, nullptr);
    1716             :     const GIntBig end = static_cast<GIntBig>(tv_end.tv_sec) * 1000000 +
    1717             :                         static_cast<GIntBig>(tv_end.tv_usec);
    1718             :     if (start_job - launch_time > 500)
    1719             :         /*ok*/ printf("Pansharpen: Delay before start=" CPL_FRMT_GIB
    1720             :                       ", completion time=" CPL_FRMT_GIB "\n",
    1721             :                       start_job - launch_time, end - start_job);
    1722             : #endif
    1723         176 : }
    1724             : 
    1725             : /************************************************************************/
    1726             : /*                           PansharpenChunk()                          */
    1727             : /************************************************************************/
    1728             : 
    1729         221 : CPLErr GDALPansharpenOperation::PansharpenChunk(
    1730             :     GDALDataType eWorkDataType, GDALDataType eBufDataType,
    1731             :     const void *pPanBuffer, const void *pUpsampledSpectralBuffer,
    1732             :     void *pDataBuf, size_t nValues, size_t nBandValues, GUInt32 nMaxValue) const
    1733             : {
    1734         221 :     CPLErr eErr = CE_None;
    1735             : 
    1736         221 :     switch (eWorkDataType)
    1737             :     {
    1738         109 :         case GDT_Byte:
    1739         109 :             eErr = WeightedBrovey(
    1740             :                 static_cast<const GByte *>(pPanBuffer),
    1741             :                 static_cast<const GByte *>(pUpsampledSpectralBuffer), pDataBuf,
    1742             :                 eBufDataType, nValues, nBandValues,
    1743             :                 static_cast<GByte>(nMaxValue));
    1744         109 :             break;
    1745             : 
    1746         107 :         case GDT_UInt16:
    1747         107 :             eErr = WeightedBrovey(
    1748             :                 static_cast<const GUInt16 *>(pPanBuffer),
    1749             :                 static_cast<const GUInt16 *>(pUpsampledSpectralBuffer),
    1750             :                 pDataBuf, eBufDataType, nValues, nBandValues,
    1751             :                 static_cast<GUInt16>(nMaxValue));
    1752         108 :             break;
    1753             : 
    1754             : #ifndef LIMIT_TYPES
    1755             :         case GDT_Int8:
    1756             :             eErr = WeightedBrovey(
    1757             :                 static_cast<const GInt8 *>(pPanBuffer),
    1758             :                 static_cast<const GInt8 *>(pUpsampledSpectralBuffer), pDataBuf,
    1759             :                 eBufDataType, nValues, nBandValues);
    1760             :             break;
    1761             : 
    1762             :         case GDT_Int16:
    1763             :             eErr = WeightedBrovey(
    1764             :                 static_cast<const GInt16 *>(pPanBuffer),
    1765             :                 static_cast<const GInt16 *>(pUpsampledSpectralBuffer), pDataBuf,
    1766             :                 eBufDataType, nValues, nBandValues);
    1767             :             break;
    1768             : 
    1769             :         case GDT_UInt32:
    1770             :             eErr = WeightedBrovey(
    1771             :                 static_cast<const GUInt32 *>(pPanBuffer),
    1772             :                 static_cast<const GUInt32 *>(pUpsampledSpectralBuffer),
    1773             :                 pDataBuf, eBufDataType, nValues, nBandValues, nMaxValue);
    1774             :             break;
    1775             : 
    1776             :         case GDT_Int32:
    1777             :             eErr = WeightedBrovey(
    1778             :                 static_cast<const GInt32 *>(pPanBuffer),
    1779             :                 static_cast<const GInt32 *>(pUpsampledSpectralBuffer), pDataBuf,
    1780             :                 eBufDataType, nValues, nBandValues);
    1781             :             break;
    1782             : 
    1783             :         case GDT_UInt64:
    1784             :             eErr = WeightedBrovey(
    1785             :                 static_cast<const std::uint64_t *>(pPanBuffer),
    1786             :                 static_cast<const std::uint64_t *>(pUpsampledSpectralBuffer),
    1787             :                 pDataBuf, eBufDataType, nValues, nBandValues, nMaxValue);
    1788             :             break;
    1789             : 
    1790             :         case GDT_Int64:
    1791             :             eErr = WeightedBrovey(
    1792             :                 static_cast<const std::int64_t *>(pPanBuffer),
    1793             :                 static_cast<const std::int64_t *>(pUpsampledSpectralBuffer),
    1794             :                 pDataBuf, eBufDataType, nValues, nBandValues);
    1795             :             break;
    1796             : 
    1797             :         case GDT_Float16:
    1798             :             eErr = WeightedBrovey(
    1799             :                 static_cast<const GFloat16 *>(pPanBuffer),
    1800             :                 static_cast<const GFloat16 *>(pUpsampledSpectralBuffer),
    1801             :                 pDataBuf, eBufDataType, nValues, nBandValues);
    1802             :             break;
    1803             : 
    1804             :         case GDT_Float32:
    1805             :             eErr = WeightedBrovey(
    1806             :                 static_cast<const float *>(pPanBuffer),
    1807             :                 static_cast<const float *>(pUpsampledSpectralBuffer), pDataBuf,
    1808             :                 eBufDataType, nValues, nBandValues);
    1809             :             break;
    1810             : #endif
    1811           6 :         case GDT_Float64:
    1812           6 :             eErr = WeightedBrovey(
    1813             :                 static_cast<const double *>(pPanBuffer),
    1814             :                 static_cast<const double *>(pUpsampledSpectralBuffer), pDataBuf,
    1815             :                 eBufDataType, nValues, nBandValues);
    1816           6 :             break;
    1817             : 
    1818           0 :         default:
    1819           0 :             CPLError(CE_Failure, CPLE_NotSupported,
    1820             :                      "eWorkDataType not supported");
    1821           0 :             eErr = CE_Failure;
    1822           0 :             break;
    1823             :     }
    1824             : 
    1825         223 :     return eErr;
    1826             : }
    1827             : 
    1828             : /************************************************************************/
    1829             : /*                             GetOptions()                             */
    1830             : /************************************************************************/
    1831             : 
    1832             : /** Return options.
    1833             :  * @return options.
    1834             :  */
    1835         151 : GDALPansharpenOptions *GDALPansharpenOperation::GetOptions()
    1836             : {
    1837         151 :     return psOptions;
    1838             : }
    1839             : 
    1840             : /************************************************************************/
    1841             : /*                     GDALCreatePansharpenOperation()                  */
    1842             : /************************************************************************/
    1843             : 
    1844             : /** Instantiate a pansharpening operation.
    1845             :  *
    1846             :  * The passed options are validated.
    1847             :  *
    1848             :  * @param psOptions a pansharpening option structure allocated with
    1849             :  * GDALCreatePansharpenOptions(). It is duplicated by this function.
    1850             :  * @return a valid pansharpening operation handle, or NULL in case of failure.
    1851             :  *
    1852             :  */
    1853             : 
    1854             : GDALPansharpenOperationH
    1855           0 : GDALCreatePansharpenOperation(const GDALPansharpenOptions *psOptions)
    1856             : {
    1857           0 :     GDALPansharpenOperation *psOperation = new GDALPansharpenOperation();
    1858           0 :     if (psOperation->Initialize(psOptions) == CE_None)
    1859           0 :         return reinterpret_cast<GDALPansharpenOperationH>(psOperation);
    1860           0 :     delete psOperation;
    1861           0 :     return nullptr;
    1862             : }
    1863             : 
    1864             : /************************************************************************/
    1865             : /*                     GDALDestroyPansharpenOperation()                 */
    1866             : /************************************************************************/
    1867             : 
    1868             : /** Destroy a pansharpening operation.
    1869             :  *
    1870             :  * @param hOperation a valid pansharpening operation.
    1871             :  *
    1872             :  */
    1873             : 
    1874           0 : void GDALDestroyPansharpenOperation(GDALPansharpenOperationH hOperation)
    1875             : {
    1876           0 :     delete reinterpret_cast<GDALPansharpenOperation *>(hOperation);
    1877           0 : }
    1878             : 
    1879             : /************************************************************************/
    1880             : /*                       GDALPansharpenProcessRegion()                  */
    1881             : /************************************************************************/
    1882             : 
    1883             : /** Executes a pansharpening operation on a rectangular region of the
    1884             :  * resulting dataset.
    1885             :  *
    1886             :  * The window is expressed with respect to the dimensions of the panchromatic
    1887             :  * band.
    1888             :  *
    1889             :  * Spectral bands are upsampled and merged with the panchromatic band according
    1890             :  * to the select algorithm and options.
    1891             :  *
    1892             :  * @param hOperation a valid pansharpening operation.
    1893             :  * @param nXOff pixel offset.
    1894             :  * @param nYOff pixel offset.
    1895             :  * @param nXSize width of the pansharpened region to compute.
    1896             :  * @param nYSize height of the pansharpened region to compute.
    1897             :  * @param pDataBuf output buffer. Must be nXSize * nYSize *
    1898             :  *                 GDALGetDataTypeSizeBytes(eBufDataType) *
    1899             :  *                 psOptions->nOutPansharpenedBands large.
    1900             :  *                 It begins with all values of the first output band, followed
    1901             :  *                 by values of the second output band, etc...
    1902             :  * @param eBufDataType data type of the output buffer
    1903             :  *
    1904             :  * @return CE_None in case of success, CE_Failure in case of failure.
    1905             :  *
    1906             :  */
    1907           0 : CPLErr GDALPansharpenProcessRegion(GDALPansharpenOperationH hOperation,
    1908             :                                    int nXOff, int nYOff, int nXSize, int nYSize,
    1909             :                                    void *pDataBuf, GDALDataType eBufDataType)
    1910             : {
    1911             :     return reinterpret_cast<GDALPansharpenOperation *>(hOperation)
    1912           0 :         ->ProcessRegion(nXOff, nYOff, nXSize, nYSize, pDataBuf, eBufDataType);
    1913             : }

Generated by: LCOV version 1.14