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 111 : GDALPansharpenOptions *GDALCreatePansharpenOptions()
54 : {
55 : GDALPansharpenOptions *psOptions = static_cast<GDALPansharpenOptions *>(
56 111 : CPLCalloc(1, sizeof(GDALPansharpenOptions)));
57 111 : psOptions->ePansharpenAlg = GDAL_PSH_WEIGHTED_BROVEY;
58 111 : psOptions->eResampleAlg = GRIORA_Cubic;
59 111 : 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 112 : void GDALDestroyPansharpenOptions(GDALPansharpenOptions *psOptions)
75 : {
76 112 : if (psOptions == nullptr)
77 1 : return;
78 111 : CPLFree(psOptions->padfWeights);
79 111 : CPLFree(psOptions->pahInputSpectralBands);
80 111 : CPLFree(psOptions->panOutPansharpenedBands);
81 111 : 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 58 : GDALClonePansharpenOptions(const GDALPansharpenOptions *psOptions)
100 : {
101 58 : GDALPansharpenOptions *psNewOptions = GDALCreatePansharpenOptions();
102 58 : psNewOptions->ePansharpenAlg = psOptions->ePansharpenAlg;
103 58 : psNewOptions->eResampleAlg = psOptions->eResampleAlg;
104 58 : psNewOptions->nBitDepth = psOptions->nBitDepth;
105 58 : psNewOptions->nWeightCount = psOptions->nWeightCount;
106 58 : if (psOptions->padfWeights)
107 : {
108 58 : psNewOptions->padfWeights = static_cast<double *>(
109 58 : CPLMalloc(sizeof(double) * psOptions->nWeightCount));
110 58 : memcpy(psNewOptions->padfWeights, psOptions->padfWeights,
111 58 : sizeof(double) * psOptions->nWeightCount);
112 : }
113 58 : psNewOptions->hPanchroBand = psOptions->hPanchroBand;
114 58 : psNewOptions->nInputSpectralBands = psOptions->nInputSpectralBands;
115 58 : if (psOptions->pahInputSpectralBands)
116 : {
117 58 : const size_t nSize =
118 58 : sizeof(GDALRasterBandH) * psOptions->nInputSpectralBands;
119 58 : psNewOptions->pahInputSpectralBands =
120 58 : static_cast<GDALRasterBandH *>(CPLMalloc(nSize));
121 58 : memcpy(psNewOptions->pahInputSpectralBands,
122 58 : psOptions->pahInputSpectralBands, nSize);
123 : }
124 58 : psNewOptions->nOutPansharpenedBands = psOptions->nOutPansharpenedBands;
125 58 : if (psOptions->panOutPansharpenedBands)
126 : {
127 57 : psNewOptions->panOutPansharpenedBands = static_cast<int *>(
128 57 : CPLMalloc(sizeof(int) * psOptions->nOutPansharpenedBands));
129 57 : memcpy(psNewOptions->panOutPansharpenedBands,
130 57 : psOptions->panOutPansharpenedBands,
131 57 : sizeof(int) * psOptions->nOutPansharpenedBands);
132 : }
133 58 : psNewOptions->bHasNoData = psOptions->bHasNoData;
134 58 : psNewOptions->dfNoData = psOptions->dfNoData;
135 58 : psNewOptions->nThreads = psOptions->nThreads;
136 58 : 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 56 : GDALPansharpenOperation::~GDALPansharpenOperation()
157 : {
158 56 : GDALDestroyPansharpenOptions(psOptions);
159 63 : for (size_t i = 0; i < aVDS.size(); i++)
160 7 : delete aVDS[i];
161 56 : delete poThreadPool;
162 56 : }
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 56 : GDALPansharpenOperation::Initialize(const GDALPansharpenOptions *psOptionsIn)
176 : {
177 56 : if (psOptionsIn->hPanchroBand == nullptr)
178 : {
179 0 : CPLError(CE_Failure, CPLE_AppDefined, "hPanchroBand not set");
180 0 : return CE_Failure;
181 : }
182 56 : 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 56 : if (psOptionsIn->padfWeights == nullptr ||
189 56 : 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 56 : auto poPanchroBand = GDALRasterBand::FromHandle(psOptionsIn->hPanchroBand);
198 56 : auto poPanchroDS = poPanchroBand->GetDataset();
199 56 : 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 56 : 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 56 : 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 56 : GDALRasterBandH hRefBand = psOptionsIn->pahInputSpectralBands[0];
222 56 : auto poRefBand = GDALRasterBand::FromHandle(hRefBand);
223 56 : auto poRefBandDS = poRefBand->GetDataset();
224 56 : 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 56 : 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 56 : 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 56 : 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 112 : m_adfPanToMSGT[1] =
258 56 : adfInvMSGT[1] * adfPanchroGT[1] + adfInvMSGT[2] * adfPanchroGT[4];
259 112 : m_adfPanToMSGT[2] =
260 56 : adfInvMSGT[1] * adfPanchroGT[2] + adfInvMSGT[2] * adfPanchroGT[5];
261 56 : m_adfPanToMSGT[0] = adfInvMSGT[1] * adfPanchroGT[0] +
262 56 : adfInvMSGT[2] * adfPanchroGT[3] + adfInvMSGT[0];
263 112 : m_adfPanToMSGT[4] =
264 56 : adfInvMSGT[4] * adfPanchroGT[1] + adfInvMSGT[5] * adfPanchroGT[4];
265 112 : m_adfPanToMSGT[5] =
266 56 : adfInvMSGT[4] * adfPanchroGT[2] + adfInvMSGT[5] * adfPanchroGT[5];
267 56 : m_adfPanToMSGT[3] = adfInvMSGT[4] * adfPanchroGT[0] +
268 56 : 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 112 : if (std::fabs(m_adfPanToMSGT[2]) > 1e-10 ||
275 56 : 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 56 : bool bSameDataset = psOptionsIn->nInputSpectralBands > 1;
284 56 : if (bSameDataset)
285 41 : anInputBands.push_back(GDALGetBandNumber(hRefBand));
286 133 : for (int i = 1; i < psOptionsIn->nInputSpectralBands; i++)
287 : {
288 78 : GDALRasterBandH hBand = psOptionsIn->pahInputSpectralBands[i];
289 155 : if (GDALGetRasterBandXSize(hBand) != GDALGetRasterBandXSize(hRefBand) ||
290 77 : 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 77 : auto poBand = GDALRasterBand::FromHandle(hBand);
300 77 : auto poBandDS = poBand->GetDataset();
301 77 : 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 77 : 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 77 : 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 77 : 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 77 : if (bSameDataset)
336 : {
337 75 : if (GDALGetBandDataset(hBand) != GDALGetBandDataset(hRefBand))
338 : {
339 4 : anInputBands.resize(0);
340 4 : bSameDataset = false;
341 : }
342 : else
343 : {
344 71 : anInputBands.push_back(GDALGetBandNumber(hBand));
345 : }
346 : }
347 : }
348 55 : if (psOptionsIn->nOutPansharpenedBands == 0)
349 : {
350 1 : CPLError(CE_Warning, CPLE_AppDefined,
351 : "No output pansharpened band defined");
352 : }
353 181 : for (int i = 0; i < psOptionsIn->nOutPansharpenedBands; i++)
354 : {
355 126 : if (psOptionsIn->panOutPansharpenedBands[i] < 0 ||
356 126 : psOptionsIn->panOutPansharpenedBands[i] >=
357 126 : 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 55 : GDALDataType eWorkDataType = poPanchroBand->GetRasterDataType();
367 55 : 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 55 : psOptions = GDALClonePansharpenOptions(psOptionsIn);
382 55 : if (psOptions->nBitDepth == GDALGetDataTypeSize(eWorkDataType))
383 6 : psOptions->nBitDepth = 0;
384 55 : 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 187 : for (int i = 0; i < psOptions->nInputSpectralBands; i++)
396 : {
397 132 : if (psOptions->padfWeights[i] < 0.0)
398 : {
399 0 : bPositiveWeights = FALSE;
400 0 : break;
401 : }
402 : }
403 :
404 187 : for (int i = 0; i < psOptions->nInputSpectralBands; i++)
405 : {
406 264 : aMSBands.push_back(
407 132 : GDALRasterBand::FromHandle(psOptions->pahInputSpectralBands[i]));
408 : }
409 :
410 55 : 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 55 : int nThreads = psOptions->nThreads;
468 55 : if (nThreads == -1)
469 9 : 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 55 : if (nThreads > 1)
483 : {
484 9 : CPLDebug("PANSHARPEN", "Using %d threads", nThreads);
485 18 : poThreadPool = new (std::nothrow) CPLWorkerThreadPool();
486 : // coverity[tainted_data]
487 18 : if (poThreadPool == nullptr ||
488 9 : !poThreadPool->Setup(nThreads, nullptr, nullptr))
489 : {
490 0 : delete poThreadPool;
491 0 : poThreadPool = nullptr;
492 : }
493 : }
494 :
495 55 : GDALRIOResampleAlg eResampleAlg = psOptions->eResampleAlg;
496 55 : if (eResampleAlg != GRIORA_NearestNeighbour)
497 : {
498 55 : const char *pszResampling =
499 110 : (eResampleAlg == GRIORA_Bilinear) ? "BILINEAR"
500 55 : : (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 55 : GDALGetResampleFunction(pszResampling, &nKernelRadius);
510 : }
511 :
512 55 : 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 19959220 : static inline double ComputeFactor(T panValue, double dfPseudoPanchro)
582 : {
583 19959220 : if (dfPseudoPanchro == 0.0)
584 2778 : return 0.0;
585 :
586 19956440 : return panValue / dfPseudoPanchro;
587 : }
588 :
589 : /************************************************************************/
590 : /* ClampAndRound() */
591 : /************************************************************************/
592 :
593 136 : template <class T> static inline T ClampAndRound(double dfVal, T nMaxValue)
594 : {
595 136 : if (dfVal > nMaxValue)
596 50 : return nMaxValue;
597 : else
598 86 : return static_cast<T>(dfVal + 0.5);
599 : }
600 :
601 : /************************************************************************/
602 : /* WeightedBrovey() */
603 : /************************************************************************/
604 :
605 : template <class WorkDataType, class OutDataType, int bHasBitDepth>
606 129 : 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 129 : if (psOptions->bHasNoData)
612 : {
613 4 : WeightedBroveyWithNoData<WorkDataType, OutDataType>(
614 : pPanBuffer, pUpsampledSpectralBuffer, pDataBuf, nValues,
615 : nBandValues, nMaxValue);
616 4 : return;
617 : }
618 :
619 19809402 : for (size_t j = 0; j < nValues; j++)
620 : {
621 19911540 : double dfFactor = 0.0;
622 : // if( pPanBuffer[j] == 0 )
623 : // dfFactor = 1.0;
624 : // else
625 : {
626 19911540 : double dfPseudoPanchro = 0.0;
627 85729500 : for (int i = 0; i < psOptions->nInputSpectralBands; i++)
628 65818000 : dfPseudoPanchro +=
629 65818000 : psOptions->padfWeights[i] *
630 65818000 : pUpsampledSpectralBuffer[i * nBandValues + j];
631 19911540 : dfFactor = ComputeFactor(pPanBuffer[j], dfPseudoPanchro);
632 : }
633 :
634 77761400 : for (int i = 0; i < psOptions->nOutPansharpenedBands; i++)
635 : {
636 57952200 : WorkDataType nRawValue =
637 57952200 : pUpsampledSpectralBuffer[psOptions->panOutPansharpenedBands[i] *
638 57952200 : nBandValues +
639 : j];
640 : WorkDataType nPansharpenedValue;
641 57952200 : GDALCopyWord(nRawValue * dfFactor, nPansharpenedValue);
642 : if constexpr (bHasBitDepth)
643 : {
644 0 : if (nPansharpenedValue > nMaxValue)
645 0 : nPansharpenedValue = nMaxValue;
646 : }
647 54579900 : 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 41 : 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 41 : const XMMReg4Double w0 =
667 41 : XMMReg4Double::Load1ValHighAndLow(psOptions->padfWeights + 0);
668 41 : const XMMReg4Double w1 =
669 41 : XMMReg4Double::Load1ValHighAndLow(psOptions->padfWeights + 1);
670 41 : const XMMReg4Double w2 =
671 41 : 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 41 : const XMMReg4Double zero = XMMReg4Double::Zero();
678 41 : double dfMaxValue = nMaxValue;
679 41 : const XMMReg4Double maxValue =
680 : XMMReg4Double::Load1ValHighAndLow(&dfMaxValue);
681 :
682 1530 : size_t j = 0; // Used after for.
683 1317533 : for (; j + 3 < nValues; j += 4)
684 : {
685 1317492 : XMMReg4Double pseudoPanchro = zero;
686 :
687 1320560 : XMMReg4Double val0 = XMMReg4Double::Load4Val(pUpsampledSpectralBuffer +
688 767133 : 0 * nBandValues + j);
689 1317653 : XMMReg4Double val1 = XMMReg4Double::Load4Val(pUpsampledSpectralBuffer +
690 1317653 : 1 * nBandValues + j);
691 1317790 : XMMReg4Double val2 = XMMReg4Double::Load4Val(pUpsampledSpectralBuffer +
692 1317790 : 2 * nBandValues + j);
693 : XMMReg4Double val3;
694 : if constexpr (NINPUT == 4 || NOUTPUT == 4)
695 : {
696 503476 : val3 = XMMReg4Double::Load4Val(pUpsampledSpectralBuffer +
697 511800 : 3 * nBandValues + j);
698 : }
699 :
700 1317206 : pseudoPanchro += w0 * val0;
701 1314749 : pseudoPanchro += w1 * val1;
702 1306639 : pseudoPanchro += w2 * val2;
703 : if constexpr (NINPUT == 4)
704 510680 : pseudoPanchro += w3 * val3;
705 :
706 : /* Little trick to avoid use of ternary operator due to one of the
707 : * branch being zero */
708 1319402 : XMMReg4Double factor = XMMReg4Double::And(
709 : XMMReg4Double::NotEquals(pseudoPanchro, zero),
710 763973 : XMMReg4Double::Load4Val(pPanBuffer + j) / pseudoPanchro);
711 :
712 1317346 : val0 = XMMReg4Double::Min(val0 * factor, maxValue);
713 1315967 : val1 = XMMReg4Double::Min(val1 * factor, maxValue);
714 1318707 : val2 = XMMReg4Double::Min(val2 * factor, maxValue);
715 : if constexpr (NOUTPUT == 4)
716 : {
717 256364 : val3 = XMMReg4Double::Min(val3 * factor, maxValue);
718 : }
719 1319768 : val0.Store4Val(pDataBuf + 0 * nBandValues + j);
720 1314662 : val1.Store4Val(pDataBuf + 1 * nBandValues + j);
721 1315983 : val2.Store4Val(pDataBuf + 2 * nBandValues + j);
722 : if constexpr (NOUTPUT == 4)
723 : {
724 254703 : 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 47 : void GDALPansharpenOperation::WeightedBroveyPositiveWeights(
792 : const T *pPanBuffer, const T *pUpsampledSpectralBuffer, T *pDataBuf,
793 : size_t nValues, size_t nBandValues, T nMaxValue) const
794 : {
795 47 : if (psOptions->bHasNoData)
796 : {
797 0 : WeightedBroveyWithNoData<T, T>(pPanBuffer, pUpsampledSpectralBuffer,
798 : pDataBuf, nValues, nBandValues,
799 : nMaxValue);
800 0 : return;
801 : }
802 :
803 47 : if (nMaxValue == 0)
804 41 : nMaxValue = cpl::NumericLimits<T>::max();
805 : size_t j;
806 47 : 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 14 : else if (psOptions->nInputSpectralBands == 4 &&
817 8 : psOptions->nOutPansharpenedBands == 4 &&
818 4 : 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 10 : else if (psOptions->nInputSpectralBands == 4 &&
828 4 : psOptions->nOutPansharpenedBands == 3 &&
829 4 : psOptions->panOutPansharpenedBands[0] == 0 &&
830 4 : psOptions->panOutPansharpenedBands[1] == 1 &&
831 4 : psOptions->panOutPansharpenedBands[2] == 2)
832 : {
833 4 : j = WeightedBroveyPositiveWeightsInternal<T, 4, 3>(
834 : pPanBuffer, pUpsampledSpectralBuffer, pDataBuf, nValues,
835 : nBandValues, nMaxValue);
836 : }
837 : else
838 : {
839 54 : for (j = 0; j + 1 < nValues; j += 2)
840 : {
841 48 : double dfFactor = 0.0;
842 48 : double dfFactor2 = 0.0;
843 48 : double dfPseudoPanchro = 0.0;
844 48 : double dfPseudoPanchro2 = 0.0;
845 96 : for (int i = 0; i < psOptions->nInputSpectralBands; i++)
846 : {
847 48 : dfPseudoPanchro +=
848 48 : psOptions->padfWeights[i] *
849 48 : pUpsampledSpectralBuffer[i * nBandValues + j];
850 48 : dfPseudoPanchro2 +=
851 48 : psOptions->padfWeights[i] *
852 48 : pUpsampledSpectralBuffer[i * nBandValues + j + 1];
853 : }
854 :
855 48 : dfFactor = ComputeFactor(pPanBuffer[j], dfPseudoPanchro);
856 48 : dfFactor2 = ComputeFactor(pPanBuffer[j + 1], dfPseudoPanchro2);
857 :
858 96 : for (int i = 0; i < psOptions->nOutPansharpenedBands; i++)
859 : {
860 48 : const T nRawValue = pUpsampledSpectralBuffer
861 48 : [psOptions->panOutPansharpenedBands[i] * nBandValues + j];
862 48 : const double dfTmp = nRawValue * dfFactor;
863 48 : pDataBuf[i * nBandValues + j] = ClampAndRound(dfTmp, nMaxValue);
864 :
865 48 : const T nRawValue2 = pUpsampledSpectralBuffer
866 48 : [psOptions->panOutPansharpenedBands[i] * nBandValues + j +
867 : 1];
868 48 : const double dfTmp2 = nRawValue2 * dfFactor2;
869 48 : pDataBuf[i * nBandValues + j + 1] =
870 48 : ClampAndRound(dfTmp2, nMaxValue);
871 : }
872 : }
873 : }
874 60 : 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 47 : void GDALPansharpenOperation::WeightedBroveyGByteOrUInt16(
915 : const T *pPanBuffer, const T *pUpsampledSpectralBuffer, T *pDataBuf,
916 : size_t nValues, size_t nBandValues, T nMaxValue) const
917 : {
918 47 : if (bPositiveWeights)
919 : {
920 47 : 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 47 : }
935 :
936 : template <>
937 32 : 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 32 : WeightedBroveyGByteOrUInt16(pPanBuffer, pUpsampledSpectralBuffer, pDataBuf,
942 : nValues, nBandValues, nMaxValue);
943 32 : }
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 173 : 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 173 : switch (eBufDataType)
963 : {
964 33 : case GDT_Byte:
965 33 : WeightedBrovey(pPanBuffer, pUpsampledSpectralBuffer,
966 : static_cast<GByte *>(pDataBuf), nValues, nBandValues,
967 : nMaxValue);
968 33 : 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 124 : case GDT_Float64:
1027 124 : WeightedBrovey(pPanBuffer, pUpsampledSpectralBuffer,
1028 : static_cast<double *>(pDataBuf), nValues,
1029 : nBandValues, nMaxValue);
1030 124 : break;
1031 :
1032 0 : default:
1033 0 : CPLError(CE_Failure, CPLE_NotSupported,
1034 : "eBufDataType not supported");
1035 0 : return CE_Failure;
1036 : break;
1037 : }
1038 :
1039 173 : 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 80 : CPLErr GDALPansharpenOperation::ProcessRegion(int nXOff, int nYOff, int nXSize,
1172 : int nYSize, void *pDataBuf,
1173 : GDALDataType eBufDataType)
1174 : {
1175 80 : if (psOptions == nullptr)
1176 0 : return CE_Failure;
1177 :
1178 : // TODO: Avoid allocating buffers each time.
1179 : GDALRasterBand *poPanchroBand =
1180 80 : GDALRasterBand::FromHandle(psOptions->hPanchroBand);
1181 80 : GDALDataType eWorkDataType = poPanchroBand->GetRasterDataType();
1182 : #ifdef LIMIT_TYPES
1183 80 : if (eWorkDataType != GDT_Byte && eWorkDataType != GDT_UInt16)
1184 6 : eWorkDataType = GDT_Float64;
1185 : #endif
1186 80 : const int nDataTypeSize = GDALGetDataTypeSizeBytes(eWorkDataType);
1187 80 : 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 80 : VSI_MALLOC3_VERBOSE(nXSize, nYSize, nDataTypeSize));
1192 80 : if (pUpsampledSpectralBuffer == nullptr || pPanBuffer == nullptr)
1193 : {
1194 0 : VSIFree(pUpsampledSpectralBuffer);
1195 0 : VSIFree(pPanBuffer);
1196 0 : return CE_Failure;
1197 : }
1198 :
1199 80 : CPLErr eErr = poPanchroBand->RasterIO(GF_Read, nXOff, nYOff, nXSize, nYSize,
1200 : pPanBuffer, nXSize, nYSize,
1201 : eWorkDataType, 0, 0, nullptr);
1202 80 : if (eErr != CE_None)
1203 : {
1204 0 : VSIFree(pUpsampledSpectralBuffer);
1205 0 : VSIFree(pPanBuffer);
1206 0 : return CE_Failure;
1207 : }
1208 :
1209 80 : int nTasks = 0;
1210 80 : if (poThreadPool)
1211 : {
1212 33 : nTasks = poThreadPool->GetThreadCount();
1213 33 : if (nTasks > nYSize)
1214 0 : nTasks = nYSize;
1215 : }
1216 :
1217 : GDALRasterIOExtraArg sExtraArg;
1218 80 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
1219 80 : const GDALRIOResampleAlg eResampleAlg = psOptions->eResampleAlg;
1220 : // cppcheck-suppress redundantAssignment
1221 80 : sExtraArg.eResampleAlg = eResampleAlg;
1222 80 : sExtraArg.bFloatingPointWindowValidity = TRUE;
1223 80 : sExtraArg.dfXOff = m_adfPanToMSGT[0] + nXOff * m_adfPanToMSGT[1];
1224 80 : sExtraArg.dfYOff = m_adfPanToMSGT[3] + nYOff * m_adfPanToMSGT[5];
1225 80 : sExtraArg.dfXSize = nXSize * m_adfPanToMSGT[1];
1226 80 : sExtraArg.dfYSize = nYSize * m_adfPanToMSGT[5];
1227 80 : if (sExtraArg.dfXOff + sExtraArg.dfXSize > aMSBands[0]->GetXSize())
1228 0 : sExtraArg.dfXOff = aMSBands[0]->GetXSize() - sExtraArg.dfXSize;
1229 80 : if (sExtraArg.dfYOff + sExtraArg.dfYSize > aMSBands[0]->GetYSize())
1230 0 : sExtraArg.dfYOff = aMSBands[0]->GetYSize() - sExtraArg.dfYSize;
1231 80 : int nSpectralXOff = static_cast<int>(sExtraArg.dfXOff);
1232 80 : int nSpectralYOff = static_cast<int>(sExtraArg.dfYOff);
1233 80 : int nSpectralXSize = static_cast<int>(0.49999 + sExtraArg.dfXSize);
1234 80 : int nSpectralYSize = static_cast<int>(0.49999 + sExtraArg.dfYSize);
1235 80 : if (nSpectralXSize == 0)
1236 10 : nSpectralXSize = 1;
1237 80 : if (nSpectralYSize == 0)
1238 10 : nSpectralYSize = 1;
1239 :
1240 : // When upsampling, extract the multispectral data at
1241 : // full resolution in a temp buffer, and then do the upsampling.
1242 80 : if (nSpectralXSize < nXSize && nSpectralYSize < nYSize &&
1243 70 : eResampleAlg != GRIORA_NearestNeighbour && nYSize > 1)
1244 : {
1245 : // Take some margin to take into account the radius of the
1246 : // resampling kernel.
1247 70 : int nXOffExtract = nSpectralXOff - nKernelRadius;
1248 70 : int nYOffExtract = nSpectralYOff - nKernelRadius;
1249 70 : int nXSizeExtract = nSpectralXSize + 1 + 2 * nKernelRadius;
1250 70 : int nYSizeExtract = nSpectralYSize + 1 + 2 * nKernelRadius;
1251 70 : if (nXOffExtract < 0)
1252 : {
1253 65 : nXSizeExtract += nXOffExtract;
1254 65 : nXOffExtract = 0;
1255 : }
1256 70 : if (nYOffExtract < 0)
1257 : {
1258 58 : nYSizeExtract += nYOffExtract;
1259 58 : nYOffExtract = 0;
1260 : }
1261 70 : if (nXOffExtract + nXSizeExtract > aMSBands[0]->GetXSize())
1262 65 : nXSizeExtract = aMSBands[0]->GetXSize() - nXOffExtract;
1263 70 : if (nYOffExtract + nYSizeExtract > aMSBands[0]->GetYSize())
1264 58 : nYSizeExtract = aMSBands[0]->GetYSize() - nYOffExtract;
1265 :
1266 70 : GByte *pSpectralBuffer = static_cast<GByte *>(VSI_MALLOC3_VERBOSE(
1267 : nXSizeExtract, nYSizeExtract,
1268 : cpl::fits_on<int>(psOptions->nInputSpectralBands * nDataTypeSize)));
1269 70 : if (pSpectralBuffer == nullptr)
1270 : {
1271 0 : VSIFree(pUpsampledSpectralBuffer);
1272 0 : VSIFree(pPanBuffer);
1273 0 : return CE_Failure;
1274 : }
1275 :
1276 70 : if (!anInputBands.empty())
1277 : {
1278 : // Use dataset RasterIO when possible.
1279 124 : eErr = aMSBands[0]->GetDataset()->RasterIO(
1280 : GF_Read, nXOffExtract, nYOffExtract, nXSizeExtract,
1281 : nYSizeExtract, pSpectralBuffer, nXSizeExtract, nYSizeExtract,
1282 62 : eWorkDataType, static_cast<int>(anInputBands.size()),
1283 62 : &anInputBands[0], 0, 0, 0, nullptr);
1284 : }
1285 : else
1286 : {
1287 8 : for (int i = 0;
1288 20 : eErr == CE_None && i < psOptions->nInputSpectralBands; i++)
1289 : {
1290 24 : eErr = aMSBands[i]->RasterIO(
1291 : GF_Read, nXOffExtract, nYOffExtract, nXSizeExtract,
1292 : nYSizeExtract,
1293 12 : pSpectralBuffer + static_cast<size_t>(i) * nXSizeExtract *
1294 12 : nYSizeExtract * nDataTypeSize,
1295 : nXSizeExtract, nYSizeExtract, eWorkDataType, 0, 0, nullptr);
1296 : }
1297 : }
1298 70 : if (eErr != CE_None)
1299 : {
1300 0 : VSIFree(pSpectralBuffer);
1301 0 : VSIFree(pUpsampledSpectralBuffer);
1302 0 : VSIFree(pPanBuffer);
1303 0 : return CE_Failure;
1304 : }
1305 :
1306 : // Create a MEM dataset that wraps the input buffer.
1307 70 : auto poMEMDS = MEMDataset::Create("", nXSizeExtract, nYSizeExtract, 0,
1308 : eWorkDataType, nullptr);
1309 :
1310 282 : for (int i = 0; i < psOptions->nInputSpectralBands; i++)
1311 : {
1312 212 : GByte *pabyBuffer =
1313 212 : pSpectralBuffer + static_cast<size_t>(i) * nDataTypeSize *
1314 212 : nXSizeExtract * nYSizeExtract;
1315 212 : GDALRasterBandH hMEMBand = MEMCreateRasterBandEx(
1316 : poMEMDS, i + 1, pabyBuffer, eWorkDataType, 0, 0, false);
1317 212 : poMEMDS->AddMEMBand(hMEMBand);
1318 :
1319 : const char *pszNBITS =
1320 212 : aMSBands[i]->GetMetadataItem("NBITS", "IMAGE_STRUCTURE");
1321 212 : if (pszNBITS)
1322 4 : poMEMDS->GetRasterBand(i + 1)->SetMetadataItem(
1323 4 : "NBITS", pszNBITS, "IMAGE_STRUCTURE");
1324 :
1325 212 : if (psOptions->bHasNoData)
1326 10 : poMEMDS->GetRasterBand(i + 1)->SetNoDataValue(
1327 10 : psOptions->dfNoData);
1328 : }
1329 :
1330 70 : if (nTasks <= 1)
1331 : {
1332 37 : nSpectralXOff -= nXOffExtract;
1333 37 : nSpectralYOff -= nYOffExtract;
1334 37 : sExtraArg.dfXOff -= nXOffExtract;
1335 37 : sExtraArg.dfYOff -= nYOffExtract;
1336 37 : CPL_IGNORE_RET_VAL(poMEMDS->RasterIO(
1337 : GF_Read, nSpectralXOff, nSpectralYOff, nSpectralXSize,
1338 : nSpectralYSize, pUpsampledSpectralBuffer, nXSize, nYSize,
1339 37 : eWorkDataType, psOptions->nInputSpectralBands, nullptr, 0, 0, 0,
1340 : &sExtraArg));
1341 : }
1342 : else
1343 : {
1344 : // We are abusing the contract of the GDAL API by using the
1345 : // MEMDataset from several threads. In this case, this is safe. In
1346 : // case, that would no longer be the case we could create as many
1347 : // MEMDataset as threads pointing to the same buffer.
1348 :
1349 : // To avoid races in threads, we query now the mask flags,
1350 : // so that implicit mask bands are created now.
1351 150 : for (int i = 0; i < poMEMDS->GetRasterCount(); i++)
1352 : {
1353 117 : poMEMDS->GetRasterBand(i + 1)->GetMaskFlags();
1354 : }
1355 :
1356 66 : std::vector<GDALPansharpenResampleJob> asJobs;
1357 33 : asJobs.resize(nTasks);
1358 33 : GDALPansharpenResampleJob *pasJobs = &(asJobs[0]);
1359 : {
1360 66 : std::vector<void *> ahJobData;
1361 33 : ahJobData.resize(nTasks);
1362 :
1363 : #ifdef DEBUG_TIMING
1364 : struct timeval tv;
1365 : #endif
1366 165 : for (int i = 0; i < nTasks; i++)
1367 : {
1368 132 : const size_t iStartLine =
1369 132 : (static_cast<size_t>(i) * nYSize) / nTasks;
1370 132 : const size_t iNextStartLine =
1371 132 : (static_cast<size_t>(i + 1) * nYSize) / nTasks;
1372 132 : pasJobs[i].poMEMDS = poMEMDS;
1373 132 : pasJobs[i].eResampleAlg = eResampleAlg;
1374 132 : pasJobs[i].dfXOff = sExtraArg.dfXOff - nXOffExtract;
1375 132 : pasJobs[i].dfYOff =
1376 132 : m_adfPanToMSGT[3] +
1377 132 : (nYOff + iStartLine) * m_adfPanToMSGT[5] - nYOffExtract;
1378 132 : pasJobs[i].dfXSize = sExtraArg.dfXSize;
1379 132 : pasJobs[i].dfYSize =
1380 132 : (iNextStartLine - iStartLine) * m_adfPanToMSGT[5];
1381 264 : if (pasJobs[i].dfXOff + pasJobs[i].dfXSize >
1382 132 : aMSBands[0]->GetXSize())
1383 : {
1384 0 : pasJobs[i].dfXOff =
1385 0 : aMSBands[0]->GetXSize() - pasJobs[i].dfXSize;
1386 : }
1387 264 : if (pasJobs[i].dfYOff + pasJobs[i].dfYSize >
1388 132 : aMSBands[0]->GetYSize())
1389 : {
1390 0 : pasJobs[i].dfYOff =
1391 0 : aMSBands[0]->GetYSize() - pasJobs[i].dfYSize;
1392 : }
1393 132 : pasJobs[i].nXOff = static_cast<int>(pasJobs[i].dfXOff);
1394 132 : pasJobs[i].nYOff = static_cast<int>(pasJobs[i].dfYOff);
1395 132 : pasJobs[i].nXSize =
1396 132 : static_cast<int>(0.4999 + pasJobs[i].dfXSize);
1397 132 : pasJobs[i].nYSize =
1398 132 : static_cast<int>(0.4999 + pasJobs[i].dfYSize);
1399 132 : if (pasJobs[i].nXSize == 0)
1400 0 : pasJobs[i].nXSize = 1;
1401 132 : if (pasJobs[i].nYSize == 0)
1402 0 : pasJobs[i].nYSize = 1;
1403 132 : pasJobs[i].pBuffer = pUpsampledSpectralBuffer +
1404 132 : static_cast<size_t>(iStartLine) *
1405 132 : nXSize * nDataTypeSize;
1406 132 : pasJobs[i].eDT = eWorkDataType;
1407 132 : pasJobs[i].nBufXSize = nXSize;
1408 132 : pasJobs[i].nBufYSize =
1409 132 : static_cast<int>(iNextStartLine - iStartLine);
1410 132 : pasJobs[i].nBandCount = psOptions->nInputSpectralBands;
1411 132 : pasJobs[i].nBandSpace =
1412 132 : static_cast<GSpacing>(nXSize) * nYSize * nDataTypeSize;
1413 : #ifdef DEBUG_TIMING
1414 : pasJobs[i].ptv = &tv;
1415 : #endif
1416 132 : ahJobData[i] = &(pasJobs[i]);
1417 : }
1418 : #ifdef DEBUG_TIMING
1419 : gettimeofday(&tv, nullptr);
1420 : #endif
1421 33 : poThreadPool->SubmitJobs(PansharpenResampleJobThreadFunc,
1422 : ahJobData);
1423 33 : poThreadPool->WaitCompletion();
1424 : }
1425 : }
1426 :
1427 70 : GDALClose(poMEMDS);
1428 :
1429 70 : VSIFree(pSpectralBuffer);
1430 : }
1431 : else
1432 : {
1433 10 : if (!anInputBands.empty())
1434 : {
1435 : // Use dataset RasterIO when possible.
1436 20 : eErr = aMSBands[0]->GetDataset()->RasterIO(
1437 : GF_Read, nSpectralXOff, nSpectralYOff, nSpectralXSize,
1438 : nSpectralYSize, pUpsampledSpectralBuffer, nXSize, nYSize,
1439 10 : eWorkDataType, static_cast<int>(anInputBands.size()),
1440 10 : &anInputBands[0], 0, 0, 0, &sExtraArg);
1441 : }
1442 : else
1443 : {
1444 0 : for (int i = 0;
1445 0 : eErr == CE_None && i < psOptions->nInputSpectralBands; i++)
1446 : {
1447 0 : eErr = aMSBands[i]->RasterIO(
1448 : GF_Read, nSpectralXOff, nSpectralYOff, nSpectralXSize,
1449 : nSpectralYSize,
1450 0 : pUpsampledSpectralBuffer + static_cast<size_t>(i) * nXSize *
1451 0 : nYSize * nDataTypeSize,
1452 : nXSize, nYSize, eWorkDataType, 0, 0, &sExtraArg);
1453 : }
1454 : }
1455 10 : if (eErr != CE_None)
1456 : {
1457 0 : VSIFree(pUpsampledSpectralBuffer);
1458 0 : VSIFree(pPanBuffer);
1459 0 : return CE_Failure;
1460 : }
1461 : }
1462 :
1463 : // In case NBITS was not set on the spectral bands, clamp the values
1464 : // if overshoot might have occurred.
1465 80 : int nBitDepth = psOptions->nBitDepth;
1466 80 : if (nBitDepth &&
1467 0 : (eResampleAlg == GRIORA_Cubic || eResampleAlg == GRIORA_CubicSpline ||
1468 : eResampleAlg == GRIORA_Lanczos))
1469 : {
1470 12 : for (int i = 0; i < psOptions->nInputSpectralBands; i++)
1471 : {
1472 6 : GDALRasterBand *poBand = aMSBands[i];
1473 6 : int nBandBitDepth = 0;
1474 : const char *pszNBITS =
1475 6 : poBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE");
1476 6 : if (pszNBITS)
1477 4 : nBandBitDepth = atoi(pszNBITS);
1478 6 : if (nBandBitDepth < nBitDepth)
1479 : {
1480 2 : if (eWorkDataType == GDT_Byte && nBitDepth >= 0 &&
1481 : nBitDepth <= 8)
1482 : {
1483 1 : ClampValues(
1484 : reinterpret_cast<GByte *>(pUpsampledSpectralBuffer) +
1485 1 : static_cast<size_t>(i) * nXSize * nYSize,
1486 1 : static_cast<size_t>(nXSize) * nYSize,
1487 1 : static_cast<GByte>((1 << nBitDepth) - 1));
1488 : }
1489 1 : else if (eWorkDataType == GDT_UInt16 && nBitDepth >= 0 &&
1490 : nBitDepth <= 16)
1491 : {
1492 1 : ClampValues(
1493 1 : reinterpret_cast<GUInt16 *>(pUpsampledSpectralBuffer) +
1494 1 : static_cast<size_t>(i) * nXSize * nYSize,
1495 1 : static_cast<size_t>(nXSize) * nYSize,
1496 1 : static_cast<GUInt16>((1 << nBitDepth) - 1));
1497 : }
1498 : #ifndef LIMIT_TYPES
1499 : else if (eWorkDataType == GDT_UInt32)
1500 : {
1501 : ClampValues(reinterpret_cast<GUInt32*>(pUpsampledSpectralBuffer) +
1502 : static_cast<size_t>(i) * nXSize * nYSize,
1503 : static_cast<size_t>(nXSize)*nYSize,
1504 : (static_cast<GUInt32>((1 << nBitDepth)-1));
1505 : }
1506 : #endif
1507 : }
1508 : }
1509 : }
1510 :
1511 80 : const GUInt32 nMaxValue = (nBitDepth >= 0 && nBitDepth <= 31)
1512 160 : ? (1U << nBitDepth) - 1
1513 : : UINT32_MAX;
1514 :
1515 80 : double *padfTempBuffer = nullptr;
1516 80 : GDALDataType eBufDataTypeOri = eBufDataType;
1517 80 : void *pDataBufOri = pDataBuf;
1518 : // CFloat64 is the query type used by gdallocationinfo...
1519 : #ifdef LIMIT_TYPES
1520 80 : if (eBufDataType != GDT_Byte && eBufDataType != GDT_UInt16)
1521 : #else
1522 : if (eBufDataType == GDT_CFloat64)
1523 : #endif
1524 : {
1525 43 : padfTempBuffer = static_cast<double *>(VSI_MALLOC3_VERBOSE(
1526 : nXSize, nYSize, psOptions->nOutPansharpenedBands * sizeof(double)));
1527 43 : if (padfTempBuffer == nullptr)
1528 : {
1529 0 : VSIFree(pUpsampledSpectralBuffer);
1530 0 : VSIFree(pPanBuffer);
1531 0 : return CE_Failure;
1532 : }
1533 43 : pDataBuf = padfTempBuffer;
1534 43 : eBufDataType = GDT_Float64;
1535 : }
1536 :
1537 80 : if (nTasks > 1)
1538 : {
1539 66 : std::vector<GDALPansharpenJob> asJobs;
1540 33 : asJobs.resize(nTasks);
1541 33 : GDALPansharpenJob *pasJobs = &(asJobs[0]);
1542 : {
1543 66 : std::vector<void *> ahJobData;
1544 33 : ahJobData.resize(nTasks);
1545 : #ifdef DEBUG_TIMING
1546 : struct timeval tv;
1547 : #endif
1548 165 : for (int i = 0; i < nTasks; i++)
1549 : {
1550 132 : const size_t iStartLine =
1551 132 : (static_cast<size_t>(i) * nYSize) / nTasks;
1552 132 : const size_t iNextStartLine =
1553 132 : (static_cast<size_t>(i + 1) * nYSize) / nTasks;
1554 132 : pasJobs[i].poPansharpenOperation = this;
1555 132 : pasJobs[i].eWorkDataType = eWorkDataType;
1556 132 : pasJobs[i].eBufDataType = eBufDataType;
1557 132 : pasJobs[i].pPanBuffer =
1558 132 : pPanBuffer + iStartLine * nXSize * nDataTypeSize;
1559 132 : pasJobs[i].pUpsampledSpectralBuffer =
1560 132 : pUpsampledSpectralBuffer +
1561 132 : iStartLine * nXSize * nDataTypeSize;
1562 132 : pasJobs[i].pDataBuf =
1563 132 : static_cast<GByte *>(pDataBuf) +
1564 264 : iStartLine * nXSize *
1565 132 : GDALGetDataTypeSizeBytes(eBufDataType);
1566 132 : pasJobs[i].nValues = (iNextStartLine - iStartLine) * nXSize;
1567 132 : pasJobs[i].nBandValues = static_cast<size_t>(nXSize) * nYSize;
1568 132 : pasJobs[i].nMaxValue = nMaxValue;
1569 : #ifdef DEBUG_TIMING
1570 : pasJobs[i].ptv = &tv;
1571 : #endif
1572 132 : ahJobData[i] = &(pasJobs[i]);
1573 : }
1574 : #ifdef DEBUG_TIMING
1575 : gettimeofday(&tv, nullptr);
1576 : #endif
1577 33 : poThreadPool->SubmitJobs(PansharpenJobThreadFunc, ahJobData);
1578 33 : poThreadPool->WaitCompletion();
1579 : }
1580 :
1581 33 : eErr = CE_None;
1582 165 : for (int i = 0; i < nTasks; i++)
1583 : {
1584 132 : if (pasJobs[i].eErr != CE_None)
1585 0 : eErr = CE_Failure;
1586 : }
1587 : }
1588 : else
1589 : {
1590 47 : eErr = PansharpenChunk(eWorkDataType, eBufDataType, pPanBuffer,
1591 : pUpsampledSpectralBuffer, pDataBuf,
1592 47 : static_cast<size_t>(nXSize) * nYSize,
1593 47 : static_cast<size_t>(nXSize) * nYSize, nMaxValue);
1594 : }
1595 :
1596 80 : if (padfTempBuffer)
1597 : {
1598 43 : GDALCopyWords64(padfTempBuffer, GDT_Float64, sizeof(double),
1599 : pDataBufOri, eBufDataTypeOri,
1600 : GDALGetDataTypeSizeBytes(eBufDataTypeOri),
1601 43 : static_cast<size_t>(nXSize) * nYSize *
1602 43 : psOptions->nOutPansharpenedBands);
1603 43 : VSIFree(padfTempBuffer);
1604 : }
1605 :
1606 80 : VSIFree(pUpsampledSpectralBuffer);
1607 80 : VSIFree(pPanBuffer);
1608 :
1609 80 : return eErr;
1610 : }
1611 :
1612 : /************************************************************************/
1613 : /* PansharpenResampleJobThreadFunc() */
1614 : /************************************************************************/
1615 :
1616 : // static int acc=0;
1617 :
1618 130 : void GDALPansharpenOperation::PansharpenResampleJobThreadFunc(void *pUserData)
1619 : {
1620 130 : GDALPansharpenResampleJob *psJob =
1621 : static_cast<GDALPansharpenResampleJob *>(pUserData);
1622 :
1623 : #ifdef DEBUG_TIMING
1624 : struct timeval tv;
1625 : gettimeofday(&tv, nullptr);
1626 : const GIntBig launch_time =
1627 : static_cast<GIntBig>(psJob->ptv->tv_sec) * 1000000 +
1628 : static_cast<GIntBig>(psJob->ptv->tv_usec);
1629 : const GIntBig start_job = static_cast<GIntBig>(tv.tv_sec) * 1000000 +
1630 : static_cast<GIntBig>(tv.tv_usec);
1631 : #endif
1632 :
1633 : #if 0
1634 : for(int i=0;i<1000000;i++)
1635 : acc += i * i;
1636 : #else
1637 : GDALRasterIOExtraArg sExtraArg;
1638 130 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
1639 : // cppcheck-suppress redundantAssignment
1640 130 : sExtraArg.eResampleAlg = psJob->eResampleAlg;
1641 130 : sExtraArg.bFloatingPointWindowValidity = TRUE;
1642 130 : sExtraArg.dfXOff = psJob->dfXOff;
1643 130 : sExtraArg.dfYOff = psJob->dfYOff;
1644 130 : sExtraArg.dfXSize = psJob->dfXSize;
1645 130 : sExtraArg.dfYSize = psJob->dfYSize;
1646 :
1647 130 : std::vector<int> anBands;
1648 596 : for (int i = 0; i < psJob->nBandCount; ++i)
1649 464 : anBands.push_back(i + 1);
1650 : // This call to RasterIO() in a thread to poMEMDS shared between several
1651 : // threads is really risky, but works given the implementation details...
1652 : // Do not do this at home!
1653 129 : CPL_IGNORE_RET_VAL(psJob->poMEMDS->RasterIO(
1654 : GF_Read, psJob->nXOff, psJob->nYOff, psJob->nXSize, psJob->nYSize,
1655 : psJob->pBuffer, psJob->nBufXSize, psJob->nBufYSize, psJob->eDT,
1656 132 : psJob->nBandCount, anBands.data(), 0, 0, psJob->nBandSpace,
1657 : &sExtraArg));
1658 : #endif
1659 :
1660 : #ifdef DEBUG_TIMING
1661 : struct timeval tv_end;
1662 : gettimeofday(&tv_end, nullptr);
1663 : const GIntBig end = static_cast<GIntBig>(tv_end.tv_sec) * 1000000 +
1664 : static_cast<GIntBig>(tv_end.tv_usec);
1665 : if (start_job - launch_time > 500)
1666 : /*ok*/ printf("Resample: Delay before start=" CPL_FRMT_GIB
1667 : ", completion time=" CPL_FRMT_GIB "\n",
1668 : start_job - launch_time, end - start_job);
1669 : #endif
1670 132 : }
1671 :
1672 : /************************************************************************/
1673 : /* PansharpenJobThreadFunc() */
1674 : /************************************************************************/
1675 :
1676 131 : void GDALPansharpenOperation::PansharpenJobThreadFunc(void *pUserData)
1677 : {
1678 131 : GDALPansharpenJob *psJob = static_cast<GDALPansharpenJob *>(pUserData);
1679 :
1680 : #ifdef DEBUG_TIMING
1681 : struct timeval tv;
1682 : gettimeofday(&tv, nullptr);
1683 : const GIntBig launch_time =
1684 : static_cast<GIntBig>(psJob->ptv->tv_sec) * 1000000 +
1685 : static_cast<GIntBig>(psJob->ptv->tv_usec);
1686 : const GIntBig start_job = static_cast<GIntBig>(tv.tv_sec) * 1000000 +
1687 : static_cast<GIntBig>(tv.tv_usec);
1688 : #endif
1689 :
1690 : #if 0
1691 : for( int i = 0; i < 1000000; i++ )
1692 : acc += i * i;
1693 : psJob->eErr = CE_None;
1694 : #else
1695 131 : psJob->eErr = psJob->poPansharpenOperation->PansharpenChunk(
1696 : psJob->eWorkDataType, psJob->eBufDataType, psJob->pPanBuffer,
1697 : psJob->pUpsampledSpectralBuffer, psJob->pDataBuf, psJob->nValues,
1698 : psJob->nBandValues, psJob->nMaxValue);
1699 : #endif
1700 :
1701 : #ifdef DEBUG_TIMING
1702 : struct timeval tv_end;
1703 : gettimeofday(&tv_end, nullptr);
1704 : const GIntBig end = static_cast<GIntBig>(tv_end.tv_sec) * 1000000 +
1705 : static_cast<GIntBig>(tv_end.tv_usec);
1706 : if (start_job - launch_time > 500)
1707 : /*ok*/ printf("Pansharpen: Delay before start=" CPL_FRMT_GIB
1708 : ", completion time=" CPL_FRMT_GIB "\n",
1709 : start_job - launch_time, end - start_job);
1710 : #endif
1711 132 : }
1712 :
1713 : /************************************************************************/
1714 : /* PansharpenChunk() */
1715 : /************************************************************************/
1716 :
1717 178 : CPLErr GDALPansharpenOperation::PansharpenChunk(
1718 : GDALDataType eWorkDataType, GDALDataType eBufDataType,
1719 : const void *pPanBuffer, const void *pUpsampledSpectralBuffer,
1720 : void *pDataBuf, size_t nValues, size_t nBandValues, GUInt32 nMaxValue) const
1721 : {
1722 178 : CPLErr eErr = CE_None;
1723 :
1724 178 : switch (eWorkDataType)
1725 : {
1726 65 : case GDT_Byte:
1727 65 : eErr = WeightedBrovey(
1728 : static_cast<const GByte *>(pPanBuffer),
1729 : static_cast<const GByte *>(pUpsampledSpectralBuffer), pDataBuf,
1730 : eBufDataType, nValues, nBandValues,
1731 : static_cast<GByte>(nMaxValue));
1732 65 : break;
1733 :
1734 107 : case GDT_UInt16:
1735 107 : eErr = WeightedBrovey(
1736 : static_cast<const GUInt16 *>(pPanBuffer),
1737 : static_cast<const GUInt16 *>(pUpsampledSpectralBuffer),
1738 : pDataBuf, eBufDataType, nValues, nBandValues,
1739 : static_cast<GUInt16>(nMaxValue));
1740 108 : break;
1741 :
1742 : #ifndef LIMIT_TYPES
1743 : case GDT_Int8:
1744 : eErr = WeightedBrovey(
1745 : static_cast<const GInt8 *>(pPanBuffer),
1746 : static_cast<const GInt8 *>(pUpsampledSpectralBuffer), pDataBuf,
1747 : eBufDataType, nValues, nBandValues);
1748 : break;
1749 :
1750 : case GDT_Int16:
1751 : eErr = WeightedBrovey(
1752 : static_cast<const GInt16 *>(pPanBuffer),
1753 : static_cast<const GInt16 *>(pUpsampledSpectralBuffer), pDataBuf,
1754 : eBufDataType, nValues, nBandValues);
1755 : break;
1756 :
1757 : case GDT_UInt32:
1758 : eErr = WeightedBrovey(
1759 : static_cast<const GUInt32 *>(pPanBuffer),
1760 : static_cast<const GUInt32 *>(pUpsampledSpectralBuffer),
1761 : pDataBuf, eBufDataType, nValues, nBandValues, nMaxValue);
1762 : break;
1763 :
1764 : case GDT_Int32:
1765 : eErr = WeightedBrovey(
1766 : static_cast<const GInt32 *>(pPanBuffer),
1767 : static_cast<const GInt32 *>(pUpsampledSpectralBuffer), pDataBuf,
1768 : eBufDataType, nValues, nBandValues);
1769 : break;
1770 :
1771 : case GDT_UInt64:
1772 : eErr = WeightedBrovey(
1773 : static_cast<const std::uint64_t *>(pPanBuffer),
1774 : static_cast<const std::uint64_t *>(pUpsampledSpectralBuffer),
1775 : pDataBuf, eBufDataType, nValues, nBandValues, nMaxValue);
1776 : break;
1777 :
1778 : case GDT_Int64:
1779 : eErr = WeightedBrovey(
1780 : static_cast<const std::int64_t *>(pPanBuffer),
1781 : static_cast<const std::int64_t *>(pUpsampledSpectralBuffer),
1782 : pDataBuf, eBufDataType, nValues, nBandValues);
1783 : break;
1784 :
1785 : case GDT_Float16:
1786 : eErr = WeightedBrovey(
1787 : static_cast<const GFloat16 *>(pPanBuffer),
1788 : static_cast<const GFloat16 *>(pUpsampledSpectralBuffer),
1789 : pDataBuf, eBufDataType, nValues, nBandValues);
1790 : break;
1791 :
1792 : case GDT_Float32:
1793 : eErr = WeightedBrovey(
1794 : static_cast<const float *>(pPanBuffer),
1795 : static_cast<const float *>(pUpsampledSpectralBuffer), pDataBuf,
1796 : eBufDataType, nValues, nBandValues);
1797 : break;
1798 : #endif
1799 6 : case GDT_Float64:
1800 6 : eErr = WeightedBrovey(
1801 : static_cast<const double *>(pPanBuffer),
1802 : static_cast<const double *>(pUpsampledSpectralBuffer), pDataBuf,
1803 : eBufDataType, nValues, nBandValues);
1804 6 : break;
1805 :
1806 0 : default:
1807 0 : CPLError(CE_Failure, CPLE_NotSupported,
1808 : "eWorkDataType not supported");
1809 0 : eErr = CE_Failure;
1810 0 : break;
1811 : }
1812 :
1813 179 : return eErr;
1814 : }
1815 :
1816 : /************************************************************************/
1817 : /* GetOptions() */
1818 : /************************************************************************/
1819 :
1820 : /** Return options.
1821 : * @return options.
1822 : */
1823 110 : GDALPansharpenOptions *GDALPansharpenOperation::GetOptions()
1824 : {
1825 110 : return psOptions;
1826 : }
1827 :
1828 : /************************************************************************/
1829 : /* GDALCreatePansharpenOperation() */
1830 : /************************************************************************/
1831 :
1832 : /** Instantiate a pansharpening operation.
1833 : *
1834 : * The passed options are validated.
1835 : *
1836 : * @param psOptions a pansharpening option structure allocated with
1837 : * GDALCreatePansharpenOptions(). It is duplicated by this function.
1838 : * @return a valid pansharpening operation handle, or NULL in case of failure.
1839 : *
1840 : * @since GDAL 2.1
1841 : */
1842 :
1843 : GDALPansharpenOperationH
1844 0 : GDALCreatePansharpenOperation(const GDALPansharpenOptions *psOptions)
1845 : {
1846 0 : GDALPansharpenOperation *psOperation = new GDALPansharpenOperation();
1847 0 : if (psOperation->Initialize(psOptions) == CE_None)
1848 0 : return reinterpret_cast<GDALPansharpenOperationH>(psOperation);
1849 0 : delete psOperation;
1850 0 : return nullptr;
1851 : }
1852 :
1853 : /************************************************************************/
1854 : /* GDALDestroyPansharpenOperation() */
1855 : /************************************************************************/
1856 :
1857 : /** Destroy a pansharpening operation.
1858 : *
1859 : * @param hOperation a valid pansharpening operation.
1860 : *
1861 : * @since GDAL 2.1
1862 : */
1863 :
1864 0 : void GDALDestroyPansharpenOperation(GDALPansharpenOperationH hOperation)
1865 : {
1866 0 : delete reinterpret_cast<GDALPansharpenOperation *>(hOperation);
1867 0 : }
1868 :
1869 : /************************************************************************/
1870 : /* GDALPansharpenProcessRegion() */
1871 : /************************************************************************/
1872 :
1873 : /** Executes a pansharpening operation on a rectangular region of the
1874 : * resulting dataset.
1875 : *
1876 : * The window is expressed with respect to the dimensions of the panchromatic
1877 : * band.
1878 : *
1879 : * Spectral bands are upsampled and merged with the panchromatic band according
1880 : * to the select algorithm and options.
1881 : *
1882 : * @param hOperation a valid pansharpening operation.
1883 : * @param nXOff pixel offset.
1884 : * @param nYOff pixel offset.
1885 : * @param nXSize width of the pansharpened region to compute.
1886 : * @param nYSize height of the pansharpened region to compute.
1887 : * @param pDataBuf output buffer. Must be nXSize * nYSize *
1888 : * GDALGetDataTypeSizeBytes(eBufDataType) *
1889 : * psOptions->nOutPansharpenedBands large.
1890 : * It begins with all values of the first output band, followed
1891 : * by values of the second output band, etc...
1892 : * @param eBufDataType data type of the output buffer
1893 : *
1894 : * @return CE_None in case of success, CE_Failure in case of failure.
1895 : *
1896 : * @since GDAL 2.1
1897 : */
1898 0 : CPLErr GDALPansharpenProcessRegion(GDALPansharpenOperationH hOperation,
1899 : int nXOff, int nYOff, int nXSize, int nYSize,
1900 : void *pDataBuf, GDALDataType eBufDataType)
1901 : {
1902 : return reinterpret_cast<GDALPansharpenOperation *>(hOperation)
1903 0 : ->ProcessRegion(nXOff, nYOff, nXSize, nYSize, pDataBuf, eBufDataType);
1904 : }
|
