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