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