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