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