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