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