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