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

Generated by: LCOV version 1.14