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 145 : GDALPansharpenOptions *GDALCreatePansharpenOptions()
54 : {
55 : GDALPansharpenOptions *psOptions = static_cast<GDALPansharpenOptions *>(
56 145 : CPLCalloc(1, sizeof(GDALPansharpenOptions)));
57 145 : psOptions->ePansharpenAlg = GDAL_PSH_WEIGHTED_BROVEY;
58 145 : psOptions->eResampleAlg = GRIORA_Cubic;
59 145 : 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 146 : void GDALDestroyPansharpenOptions(GDALPansharpenOptions *psOptions)
74 : {
75 146 : if (psOptions == nullptr)
76 1 : return;
77 145 : CPLFree(psOptions->padfWeights);
78 145 : CPLFree(psOptions->pahInputSpectralBands);
79 145 : CPLFree(psOptions->panOutPansharpenedBands);
80 145 : 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 75 : GDALClonePansharpenOptions(const GDALPansharpenOptions *psOptions)
98 : {
99 75 : GDALPansharpenOptions *psNewOptions = GDALCreatePansharpenOptions();
100 75 : psNewOptions->ePansharpenAlg = psOptions->ePansharpenAlg;
101 75 : psNewOptions->eResampleAlg = psOptions->eResampleAlg;
102 75 : psNewOptions->nBitDepth = psOptions->nBitDepth;
103 75 : psNewOptions->nWeightCount = psOptions->nWeightCount;
104 75 : if (psOptions->padfWeights)
105 : {
106 75 : psNewOptions->padfWeights = static_cast<double *>(
107 75 : CPLMalloc(sizeof(double) * psOptions->nWeightCount));
108 75 : memcpy(psNewOptions->padfWeights, psOptions->padfWeights,
109 75 : sizeof(double) * psOptions->nWeightCount);
110 : }
111 75 : psNewOptions->hPanchroBand = psOptions->hPanchroBand;
112 75 : psNewOptions->nInputSpectralBands = psOptions->nInputSpectralBands;
113 75 : if (psOptions->pahInputSpectralBands)
114 : {
115 75 : const size_t nSize =
116 75 : sizeof(GDALRasterBandH) * psOptions->nInputSpectralBands;
117 75 : psNewOptions->pahInputSpectralBands =
118 75 : static_cast<GDALRasterBandH *>(CPLMalloc(nSize));
119 75 : memcpy(psNewOptions->pahInputSpectralBands,
120 75 : psOptions->pahInputSpectralBands, nSize);
121 : }
122 75 : psNewOptions->nOutPansharpenedBands = psOptions->nOutPansharpenedBands;
123 75 : if (psOptions->panOutPansharpenedBands)
124 : {
125 74 : psNewOptions->panOutPansharpenedBands = static_cast<int *>(
126 74 : CPLMalloc(sizeof(int) * psOptions->nOutPansharpenedBands));
127 74 : memcpy(psNewOptions->panOutPansharpenedBands,
128 74 : psOptions->panOutPansharpenedBands,
129 74 : sizeof(int) * psOptions->nOutPansharpenedBands);
130 : }
131 75 : psNewOptions->bHasNoData = psOptions->bHasNoData;
132 75 : psNewOptions->dfNoData = psOptions->dfNoData;
133 75 : psNewOptions->nThreads = psOptions->nThreads;
134 75 : 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 73 : GDALPansharpenOperation::~GDALPansharpenOperation()
155 : {
156 73 : GDALDestroyPansharpenOptions(psOptions);
157 82 : for (size_t i = 0; i < aVDS.size(); i++)
158 9 : delete aVDS[i];
159 73 : delete poThreadPool;
160 73 : }
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 73 : GDALPansharpenOperation::Initialize(const GDALPansharpenOptions *psOptionsIn)
174 : {
175 73 : if (psOptionsIn->hPanchroBand == nullptr)
176 : {
177 0 : CPLError(CE_Failure, CPLE_AppDefined, "hPanchroBand not set");
178 0 : return CE_Failure;
179 : }
180 73 : 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 73 : if (psOptionsIn->padfWeights == nullptr ||
187 73 : 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 73 : auto poPanchroBand = GDALRasterBand::FromHandle(psOptionsIn->hPanchroBand);
196 73 : auto poPanchroDS = poPanchroBand->GetDataset();
197 73 : 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 73 : 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 73 : GDALGeoTransform panchroGT;
212 73 : 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 73 : GDALRasterBandH hRefBand = psOptionsIn->pahInputSpectralBands[0];
220 73 : auto poRefBand = GDALRasterBand::FromHandle(hRefBand);
221 73 : auto poRefBandDS = poRefBand->GetDataset();
222 73 : 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 73 : 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 73 : GDALGeoTransform refMSGT;
239 73 : 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 73 : GDALGeoTransform invMSGT;
247 73 : 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 73 : m_panToMSGT[1] = invMSGT[1] * panchroGT[1] + invMSGT[2] * panchroGT[4];
256 73 : m_panToMSGT[2] = invMSGT[1] * panchroGT[2] + invMSGT[2] * panchroGT[5];
257 146 : m_panToMSGT[0] =
258 73 : invMSGT[1] * panchroGT[0] + invMSGT[2] * panchroGT[3] + invMSGT[0];
259 73 : m_panToMSGT[4] = invMSGT[4] * panchroGT[1] + invMSGT[5] * panchroGT[4];
260 73 : m_panToMSGT[5] = invMSGT[4] * panchroGT[2] + invMSGT[5] * panchroGT[5];
261 146 : m_panToMSGT[3] =
262 73 : 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 73 : 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 73 : bool bSameDataset = psOptionsIn->nInputSpectralBands > 1;
277 73 : if (bSameDataset)
278 55 : anInputBands.push_back(GDALGetBandNumber(hRefBand));
279 177 : 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 72 : if (psOptionsIn->nOutPansharpenedBands == 0)
342 : {
343 1 : CPLError(CE_Warning, CPLE_AppDefined,
344 : "No output pansharpened band defined");
345 : }
346 242 : for (int i = 0; i < psOptionsIn->nOutPansharpenedBands; i++)
347 : {
348 170 : if (psOptionsIn->panOutPansharpenedBands[i] < 0 ||
349 170 : psOptionsIn->panOutPansharpenedBands[i] >=
350 170 : 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 72 : GDALDataType eWorkDataType = poPanchroBand->GetRasterDataType();
360 72 : 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 72 : psOptions = GDALClonePansharpenOptions(psOptionsIn);
375 72 : if (psOptions->nBitDepth == GDALGetDataTypeSizeBits(eWorkDataType))
376 7 : psOptions->nBitDepth = 0;
377 72 : 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 248 : for (int i = 0; i < psOptions->nInputSpectralBands; i++)
389 : {
390 176 : if (psOptions->padfWeights[i] < 0.0)
391 : {
392 0 : bPositiveWeights = FALSE;
393 0 : break;
394 : }
395 : }
396 :
397 248 : for (int i = 0; i < psOptions->nInputSpectralBands; i++)
398 : {
399 352 : aMSBands.push_back(
400 176 : GDALRasterBand::FromHandle(psOptions->pahInputSpectralBands[i]));
401 : }
402 :
403 72 : 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 : const char *pszNBITS =
444 17 : poSrcBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE");
445 17 : if (pszNBITS)
446 0 : poVRTBand->SetMetadataItem("NBITS", pszNBITS,
447 0 : "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 72 : int nThreads = psOptions->nThreads;
461 72 : if (nThreads == -1)
462 26 : nThreads = CPLGetNumCPUs();
463 46 : else if (nThreads == 0)
464 : {
465 46 : nThreads = GDALGetNumThreads(GDAL_DEFAULT_MAX_THREAD_COUNT,
466 : /* bDefaultAllCPUs = */ false);
467 : }
468 72 : if (nThreads > 1)
469 : {
470 26 : CPLDebug("PANSHARPEN", "Using %d threads", nThreads);
471 52 : poThreadPool = new (std::nothrow) CPLWorkerThreadPool();
472 : // coverity[tainted_data]
473 52 : if (poThreadPool == nullptr ||
474 26 : !poThreadPool->Setup(nThreads, nullptr, nullptr))
475 : {
476 0 : delete poThreadPool;
477 0 : poThreadPool = nullptr;
478 : }
479 : }
480 :
481 72 : GDALRIOResampleAlg eResampleAlg = psOptions->eResampleAlg;
482 72 : if (eResampleAlg != GRIORA_NearestNeighbour)
483 : {
484 72 : const char *pszResampling =
485 144 : (eResampleAlg == GRIORA_Bilinear) ? "BILINEAR"
486 72 : : (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 72 : GDALGetResampleFunction(pszResampling, &nKernelRadius);
496 : }
497 :
498 72 : 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 5170 : 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 23964960 : static inline double ComputeFactor(T panValue, double dfPseudoPanchro)
568 : {
569 23964960 : if (dfPseudoPanchro == 0.0)
570 4163 : return 0.0;
571 :
572 23960800 : return panValue / dfPseudoPanchro;
573 : }
574 :
575 : /************************************************************************/
576 : /* ClampAndRound() */
577 : /************************************************************************/
578 :
579 2000138 : template <class T> static inline T ClampAndRound(double dfVal, T nMaxValue)
580 : {
581 2000138 : if (dfVal > nMaxValue)
582 50 : return nMaxValue;
583 : else
584 2000084 : 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 63 : void GDALPansharpenOperation::WeightedBroveyPositiveWeights(
778 : const T *pPanBuffer, const T *pUpsampledSpectralBuffer, T *pDataBuf,
779 : size_t nValues, size_t nBandValues, T nMaxValue) const
780 : {
781 63 : if (psOptions->bHasNoData)
782 : {
783 0 : WeightedBroveyWithNoData<T, T>(pPanBuffer, pUpsampledSpectralBuffer,
784 : pDataBuf, nValues, nBandValues,
785 : nMaxValue);
786 0 : return;
787 : }
788 :
789 63 : if (nMaxValue == 0)
790 57 : nMaxValue = cpl::NumericLimits<T>::max();
791 : size_t j;
792 63 : 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 30 : 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 26 : 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 500070 : for (j = 0; j + 1 < nValues; j += 2)
826 : {
827 500048 : double dfFactor = 0.0;
828 500048 : double dfFactor2 = 0.0;
829 500048 : double dfPseudoPanchro = 0.0;
830 500048 : double dfPseudoPanchro2 = 0.0;
831 1500098 : for (int i = 0; i < psOptions->nInputSpectralBands; i++)
832 : {
833 1000044 : dfPseudoPanchro +=
834 1000044 : psOptions->padfWeights[i] *
835 1000044 : pUpsampledSpectralBuffer[i * nBandValues + j];
836 1000044 : dfPseudoPanchro2 +=
837 1000044 : psOptions->padfWeights[i] *
838 1000044 : pUpsampledSpectralBuffer[i * nBandValues + j + 1];
839 : }
840 :
841 500048 : dfFactor = ComputeFactor(pPanBuffer[j], dfPseudoPanchro);
842 500048 : dfFactor2 = ComputeFactor(pPanBuffer[j + 1], dfPseudoPanchro2);
843 :
844 1500098 : for (int i = 0; i < psOptions->nOutPansharpenedBands; i++)
845 : {
846 1000044 : const T nRawValue = pUpsampledSpectralBuffer
847 1000044 : [psOptions->panOutPansharpenedBands[i] * nBandValues + j];
848 1000044 : const double dfTmp = nRawValue * dfFactor;
849 1000044 : pDataBuf[i * nBandValues + j] = ClampAndRound(dfTmp, nMaxValue);
850 :
851 1000044 : const T nRawValue2 = pUpsampledSpectralBuffer
852 1000044 : [psOptions->panOutPansharpenedBands[i] * nBandValues + j +
853 : 1];
854 1000044 : const double dfTmp2 = nRawValue2 * dfFactor2;
855 1000044 : pDataBuf[i * nBandValues + j + 1] =
856 1000044 : ClampAndRound(dfTmp2, nMaxValue);
857 : }
858 : }
859 : }
860 76 : 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 63 : void GDALPansharpenOperation::WeightedBroveyGByteOrUInt16(
901 : const T *pPanBuffer, const T *pUpsampledSpectralBuffer, T *pDataBuf,
902 : size_t nValues, size_t nBandValues, T nMaxValue) const
903 : {
904 63 : if (bPositiveWeights)
905 : {
906 63 : 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 63 : }
921 :
922 : template <>
923 48 : 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 48 : WeightedBroveyGByteOrUInt16(pPanBuffer, pUpsampledSpectralBuffer, pDataBuf,
928 : nValues, nBandValues, nMaxValue);
929 48 : }
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 217 : 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 217 : switch (eBufDataType)
949 : {
950 49 : case GDT_UInt8:
951 49 : WeightedBrovey(pPanBuffer, pUpsampledSpectralBuffer,
952 : static_cast<GByte *>(pDataBuf), nValues, nBandValues,
953 : nMaxValue);
954 49 : 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 217 : 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 91 : CPLErr GDALPansharpenOperation::ProcessRegion(int nXOff, int nYOff, int nXSize,
1157 : int nYSize, void *pDataBuf,
1158 : GDALDataType eBufDataType)
1159 : {
1160 91 : if (psOptions == nullptr)
1161 0 : return CE_Failure;
1162 :
1163 : // TODO: Avoid allocating buffers each time.
1164 : GDALRasterBand *poPanchroBand =
1165 91 : GDALRasterBand::FromHandle(psOptions->hPanchroBand);
1166 91 : GDALDataType eWorkDataType = poPanchroBand->GetRasterDataType();
1167 : #ifdef LIMIT_TYPES
1168 91 : if (eWorkDataType != GDT_UInt8 && eWorkDataType != GDT_UInt16)
1169 6 : eWorkDataType = GDT_Float64;
1170 : #endif
1171 91 : const int nDataTypeSize = GDALGetDataTypeSizeBytes(eWorkDataType);
1172 91 : 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 91 : VSI_MALLOC3_VERBOSE(nXSize, nYSize, nDataTypeSize));
1177 91 : if (pUpsampledSpectralBuffer == nullptr || pPanBuffer == nullptr)
1178 : {
1179 0 : VSIFree(pUpsampledSpectralBuffer);
1180 0 : VSIFree(pPanBuffer);
1181 0 : return CE_Failure;
1182 : }
1183 :
1184 91 : CPLErr eErr = poPanchroBand->RasterIO(GF_Read, nXOff, nYOff, nXSize, nYSize,
1185 : pPanBuffer, nXSize, nYSize,
1186 : eWorkDataType, 0, 0, nullptr);
1187 91 : if (eErr != CE_None)
1188 : {
1189 0 : VSIFree(pUpsampledSpectralBuffer);
1190 0 : VSIFree(pPanBuffer);
1191 0 : return CE_Failure;
1192 : }
1193 :
1194 91 : int nTasks = 0;
1195 91 : if (poThreadPool)
1196 : {
1197 44 : nTasks = poThreadPool->GetThreadCount();
1198 44 : if (nTasks > nYSize)
1199 0 : nTasks = nYSize;
1200 : }
1201 :
1202 : GDALRasterIOExtraArg sExtraArg;
1203 91 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
1204 91 : const GDALRIOResampleAlg eResampleAlg = psOptions->eResampleAlg;
1205 : // cppcheck-suppress redundantAssignment
1206 91 : sExtraArg.eResampleAlg = eResampleAlg;
1207 91 : sExtraArg.bFloatingPointWindowValidity = TRUE;
1208 91 : sExtraArg.dfXOff = m_panToMSGT[0] + nXOff * m_panToMSGT[1];
1209 91 : sExtraArg.dfYOff = m_panToMSGT[3] + nYOff * m_panToMSGT[5];
1210 91 : sExtraArg.dfXSize = nXSize * m_panToMSGT[1];
1211 91 : sExtraArg.dfYSize = nYSize * m_panToMSGT[5];
1212 91 : if (sExtraArg.dfXOff + sExtraArg.dfXSize > aMSBands[0]->GetXSize())
1213 2 : sExtraArg.dfXSize = aMSBands[0]->GetXSize() - sExtraArg.dfXOff;
1214 91 : if (sExtraArg.dfYOff + sExtraArg.dfYSize > aMSBands[0]->GetYSize())
1215 2 : sExtraArg.dfYSize = aMSBands[0]->GetYSize() - sExtraArg.dfYOff;
1216 91 : int nSpectralXOff = static_cast<int>(sExtraArg.dfXOff);
1217 91 : int nSpectralYOff = static_cast<int>(sExtraArg.dfYOff);
1218 91 : int nSpectralXSize = static_cast<int>(0.49999 + sExtraArg.dfXSize);
1219 91 : int nSpectralYSize = static_cast<int>(0.49999 + sExtraArg.dfYSize);
1220 91 : if (nSpectralXOff + nSpectralXSize > aMSBands[0]->GetXSize())
1221 0 : nSpectralXSize = aMSBands[0]->GetXSize() - nSpectralXOff;
1222 91 : if (nSpectralYOff + nSpectralYSize > aMSBands[0]->GetYSize())
1223 0 : nSpectralYSize = aMSBands[0]->GetYSize() - nSpectralYOff;
1224 91 : if (nSpectralXSize == 0)
1225 10 : nSpectralXSize = 1;
1226 91 : if (nSpectralYSize == 0)
1227 10 : nSpectralYSize = 1;
1228 :
1229 : // When upsampling, extract the multispectral data at
1230 : // full resolution in a temp buffer, and then do the upsampling.
1231 91 : if (nSpectralXSize < nXSize && nSpectralYSize < nYSize &&
1232 81 : eResampleAlg != GRIORA_NearestNeighbour && nYSize > 1)
1233 : {
1234 : // Take some margin to take into account the radius of the
1235 : // resampling kernel.
1236 81 : int nXOffExtract = nSpectralXOff - nKernelRadius;
1237 81 : int nYOffExtract = nSpectralYOff - nKernelRadius;
1238 81 : int nXSizeExtract = nSpectralXSize + 1 + 2 * nKernelRadius;
1239 81 : int nYSizeExtract = nSpectralYSize + 1 + 2 * nKernelRadius;
1240 81 : if (nXOffExtract < 0)
1241 : {
1242 74 : nXSizeExtract += nXOffExtract;
1243 74 : nXOffExtract = 0;
1244 : }
1245 81 : if (nYOffExtract < 0)
1246 : {
1247 67 : nYSizeExtract += nYOffExtract;
1248 67 : nYOffExtract = 0;
1249 : }
1250 81 : if (nXOffExtract + nXSizeExtract > aMSBands[0]->GetXSize())
1251 74 : nXSizeExtract = aMSBands[0]->GetXSize() - nXOffExtract;
1252 81 : if (nYOffExtract + nYSizeExtract > aMSBands[0]->GetYSize())
1253 67 : nYSizeExtract = aMSBands[0]->GetYSize() - nYOffExtract;
1254 :
1255 81 : GByte *pSpectralBuffer = static_cast<GByte *>(VSI_MALLOC3_VERBOSE(
1256 : nXSizeExtract, nYSizeExtract,
1257 : cpl::fits_on<int>(psOptions->nInputSpectralBands * nDataTypeSize)));
1258 81 : if (pSpectralBuffer == nullptr)
1259 : {
1260 0 : VSIFree(pUpsampledSpectralBuffer);
1261 0 : VSIFree(pPanBuffer);
1262 0 : return CE_Failure;
1263 : }
1264 :
1265 81 : 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 8 : for (int i = 0;
1277 20 : eErr == CE_None && i < psOptions->nInputSpectralBands; i++)
1278 : {
1279 24 : eErr = aMSBands[i]->RasterIO(
1280 : GF_Read, nXOffExtract, nYOffExtract, nXSizeExtract,
1281 : nYSizeExtract,
1282 12 : pSpectralBuffer + static_cast<size_t>(i) * nXSizeExtract *
1283 12 : nYSizeExtract * nDataTypeSize,
1284 : nXSizeExtract, nYSizeExtract, eWorkDataType, 0, 0, nullptr);
1285 : }
1286 : }
1287 81 : 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 81 : auto poMEMDS = MEMDataset::Create("", nXSizeExtract, nYSizeExtract, 0,
1297 : eWorkDataType, nullptr);
1298 :
1299 322 : for (int i = 0; i < psOptions->nInputSpectralBands; i++)
1300 : {
1301 241 : GByte *pabyBuffer =
1302 241 : pSpectralBuffer + static_cast<size_t>(i) * nDataTypeSize *
1303 241 : nXSizeExtract * nYSizeExtract;
1304 241 : GDALRasterBandH hMEMBand = MEMCreateRasterBandEx(
1305 : poMEMDS, i + 1, pabyBuffer, eWorkDataType, 0, 0, false);
1306 241 : poMEMDS->AddMEMBand(hMEMBand);
1307 :
1308 : const char *pszNBITS =
1309 241 : aMSBands[i]->GetMetadataItem("NBITS", "IMAGE_STRUCTURE");
1310 241 : if (pszNBITS)
1311 4 : poMEMDS->GetRasterBand(i + 1)->SetMetadataItem(
1312 4 : "NBITS", pszNBITS, "IMAGE_STRUCTURE");
1313 :
1314 241 : if (psOptions->bHasNoData)
1315 13 : poMEMDS->GetRasterBand(i + 1)->SetNoDataValue(
1316 13 : psOptions->dfNoData);
1317 : }
1318 :
1319 81 : 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 190 : for (int i = 0; i < poMEMDS->GetRasterCount(); i++)
1341 : {
1342 146 : poMEMDS->GetRasterBand(i + 1)->GetMaskFlags();
1343 : }
1344 :
1345 44 : std::vector<GDALPansharpenResampleJob> asJobs;
1346 44 : asJobs.resize(nTasks);
1347 44 : GDALPansharpenResampleJob *pasJobs = &(asJobs[0]);
1348 : {
1349 44 : std::vector<void *> ahJobData;
1350 44 : ahJobData.resize(nTasks);
1351 :
1352 : #ifdef DEBUG_TIMING
1353 : struct timeval tv;
1354 : #endif
1355 220 : for (int i = 0; i < nTasks; i++)
1356 : {
1357 176 : const size_t iStartLine =
1358 176 : (static_cast<size_t>(i) * nYSize) / nTasks;
1359 176 : const size_t iNextStartLine =
1360 176 : (static_cast<size_t>(i + 1) * nYSize) / nTasks;
1361 176 : pasJobs[i].poMEMDS = poMEMDS;
1362 176 : pasJobs[i].eResampleAlg = eResampleAlg;
1363 176 : pasJobs[i].dfXOff = sExtraArg.dfXOff - nXOffExtract;
1364 176 : pasJobs[i].dfYOff = m_panToMSGT[3] +
1365 176 : (nYOff + iStartLine) * m_panToMSGT[5] -
1366 : nYOffExtract;
1367 176 : pasJobs[i].dfXSize = sExtraArg.dfXSize;
1368 176 : pasJobs[i].dfYSize =
1369 176 : (iNextStartLine - iStartLine) * m_panToMSGT[5];
1370 176 : if (pasJobs[i].dfXOff + pasJobs[i].dfXSize > nXSizeExtract)
1371 : {
1372 0 : pasJobs[i].dfXSize = nYSizeExtract - pasJobs[i].dfXOff;
1373 : }
1374 352 : if (pasJobs[i].dfYOff + pasJobs[i].dfYSize >
1375 176 : aMSBands[0]->GetYSize())
1376 : {
1377 0 : pasJobs[i].dfYSize =
1378 0 : aMSBands[0]->GetYSize() - pasJobs[i].dfYOff;
1379 : }
1380 176 : pasJobs[i].nXOff = static_cast<int>(pasJobs[i].dfXOff);
1381 176 : pasJobs[i].nYOff = static_cast<int>(pasJobs[i].dfYOff);
1382 176 : pasJobs[i].nXSize =
1383 176 : static_cast<int>(0.4999 + pasJobs[i].dfXSize);
1384 176 : pasJobs[i].nYSize =
1385 176 : static_cast<int>(0.4999 + pasJobs[i].dfYSize);
1386 176 : if (pasJobs[i].nXOff + pasJobs[i].nXSize > nXSizeExtract)
1387 : {
1388 0 : pasJobs[i].nXSize = nXSizeExtract - pasJobs[i].nXOff;
1389 : }
1390 176 : if (pasJobs[i].nYOff + pasJobs[i].nYSize > nYSizeExtract)
1391 : {
1392 0 : pasJobs[i].nYSize = nYSizeExtract - pasJobs[i].nYOff;
1393 : }
1394 176 : if (pasJobs[i].nXSize == 0)
1395 0 : pasJobs[i].nXSize = 1;
1396 176 : if (pasJobs[i].nYSize == 0)
1397 0 : pasJobs[i].nYSize = 1;
1398 176 : pasJobs[i].pBuffer = pUpsampledSpectralBuffer +
1399 176 : static_cast<size_t>(iStartLine) *
1400 176 : nXSize * nDataTypeSize;
1401 176 : pasJobs[i].eDT = eWorkDataType;
1402 176 : pasJobs[i].nBufXSize = nXSize;
1403 176 : pasJobs[i].nBufYSize =
1404 176 : static_cast<int>(iNextStartLine - iStartLine);
1405 176 : pasJobs[i].nBandCount = psOptions->nInputSpectralBands;
1406 176 : pasJobs[i].nBandSpace =
1407 176 : static_cast<GSpacing>(nXSize) * nYSize * nDataTypeSize;
1408 : #ifdef DEBUG_TIMING
1409 : pasJobs[i].ptv = &tv;
1410 : #endif
1411 176 : pasJobs[i].eErr = CE_Failure;
1412 :
1413 176 : ahJobData[i] = &(pasJobs[i]);
1414 : }
1415 : #ifdef DEBUG_TIMING
1416 : gettimeofday(&tv, nullptr);
1417 : #endif
1418 44 : poThreadPool->SubmitJobs(PansharpenResampleJobThreadFunc,
1419 : ahJobData);
1420 44 : poThreadPool->WaitCompletion();
1421 :
1422 220 : for (int i = 0; i < nTasks; i++)
1423 : {
1424 176 : if (pasJobs[i].eErr == CE_Failure)
1425 : {
1426 0 : CPLError(CE_Failure, CPLE_AppDefined, "%s",
1427 0 : pasJobs[i].osLastErrorMsg.c_str());
1428 0 : GDALClose(poMEMDS);
1429 0 : VSIFree(pSpectralBuffer);
1430 0 : VSIFree(pUpsampledSpectralBuffer);
1431 0 : VSIFree(pPanBuffer);
1432 :
1433 0 : return CE_Failure;
1434 : }
1435 : }
1436 : }
1437 : }
1438 :
1439 81 : GDALClose(poMEMDS);
1440 :
1441 81 : VSIFree(pSpectralBuffer);
1442 : }
1443 : else
1444 : {
1445 10 : if (!anInputBands.empty())
1446 : {
1447 : // Use dataset RasterIO when possible.
1448 20 : eErr = aMSBands[0]->GetDataset()->RasterIO(
1449 : GF_Read, nSpectralXOff, nSpectralYOff, nSpectralXSize,
1450 : nSpectralYSize, pUpsampledSpectralBuffer, nXSize, nYSize,
1451 10 : eWorkDataType, static_cast<int>(anInputBands.size()),
1452 10 : &anInputBands[0], 0, 0, 0, &sExtraArg);
1453 : }
1454 : else
1455 : {
1456 0 : for (int i = 0;
1457 0 : eErr == CE_None && i < psOptions->nInputSpectralBands; i++)
1458 : {
1459 0 : eErr = aMSBands[i]->RasterIO(
1460 : GF_Read, nSpectralXOff, nSpectralYOff, nSpectralXSize,
1461 : nSpectralYSize,
1462 0 : pUpsampledSpectralBuffer + static_cast<size_t>(i) * nXSize *
1463 0 : nYSize * nDataTypeSize,
1464 : nXSize, nYSize, eWorkDataType, 0, 0, &sExtraArg);
1465 : }
1466 : }
1467 10 : if (eErr != CE_None)
1468 : {
1469 0 : VSIFree(pUpsampledSpectralBuffer);
1470 0 : VSIFree(pPanBuffer);
1471 0 : return CE_Failure;
1472 : }
1473 : }
1474 :
1475 : // In case NBITS was not set on the spectral bands, clamp the values
1476 : // if overshoot might have occurred.
1477 91 : int nBitDepth = psOptions->nBitDepth;
1478 91 : if (nBitDepth &&
1479 0 : (eResampleAlg == GRIORA_Cubic || eResampleAlg == GRIORA_CubicSpline ||
1480 : eResampleAlg == GRIORA_Lanczos))
1481 : {
1482 12 : for (int i = 0; i < psOptions->nInputSpectralBands; i++)
1483 : {
1484 6 : GDALRasterBand *poBand = aMSBands[i];
1485 6 : int nBandBitDepth = 0;
1486 : const char *pszNBITS =
1487 6 : poBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE");
1488 6 : if (pszNBITS)
1489 4 : nBandBitDepth = atoi(pszNBITS);
1490 6 : if (nBandBitDepth < nBitDepth)
1491 : {
1492 2 : if (eWorkDataType == GDT_UInt8 && nBitDepth >= 0 &&
1493 : nBitDepth <= 8)
1494 : {
1495 1 : ClampValues(
1496 : reinterpret_cast<GByte *>(pUpsampledSpectralBuffer) +
1497 1 : static_cast<size_t>(i) * nXSize * nYSize,
1498 1 : static_cast<size_t>(nXSize) * nYSize,
1499 1 : static_cast<GByte>((1 << nBitDepth) - 1));
1500 : }
1501 1 : else if (eWorkDataType == GDT_UInt16 && nBitDepth >= 0 &&
1502 : nBitDepth <= 16)
1503 : {
1504 1 : ClampValues(
1505 1 : reinterpret_cast<GUInt16 *>(pUpsampledSpectralBuffer) +
1506 1 : static_cast<size_t>(i) * nXSize * nYSize,
1507 1 : static_cast<size_t>(nXSize) * nYSize,
1508 1 : static_cast<GUInt16>((1 << nBitDepth) - 1));
1509 : }
1510 : #ifndef LIMIT_TYPES
1511 : else if (eWorkDataType == GDT_UInt32)
1512 : {
1513 : ClampValues(reinterpret_cast<GUInt32*>(pUpsampledSpectralBuffer) +
1514 : static_cast<size_t>(i) * nXSize * nYSize,
1515 : static_cast<size_t>(nXSize)*nYSize,
1516 : (static_cast<GUInt32>((1 << nBitDepth)-1));
1517 : }
1518 : #endif
1519 : }
1520 : }
1521 : }
1522 :
1523 91 : const GUInt32 nMaxValue = (nBitDepth >= 0 && nBitDepth <= 31)
1524 182 : ? (1U << nBitDepth) - 1
1525 : : UINT32_MAX;
1526 :
1527 91 : double *padfTempBuffer = nullptr;
1528 91 : GDALDataType eBufDataTypeOri = eBufDataType;
1529 91 : void *pDataBufOri = pDataBuf;
1530 : // CFloat64 is the query type used by gdallocationinfo...
1531 : #ifdef LIMIT_TYPES
1532 91 : if (eBufDataType != GDT_UInt8 && eBufDataType != GDT_UInt16)
1533 : #else
1534 : if (eBufDataType == GDT_CFloat64)
1535 : #endif
1536 : {
1537 50 : padfTempBuffer = static_cast<double *>(VSI_MALLOC3_VERBOSE(
1538 : nXSize, nYSize, psOptions->nOutPansharpenedBands * sizeof(double)));
1539 50 : if (padfTempBuffer == nullptr)
1540 : {
1541 0 : VSIFree(pUpsampledSpectralBuffer);
1542 0 : VSIFree(pPanBuffer);
1543 0 : return CE_Failure;
1544 : }
1545 50 : pDataBuf = padfTempBuffer;
1546 50 : eBufDataType = GDT_Float64;
1547 : }
1548 :
1549 91 : if (nTasks > 1)
1550 : {
1551 88 : std::vector<GDALPansharpenJob> asJobs;
1552 44 : asJobs.resize(nTasks);
1553 44 : GDALPansharpenJob *pasJobs = &(asJobs[0]);
1554 : {
1555 88 : std::vector<void *> ahJobData;
1556 44 : ahJobData.resize(nTasks);
1557 : #ifdef DEBUG_TIMING
1558 : struct timeval tv;
1559 : #endif
1560 220 : for (int i = 0; i < nTasks; i++)
1561 : {
1562 176 : const size_t iStartLine =
1563 176 : (static_cast<size_t>(i) * nYSize) / nTasks;
1564 176 : const size_t iNextStartLine =
1565 176 : (static_cast<size_t>(i + 1) * nYSize) / nTasks;
1566 176 : pasJobs[i].poPansharpenOperation = this;
1567 176 : pasJobs[i].eWorkDataType = eWorkDataType;
1568 176 : pasJobs[i].eBufDataType = eBufDataType;
1569 176 : pasJobs[i].pPanBuffer =
1570 176 : pPanBuffer + iStartLine * nXSize * nDataTypeSize;
1571 176 : pasJobs[i].pUpsampledSpectralBuffer =
1572 176 : pUpsampledSpectralBuffer +
1573 176 : iStartLine * nXSize * nDataTypeSize;
1574 176 : pasJobs[i].pDataBuf =
1575 176 : static_cast<GByte *>(pDataBuf) +
1576 352 : iStartLine * nXSize *
1577 176 : GDALGetDataTypeSizeBytes(eBufDataType);
1578 176 : pasJobs[i].nValues = (iNextStartLine - iStartLine) * nXSize;
1579 176 : pasJobs[i].nBandValues = static_cast<size_t>(nXSize) * nYSize;
1580 176 : pasJobs[i].nMaxValue = nMaxValue;
1581 : #ifdef DEBUG_TIMING
1582 : pasJobs[i].ptv = &tv;
1583 : #endif
1584 176 : ahJobData[i] = &(pasJobs[i]);
1585 : }
1586 : #ifdef DEBUG_TIMING
1587 : gettimeofday(&tv, nullptr);
1588 : #endif
1589 44 : poThreadPool->SubmitJobs(PansharpenJobThreadFunc, ahJobData);
1590 44 : poThreadPool->WaitCompletion();
1591 : }
1592 :
1593 44 : eErr = CE_None;
1594 220 : for (int i = 0; i < nTasks; i++)
1595 : {
1596 176 : if (pasJobs[i].eErr != CE_None)
1597 0 : eErr = CE_Failure;
1598 : }
1599 : }
1600 : else
1601 : {
1602 47 : eErr = PansharpenChunk(eWorkDataType, eBufDataType, pPanBuffer,
1603 : pUpsampledSpectralBuffer, pDataBuf,
1604 47 : static_cast<size_t>(nXSize) * nYSize,
1605 47 : static_cast<size_t>(nXSize) * nYSize, nMaxValue);
1606 : }
1607 :
1608 91 : if (padfTempBuffer)
1609 : {
1610 50 : GDALCopyWords64(padfTempBuffer, GDT_Float64, sizeof(double),
1611 : pDataBufOri, eBufDataTypeOri,
1612 : GDALGetDataTypeSizeBytes(eBufDataTypeOri),
1613 50 : static_cast<size_t>(nXSize) * nYSize *
1614 50 : psOptions->nOutPansharpenedBands);
1615 50 : VSIFree(padfTempBuffer);
1616 : }
1617 :
1618 91 : VSIFree(pUpsampledSpectralBuffer);
1619 91 : VSIFree(pPanBuffer);
1620 :
1621 91 : return eErr;
1622 : }
1623 :
1624 : /************************************************************************/
1625 : /* PansharpenResampleJobThreadFunc() */
1626 : /************************************************************************/
1627 :
1628 176 : void GDALPansharpenOperation::PansharpenResampleJobThreadFunc(void *pUserData)
1629 : {
1630 176 : GDALPansharpenResampleJob *psJob =
1631 : static_cast<GDALPansharpenResampleJob *>(pUserData);
1632 :
1633 : #ifdef DEBUG_TIMING
1634 : struct timeval tv;
1635 : gettimeofday(&tv, nullptr);
1636 : const GIntBig launch_time =
1637 : static_cast<GIntBig>(psJob->ptv->tv_sec) * 1000000 +
1638 : static_cast<GIntBig>(psJob->ptv->tv_usec);
1639 : const GIntBig start_job = static_cast<GIntBig>(tv.tv_sec) * 1000000 +
1640 : static_cast<GIntBig>(tv.tv_usec);
1641 : #endif
1642 :
1643 : GDALRasterIOExtraArg sExtraArg;
1644 176 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
1645 : // cppcheck-suppress redundantAssignment
1646 176 : sExtraArg.eResampleAlg = psJob->eResampleAlg;
1647 176 : sExtraArg.bFloatingPointWindowValidity = TRUE;
1648 176 : sExtraArg.dfXOff = psJob->dfXOff;
1649 176 : sExtraArg.dfYOff = psJob->dfYOff;
1650 176 : sExtraArg.dfXSize = psJob->dfXSize;
1651 176 : sExtraArg.dfYSize = psJob->dfYSize;
1652 :
1653 352 : std::vector<int> anBands;
1654 760 : for (int i = 0; i < psJob->nBandCount; ++i)
1655 584 : anBands.push_back(i + 1);
1656 : // This call to RasterIO() in a thread to poMEMDS shared between several
1657 : // threads is really risky, but works given the implementation details...
1658 : // Do not do this at home!
1659 176 : psJob->eErr = psJob->poMEMDS->RasterIO(
1660 : GF_Read, psJob->nXOff, psJob->nYOff, psJob->nXSize, psJob->nYSize,
1661 : psJob->pBuffer, psJob->nBufXSize, psJob->nBufYSize, psJob->eDT,
1662 176 : psJob->nBandCount, anBands.data(), 0, 0, psJob->nBandSpace, &sExtraArg);
1663 176 : if (CPLGetLastErrorType() == CE_Failure)
1664 0 : psJob->osLastErrorMsg = CPLGetLastErrorMsg();
1665 :
1666 : #ifdef DEBUG_TIMING
1667 : struct timeval tv_end;
1668 : gettimeofday(&tv_end, nullptr);
1669 : const GIntBig end = static_cast<GIntBig>(tv_end.tv_sec) * 1000000 +
1670 : static_cast<GIntBig>(tv_end.tv_usec);
1671 : if (start_job - launch_time > 500)
1672 : /*ok*/ printf("Resample: Delay before start=" CPL_FRMT_GIB
1673 : ", completion time=" CPL_FRMT_GIB "\n",
1674 : start_job - launch_time, end - start_job);
1675 : #endif
1676 176 : }
1677 :
1678 : /************************************************************************/
1679 : /* PansharpenJobThreadFunc() */
1680 : /************************************************************************/
1681 :
1682 176 : void GDALPansharpenOperation::PansharpenJobThreadFunc(void *pUserData)
1683 : {
1684 176 : GDALPansharpenJob *psJob = static_cast<GDALPansharpenJob *>(pUserData);
1685 :
1686 : #ifdef DEBUG_TIMING
1687 : struct timeval tv;
1688 : gettimeofday(&tv, nullptr);
1689 : const GIntBig launch_time =
1690 : static_cast<GIntBig>(psJob->ptv->tv_sec) * 1000000 +
1691 : static_cast<GIntBig>(psJob->ptv->tv_usec);
1692 : const GIntBig start_job = static_cast<GIntBig>(tv.tv_sec) * 1000000 +
1693 : static_cast<GIntBig>(tv.tv_usec);
1694 : #endif
1695 :
1696 : #if 0
1697 : for( int i = 0; i < 1000000; i++ )
1698 : acc += i * i;
1699 : psJob->eErr = CE_None;
1700 : #else
1701 176 : psJob->eErr = psJob->poPansharpenOperation->PansharpenChunk(
1702 : psJob->eWorkDataType, psJob->eBufDataType, psJob->pPanBuffer,
1703 : psJob->pUpsampledSpectralBuffer, psJob->pDataBuf, psJob->nValues,
1704 : psJob->nBandValues, psJob->nMaxValue);
1705 : #endif
1706 :
1707 : #ifdef DEBUG_TIMING
1708 : struct timeval tv_end;
1709 : gettimeofday(&tv_end, nullptr);
1710 : const GIntBig end = static_cast<GIntBig>(tv_end.tv_sec) * 1000000 +
1711 : static_cast<GIntBig>(tv_end.tv_usec);
1712 : if (start_job - launch_time > 500)
1713 : /*ok*/ printf("Pansharpen: Delay before start=" CPL_FRMT_GIB
1714 : ", completion time=" CPL_FRMT_GIB "\n",
1715 : start_job - launch_time, end - start_job);
1716 : #endif
1717 176 : }
1718 :
1719 : /************************************************************************/
1720 : /* PansharpenChunk() */
1721 : /************************************************************************/
1722 :
1723 223 : CPLErr GDALPansharpenOperation::PansharpenChunk(
1724 : GDALDataType eWorkDataType, GDALDataType eBufDataType,
1725 : const void *pPanBuffer, const void *pUpsampledSpectralBuffer,
1726 : void *pDataBuf, size_t nValues, size_t nBandValues, GUInt32 nMaxValue) const
1727 : {
1728 223 : CPLErr eErr = CE_None;
1729 :
1730 223 : switch (eWorkDataType)
1731 : {
1732 109 : case GDT_UInt8:
1733 109 : eErr = WeightedBrovey(
1734 : static_cast<const GByte *>(pPanBuffer),
1735 : static_cast<const GByte *>(pUpsampledSpectralBuffer), pDataBuf,
1736 : eBufDataType, nValues, nBandValues,
1737 : static_cast<GByte>(nMaxValue));
1738 109 : break;
1739 :
1740 108 : case GDT_UInt16:
1741 108 : eErr = WeightedBrovey(
1742 : static_cast<const GUInt16 *>(pPanBuffer),
1743 : static_cast<const GUInt16 *>(pUpsampledSpectralBuffer),
1744 : pDataBuf, eBufDataType, nValues, nBandValues,
1745 : static_cast<GUInt16>(nMaxValue));
1746 108 : break;
1747 :
1748 : #ifndef LIMIT_TYPES
1749 : case GDT_Int8:
1750 : eErr = WeightedBrovey(
1751 : static_cast<const GInt8 *>(pPanBuffer),
1752 : static_cast<const GInt8 *>(pUpsampledSpectralBuffer), pDataBuf,
1753 : eBufDataType, nValues, nBandValues);
1754 : break;
1755 :
1756 : case GDT_Int16:
1757 : eErr = WeightedBrovey(
1758 : static_cast<const GInt16 *>(pPanBuffer),
1759 : static_cast<const GInt16 *>(pUpsampledSpectralBuffer), pDataBuf,
1760 : eBufDataType, nValues, nBandValues);
1761 : break;
1762 :
1763 : case GDT_UInt32:
1764 : eErr = WeightedBrovey(
1765 : static_cast<const GUInt32 *>(pPanBuffer),
1766 : static_cast<const GUInt32 *>(pUpsampledSpectralBuffer),
1767 : pDataBuf, eBufDataType, nValues, nBandValues, nMaxValue);
1768 : break;
1769 :
1770 : case GDT_Int32:
1771 : eErr = WeightedBrovey(
1772 : static_cast<const GInt32 *>(pPanBuffer),
1773 : static_cast<const GInt32 *>(pUpsampledSpectralBuffer), pDataBuf,
1774 : eBufDataType, nValues, nBandValues);
1775 : break;
1776 :
1777 : case GDT_UInt64:
1778 : eErr = WeightedBrovey(
1779 : static_cast<const std::uint64_t *>(pPanBuffer),
1780 : static_cast<const std::uint64_t *>(pUpsampledSpectralBuffer),
1781 : pDataBuf, eBufDataType, nValues, nBandValues, nMaxValue);
1782 : break;
1783 :
1784 : case GDT_Int64:
1785 : eErr = WeightedBrovey(
1786 : static_cast<const std::int64_t *>(pPanBuffer),
1787 : static_cast<const std::int64_t *>(pUpsampledSpectralBuffer),
1788 : pDataBuf, eBufDataType, nValues, nBandValues);
1789 : break;
1790 :
1791 : case GDT_Float16:
1792 : eErr = WeightedBrovey(
1793 : static_cast<const GFloat16 *>(pPanBuffer),
1794 : static_cast<const GFloat16 *>(pUpsampledSpectralBuffer),
1795 : pDataBuf, eBufDataType, nValues, nBandValues);
1796 : break;
1797 :
1798 : case GDT_Float32:
1799 : eErr = WeightedBrovey(
1800 : static_cast<const float *>(pPanBuffer),
1801 : static_cast<const float *>(pUpsampledSpectralBuffer), pDataBuf,
1802 : eBufDataType, nValues, nBandValues);
1803 : break;
1804 : #endif
1805 6 : case GDT_Float64:
1806 6 : eErr = WeightedBrovey(
1807 : static_cast<const double *>(pPanBuffer),
1808 : static_cast<const double *>(pUpsampledSpectralBuffer), pDataBuf,
1809 : eBufDataType, nValues, nBandValues);
1810 6 : break;
1811 :
1812 0 : default:
1813 0 : CPLError(CE_Failure, CPLE_NotSupported,
1814 : "eWorkDataType not supported");
1815 0 : eErr = CE_Failure;
1816 0 : break;
1817 : }
1818 :
1819 223 : return eErr;
1820 : }
1821 :
1822 : /************************************************************************/
1823 : /* GetOptions() */
1824 : /************************************************************************/
1825 :
1826 : /** Return options.
1827 : * @return options.
1828 : */
1829 154 : GDALPansharpenOptions *GDALPansharpenOperation::GetOptions()
1830 : {
1831 154 : return psOptions;
1832 : }
1833 :
1834 : /************************************************************************/
1835 : /* GDALCreatePansharpenOperation() */
1836 : /************************************************************************/
1837 :
1838 : /** Instantiate a pansharpening operation.
1839 : *
1840 : * The passed options are validated.
1841 : *
1842 : * @param psOptions a pansharpening option structure allocated with
1843 : * GDALCreatePansharpenOptions(). It is duplicated by this function.
1844 : * @return a valid pansharpening operation handle, or NULL in case of failure.
1845 : *
1846 : */
1847 :
1848 : GDALPansharpenOperationH
1849 0 : GDALCreatePansharpenOperation(const GDALPansharpenOptions *psOptions)
1850 : {
1851 0 : GDALPansharpenOperation *psOperation = new GDALPansharpenOperation();
1852 0 : if (psOperation->Initialize(psOptions) == CE_None)
1853 0 : return reinterpret_cast<GDALPansharpenOperationH>(psOperation);
1854 0 : delete psOperation;
1855 0 : return nullptr;
1856 : }
1857 :
1858 : /************************************************************************/
1859 : /* GDALDestroyPansharpenOperation() */
1860 : /************************************************************************/
1861 :
1862 : /** Destroy a pansharpening operation.
1863 : *
1864 : * @param hOperation a valid pansharpening operation.
1865 : *
1866 : */
1867 :
1868 0 : void GDALDestroyPansharpenOperation(GDALPansharpenOperationH hOperation)
1869 : {
1870 0 : delete reinterpret_cast<GDALPansharpenOperation *>(hOperation);
1871 0 : }
1872 :
1873 : /************************************************************************/
1874 : /* GDALPansharpenProcessRegion() */
1875 : /************************************************************************/
1876 :
1877 : /** Executes a pansharpening operation on a rectangular region of the
1878 : * resulting dataset.
1879 : *
1880 : * The window is expressed with respect to the dimensions of the panchromatic
1881 : * band.
1882 : *
1883 : * Spectral bands are upsampled and merged with the panchromatic band according
1884 : * to the select algorithm and options.
1885 : *
1886 : * @param hOperation a valid pansharpening operation.
1887 : * @param nXOff pixel offset.
1888 : * @param nYOff pixel offset.
1889 : * @param nXSize width of the pansharpened region to compute.
1890 : * @param nYSize height of the pansharpened region to compute.
1891 : * @param pDataBuf output buffer. Must be nXSize * nYSize *
1892 : * GDALGetDataTypeSizeBytes(eBufDataType) *
1893 : * psOptions->nOutPansharpenedBands large.
1894 : * It begins with all values of the first output band, followed
1895 : * by values of the second output band, etc...
1896 : * @param eBufDataType data type of the output buffer
1897 : *
1898 : * @return CE_None in case of success, CE_Failure in case of failure.
1899 : *
1900 : */
1901 0 : CPLErr GDALPansharpenProcessRegion(GDALPansharpenOperationH hOperation,
1902 : int nXOff, int nYOff, int nXSize, int nYSize,
1903 : void *pDataBuf, GDALDataType eBufDataType)
1904 : {
1905 : return reinterpret_cast<GDALPansharpenOperation *>(hOperation)
1906 0 : ->ProcessRegion(nXOff, nYOff, nXSize, nYSize, pDataBuf, eBufDataType);
1907 : }
|
