Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: Implementation of a set of GDALDerivedPixelFunc(s) to be used
5 : * with source raster band of virtual GDAL datasets.
6 : * Author: Antonio Valentino <antonio.valentino@tiscali.it>
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2008-2014,2022 Antonio Valentino <antonio.valentino@tiscali.it>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : *****************************************************************************/
13 :
14 : #include <array>
15 : #include <cmath>
16 : #include "gdal.h"
17 : #include "vrtdataset.h"
18 : #include "vrtexpression.h"
19 : #include "vrtreclassifier.h"
20 : #include "cpl_float.h"
21 :
22 : #if defined(__x86_64) || defined(_M_X64) || defined(USE_NEON_OPTIMIZATIONS)
23 : #define USE_SSE2
24 : #include "gdalsse_priv.h"
25 :
26 : #if !defined(USE_NEON_OPTIMIZATIONS)
27 : #define LIBDIVIDE_SSE2
28 : #ifdef __GNUC__
29 : #pragma GCC diagnostic push
30 : #pragma GCC diagnostic ignored "-Wold-style-cast"
31 : #pragma GCC diagnostic ignored "-Weffc++"
32 : #endif
33 : #include "../../third_party/libdivide/libdivide.h"
34 : #ifdef __GNUC__
35 : #pragma GCC diagnostic pop
36 : #endif
37 : #endif
38 :
39 : #endif
40 :
41 : #include "gdal_priv_templates.hpp"
42 :
43 : #include <algorithm>
44 : #include <cassert>
45 : #include <limits>
46 :
47 : namespace gdal
48 : {
49 : MathExpression::~MathExpression() = default;
50 : }
51 :
52 : template <typename T>
53 26245100 : inline double GetSrcVal(const void *pSource, GDALDataType eSrcType, T ii)
54 : {
55 26245100 : switch (eSrcType)
56 : {
57 0 : case GDT_Unknown:
58 0 : return 0;
59 26092700 : case GDT_Byte:
60 26092700 : return static_cast<const GByte *>(pSource)[ii];
61 1876 : case GDT_Int8:
62 1876 : return static_cast<const GInt8 *>(pSource)[ii];
63 4556 : case GDT_UInt16:
64 4556 : return static_cast<const GUInt16 *>(pSource)[ii];
65 4568 : case GDT_Int16:
66 4568 : return static_cast<const GInt16 *>(pSource)[ii];
67 31892 : case GDT_UInt32:
68 31892 : return static_cast<const GUInt32 *>(pSource)[ii];
69 5806 : case GDT_Int32:
70 5806 : return static_cast<const GInt32 *>(pSource)[ii];
71 : // Precision loss currently for int64/uint64
72 1876 : case GDT_UInt64:
73 : return static_cast<double>(
74 1876 : static_cast<const uint64_t *>(pSource)[ii]);
75 1876 : case GDT_Int64:
76 : return static_cast<double>(
77 1876 : static_cast<const int64_t *>(pSource)[ii]);
78 0 : case GDT_Float16:
79 0 : return static_cast<const GFloat16 *>(pSource)[ii];
80 9924 : case GDT_Float32:
81 9924 : return static_cast<const float *>(pSource)[ii];
82 65793 : case GDT_Float64:
83 65793 : return static_cast<const double *>(pSource)[ii];
84 4352 : case GDT_CInt16:
85 4352 : return static_cast<const GInt16 *>(pSource)[2 * ii];
86 3752 : case GDT_CInt32:
87 3752 : return static_cast<const GInt32 *>(pSource)[2 * ii];
88 0 : case GDT_CFloat16:
89 0 : return static_cast<const GFloat16 *>(pSource)[2 * ii];
90 11504 : case GDT_CFloat32:
91 11504 : return static_cast<const float *>(pSource)[2 * ii];
92 4654 : case GDT_CFloat64:
93 4654 : return static_cast<const double *>(pSource)[2 * ii];
94 0 : case GDT_TypeCount:
95 0 : break;
96 : }
97 0 : return 0;
98 : }
99 :
100 10289 : static bool IsNoData(double dfVal, double dfNoData)
101 : {
102 10289 : return dfVal == dfNoData || (std::isnan(dfVal) && std::isnan(dfNoData));
103 : }
104 :
105 2324 : static CPLErr FetchDoubleArg(CSLConstList papszArgs, const char *pszName,
106 : double *pdfX, double *pdfDefault = nullptr)
107 : {
108 2324 : const char *pszVal = CSLFetchNameValue(papszArgs, pszName);
109 :
110 2324 : if (pszVal == nullptr)
111 : {
112 1134 : if (pdfDefault == nullptr)
113 : {
114 0 : CPLError(CE_Failure, CPLE_AppDefined,
115 : "Missing pixel function argument: %s", pszName);
116 0 : return CE_Failure;
117 : }
118 : else
119 : {
120 1134 : *pdfX = *pdfDefault;
121 1134 : return CE_None;
122 : }
123 : }
124 :
125 1190 : char *pszEnd = nullptr;
126 1190 : *pdfX = std::strtod(pszVal, &pszEnd);
127 1190 : if (pszEnd == pszVal)
128 : {
129 0 : CPLError(CE_Failure, CPLE_AppDefined,
130 : "Failed to parse pixel function argument: %s", pszName);
131 0 : return CE_Failure;
132 : }
133 :
134 1190 : return CE_None;
135 : }
136 :
137 7 : static CPLErr RealPixelFunc(void **papoSources, int nSources, void *pData,
138 : int nXSize, int nYSize, GDALDataType eSrcType,
139 : GDALDataType eBufType, int nPixelSpace,
140 : int nLineSpace)
141 : {
142 : /* ---- Init ---- */
143 7 : if (nSources != 1)
144 1 : return CE_Failure;
145 :
146 6 : const int nPixelSpaceSrc = GDALGetDataTypeSizeBytes(eSrcType);
147 6 : const size_t nLineSpaceSrc = static_cast<size_t>(nPixelSpaceSrc) * nXSize;
148 :
149 : /* ---- Set pixels ---- */
150 98 : for (int iLine = 0; iLine < nYSize; ++iLine)
151 : {
152 92 : GDALCopyWords(static_cast<GByte *>(papoSources[0]) +
153 92 : nLineSpaceSrc * iLine,
154 : eSrcType, nPixelSpaceSrc,
155 : static_cast<GByte *>(pData) +
156 92 : static_cast<GSpacing>(nLineSpace) * iLine,
157 : eBufType, nPixelSpace, nXSize);
158 : }
159 :
160 : /* ---- Return success ---- */
161 6 : return CE_None;
162 : } // RealPixelFunc
163 :
164 8 : static CPLErr ImagPixelFunc(void **papoSources, int nSources, void *pData,
165 : int nXSize, int nYSize, GDALDataType eSrcType,
166 : GDALDataType eBufType, int nPixelSpace,
167 : int nLineSpace)
168 : {
169 : /* ---- Init ---- */
170 8 : if (nSources != 1)
171 1 : return CE_Failure;
172 :
173 7 : if (GDALDataTypeIsComplex(eSrcType))
174 : {
175 6 : const GDALDataType eSrcBaseType = GDALGetNonComplexDataType(eSrcType);
176 6 : const int nPixelSpaceSrc = GDALGetDataTypeSizeBytes(eSrcType);
177 6 : const size_t nLineSpaceSrc =
178 6 : static_cast<size_t>(nPixelSpaceSrc) * nXSize;
179 :
180 6 : const void *const pImag = static_cast<GByte *>(papoSources[0]) +
181 6 : GDALGetDataTypeSizeBytes(eSrcType) / 2;
182 :
183 : /* ---- Set pixels ---- */
184 56 : for (int iLine = 0; iLine < nYSize; ++iLine)
185 : {
186 50 : GDALCopyWords(static_cast<const GByte *>(pImag) +
187 50 : nLineSpaceSrc * iLine,
188 : eSrcBaseType, nPixelSpaceSrc,
189 : static_cast<GByte *>(pData) +
190 50 : static_cast<GSpacing>(nLineSpace) * iLine,
191 : eBufType, nPixelSpace, nXSize);
192 : }
193 : }
194 : else
195 : {
196 1 : const double dfImag = 0;
197 :
198 : /* ---- Set pixels ---- */
199 21 : for (int iLine = 0; iLine < nYSize; ++iLine)
200 : {
201 : // Always copy from the same location.
202 20 : GDALCopyWords(&dfImag, eSrcType, 0,
203 : static_cast<GByte *>(pData) +
204 20 : static_cast<GSpacing>(nLineSpace) * iLine,
205 : eBufType, nPixelSpace, nXSize);
206 : }
207 : }
208 :
209 : /* ---- Return success ---- */
210 7 : return CE_None;
211 : } // ImagPixelFunc
212 :
213 6 : static CPLErr ComplexPixelFunc(void **papoSources, int nSources, void *pData,
214 : int nXSize, int nYSize, GDALDataType eSrcType,
215 : GDALDataType eBufType, int nPixelSpace,
216 : int nLineSpace)
217 : {
218 : /* ---- Init ---- */
219 6 : if (nSources != 2)
220 1 : return CE_Failure;
221 :
222 5 : const void *const pReal = papoSources[0];
223 5 : const void *const pImag = papoSources[1];
224 :
225 : /* ---- Set pixels ---- */
226 5 : size_t ii = 0;
227 281 : for (int iLine = 0; iLine < nYSize; ++iLine)
228 : {
229 17060 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
230 : {
231 : // Source raster pixels may be obtained with GetSrcVal macro.
232 : const double adfPixVal[2] = {
233 16784 : GetSrcVal(pReal, eSrcType, ii), // re
234 33568 : GetSrcVal(pImag, eSrcType, ii) // im
235 16784 : };
236 :
237 16784 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
238 : static_cast<GByte *>(pData) +
239 16784 : static_cast<GSpacing>(nLineSpace) * iLine +
240 16784 : iCol * nPixelSpace,
241 : eBufType, nPixelSpace, 1);
242 : }
243 : }
244 :
245 : /* ---- Return success ---- */
246 5 : return CE_None;
247 : } // ComplexPixelFunc
248 :
249 : typedef enum
250 : {
251 : GAT_amplitude,
252 : GAT_intensity,
253 : GAT_dB
254 : } PolarAmplitudeType;
255 :
256 : static const char pszPolarPixelFuncMetadata[] =
257 : "<PixelFunctionArgumentsList>"
258 : " <Argument name='amplitude_type' description='Amplitude Type' "
259 : "type='string-select' default='AMPLITUDE'>"
260 : " <Value>INTENSITY</Value>"
261 : " <Value>dB</Value>"
262 : " <Value>AMPLITUDE</Value>"
263 : " </Argument>"
264 : "</PixelFunctionArgumentsList>";
265 :
266 4 : static CPLErr PolarPixelFunc(void **papoSources, int nSources, void *pData,
267 : int nXSize, int nYSize, GDALDataType eSrcType,
268 : GDALDataType eBufType, int nPixelSpace,
269 : int nLineSpace, CSLConstList papszArgs)
270 : {
271 : /* ---- Init ---- */
272 4 : if (nSources != 2)
273 0 : return CE_Failure;
274 :
275 4 : const char pszName[] = "amplitude_type";
276 4 : const char *pszVal = CSLFetchNameValue(papszArgs, pszName);
277 4 : PolarAmplitudeType amplitudeType = GAT_amplitude;
278 4 : if (pszVal != nullptr)
279 : {
280 3 : if (strcmp(pszVal, "INTENSITY") == 0)
281 1 : amplitudeType = GAT_intensity;
282 2 : else if (strcmp(pszVal, "dB") == 0)
283 1 : amplitudeType = GAT_dB;
284 1 : else if (strcmp(pszVal, "AMPLITUDE") != 0)
285 : {
286 0 : CPLError(CE_Failure, CPLE_AppDefined,
287 : "Invalid value for pixel function argument '%s': %s",
288 : pszName, pszVal);
289 0 : return CE_Failure;
290 : }
291 : }
292 :
293 4 : const void *const pAmp = papoSources[0];
294 4 : const void *const pPhase = papoSources[1];
295 :
296 : /* ---- Set pixels ---- */
297 4 : size_t ii = 0;
298 84 : for (int iLine = 0; iLine < nYSize; ++iLine)
299 : {
300 1680 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
301 : {
302 : // Source raster pixels may be obtained with GetSrcVal macro.
303 1600 : double dfAmp = GetSrcVal(pAmp, eSrcType, ii);
304 1600 : switch (amplitudeType)
305 : {
306 400 : case GAT_intensity:
307 : // clip to zero
308 400 : dfAmp = dfAmp <= 0 ? 0 : std::sqrt(dfAmp);
309 400 : break;
310 400 : case GAT_dB:
311 400 : dfAmp = dfAmp <= 0
312 400 : ? -std::numeric_limits<double>::infinity()
313 400 : : pow(10, dfAmp / 20.);
314 400 : break;
315 800 : case GAT_amplitude:
316 800 : break;
317 : }
318 1600 : const double dfPhase = GetSrcVal(pPhase, eSrcType, ii);
319 : const double adfPixVal[2] = {
320 1600 : dfAmp * std::cos(dfPhase), // re
321 1600 : dfAmp * std::sin(dfPhase) // im
322 1600 : };
323 :
324 1600 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
325 : static_cast<GByte *>(pData) +
326 1600 : static_cast<GSpacing>(nLineSpace) * iLine +
327 1600 : iCol * nPixelSpace,
328 : eBufType, nPixelSpace, 1);
329 : }
330 : }
331 :
332 : /* ---- Return success ---- */
333 4 : return CE_None;
334 : } // PolarPixelFunc
335 :
336 11 : static CPLErr ModulePixelFunc(void **papoSources, int nSources, void *pData,
337 : int nXSize, int nYSize, GDALDataType eSrcType,
338 : GDALDataType eBufType, int nPixelSpace,
339 : int nLineSpace)
340 : {
341 : /* ---- Init ---- */
342 11 : if (nSources != 1)
343 1 : return CE_Failure;
344 :
345 10 : if (GDALDataTypeIsComplex(eSrcType))
346 : {
347 2 : const void *pReal = papoSources[0];
348 2 : const void *pImag = static_cast<GByte *>(papoSources[0]) +
349 2 : GDALGetDataTypeSizeBytes(eSrcType) / 2;
350 :
351 : /* ---- Set pixels ---- */
352 2 : size_t ii = 0;
353 14 : for (int iLine = 0; iLine < nYSize; ++iLine)
354 : {
355 72 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
356 : {
357 : // Source raster pixels may be obtained with GetSrcVal macro.
358 60 : const double dfReal = GetSrcVal(pReal, eSrcType, ii);
359 60 : const double dfImag = GetSrcVal(pImag, eSrcType, ii);
360 :
361 60 : const double dfPixVal = sqrt(dfReal * dfReal + dfImag * dfImag);
362 :
363 60 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
364 : static_cast<GByte *>(pData) +
365 60 : static_cast<GSpacing>(nLineSpace) * iLine +
366 60 : iCol * nPixelSpace,
367 : eBufType, nPixelSpace, 1);
368 : }
369 : }
370 : }
371 : else
372 : {
373 : /* ---- Set pixels ---- */
374 8 : size_t ii = 0;
375 35 : for (int iLine = 0; iLine < nYSize; ++iLine)
376 : {
377 434 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
378 : {
379 : // Source raster pixels may be obtained with GetSrcVal macro.
380 : const double dfPixVal =
381 407 : fabs(GetSrcVal(papoSources[0], eSrcType, ii));
382 :
383 407 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
384 : static_cast<GByte *>(pData) +
385 407 : static_cast<GSpacing>(nLineSpace) * iLine +
386 407 : iCol * nPixelSpace,
387 : eBufType, nPixelSpace, 1);
388 : }
389 : }
390 : }
391 :
392 : /* ---- Return success ---- */
393 10 : return CE_None;
394 : } // ModulePixelFunc
395 :
396 5 : static CPLErr PhasePixelFunc(void **papoSources, int nSources, void *pData,
397 : int nXSize, int nYSize, GDALDataType eSrcType,
398 : GDALDataType eBufType, int nPixelSpace,
399 : int nLineSpace)
400 : {
401 : /* ---- Init ---- */
402 5 : if (nSources != 1)
403 1 : return CE_Failure;
404 :
405 4 : if (GDALDataTypeIsComplex(eSrcType))
406 : {
407 2 : const void *const pReal = papoSources[0];
408 2 : const void *const pImag = static_cast<GByte *>(papoSources[0]) +
409 2 : GDALGetDataTypeSizeBytes(eSrcType) / 2;
410 :
411 : /* ---- Set pixels ---- */
412 2 : size_t ii = 0;
413 14 : for (int iLine = 0; iLine < nYSize; ++iLine)
414 : {
415 72 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
416 : {
417 : // Source raster pixels may be obtained with GetSrcVal macro.
418 60 : const double dfReal = GetSrcVal(pReal, eSrcType, ii);
419 60 : const double dfImag = GetSrcVal(pImag, eSrcType, ii);
420 :
421 60 : const double dfPixVal = atan2(dfImag, dfReal);
422 :
423 60 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
424 : static_cast<GByte *>(pData) +
425 60 : static_cast<GSpacing>(nLineSpace) * iLine +
426 60 : iCol * nPixelSpace,
427 : eBufType, nPixelSpace, 1);
428 : }
429 : }
430 : }
431 2 : else if (GDALDataTypeIsInteger(eSrcType) && !GDALDataTypeIsSigned(eSrcType))
432 : {
433 1 : constexpr double dfZero = 0;
434 7 : for (int iLine = 0; iLine < nYSize; ++iLine)
435 : {
436 6 : GDALCopyWords(&dfZero, GDT_Float64, 0,
437 : static_cast<GByte *>(pData) +
438 6 : static_cast<GSpacing>(nLineSpace) * iLine,
439 : eBufType, nPixelSpace, nXSize);
440 : }
441 : }
442 : else
443 : {
444 : /* ---- Set pixels ---- */
445 1 : size_t ii = 0;
446 7 : for (int iLine = 0; iLine < nYSize; ++iLine)
447 : {
448 36 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
449 : {
450 30 : const void *const pReal = papoSources[0];
451 :
452 : // Source raster pixels may be obtained with GetSrcVal macro.
453 30 : const double dfReal = GetSrcVal(pReal, eSrcType, ii);
454 30 : const double dfPixVal = (dfReal < 0) ? M_PI : 0.0;
455 :
456 30 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
457 : static_cast<GByte *>(pData) +
458 30 : static_cast<GSpacing>(nLineSpace) * iLine +
459 30 : iCol * nPixelSpace,
460 : eBufType, nPixelSpace, 1);
461 : }
462 : }
463 : }
464 :
465 : /* ---- Return success ---- */
466 4 : return CE_None;
467 : } // PhasePixelFunc
468 :
469 4 : static CPLErr ConjPixelFunc(void **papoSources, int nSources, void *pData,
470 : int nXSize, int nYSize, GDALDataType eSrcType,
471 : GDALDataType eBufType, int nPixelSpace,
472 : int nLineSpace)
473 : {
474 : /* ---- Init ---- */
475 4 : if (nSources != 1)
476 1 : return CE_Failure;
477 :
478 3 : if (GDALDataTypeIsComplex(eSrcType) && GDALDataTypeIsComplex(eBufType))
479 : {
480 2 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
481 2 : const void *const pReal = papoSources[0];
482 2 : const void *const pImag =
483 2 : static_cast<GByte *>(papoSources[0]) + nOffset;
484 :
485 : /* ---- Set pixels ---- */
486 2 : size_t ii = 0;
487 14 : for (int iLine = 0; iLine < nYSize; ++iLine)
488 : {
489 72 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
490 : {
491 : // Source raster pixels may be obtained with GetSrcVal macro.
492 : const double adfPixVal[2] = {
493 60 : +GetSrcVal(pReal, eSrcType, ii), // re
494 120 : -GetSrcVal(pImag, eSrcType, ii) // im
495 60 : };
496 :
497 60 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
498 : static_cast<GByte *>(pData) +
499 60 : static_cast<GSpacing>(nLineSpace) * iLine +
500 60 : iCol * nPixelSpace,
501 : eBufType, nPixelSpace, 1);
502 : }
503 : }
504 : }
505 : else
506 : {
507 : // No complex data type.
508 1 : return RealPixelFunc(papoSources, nSources, pData, nXSize, nYSize,
509 1 : eSrcType, eBufType, nPixelSpace, nLineSpace);
510 : }
511 :
512 : /* ---- Return success ---- */
513 2 : return CE_None;
514 : } // ConjPixelFunc
515 :
516 : #ifdef USE_SSE2
517 :
518 : /************************************************************************/
519 : /* OptimizedSumToFloat_SSE2() */
520 : /************************************************************************/
521 :
522 : template <typename Tsrc>
523 87 : static void OptimizedSumToFloat_SSE2(double dfK, void *pOutBuffer,
524 : int nLineSpace, int nXSize, int nYSize,
525 : int nSources,
526 : const void *const *papoSources)
527 : {
528 87 : const XMMReg4Float cst = XMMReg4Float::Set1(static_cast<float>(dfK));
529 :
530 279 : for (int iLine = 0; iLine < nYSize; ++iLine)
531 : {
532 192 : float *CPL_RESTRICT const pDest = reinterpret_cast<float *>(
533 : static_cast<GByte *>(pOutBuffer) +
534 192 : static_cast<GSpacing>(nLineSpace) * iLine);
535 192 : const size_t iOffsetLine = static_cast<size_t>(iLine) * nXSize;
536 :
537 192 : constexpr int VALUES_PER_REG = 4;
538 192 : constexpr int UNROLLING = 4 * VALUES_PER_REG;
539 192 : int iCol = 0;
540 900 : for (; iCol < nXSize - (UNROLLING - 1); iCol += UNROLLING)
541 : {
542 708 : XMMReg4Float d0(cst);
543 708 : XMMReg4Float d1(cst);
544 708 : XMMReg4Float d2(cst);
545 708 : XMMReg4Float d3(cst);
546 2104 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
547 : {
548 1396 : XMMReg4Float t0, t1, t2, t3;
549 1396 : XMMReg4Float::Load16Val(
550 1396 : static_cast<const Tsrc * CPL_RESTRICT>(papoSources[iSrc]) +
551 1396 : iOffsetLine + iCol,
552 : t0, t1, t2, t3);
553 1396 : d0 += t0;
554 1396 : d1 += t1;
555 1396 : d2 += t2;
556 1396 : d3 += t3;
557 : }
558 708 : d0.Store4Val(pDest + iCol + VALUES_PER_REG * 0);
559 708 : d1.Store4Val(pDest + iCol + VALUES_PER_REG * 1);
560 708 : d2.Store4Val(pDest + iCol + VALUES_PER_REG * 2);
561 708 : d3.Store4Val(pDest + iCol + VALUES_PER_REG * 3);
562 : }
563 :
564 788 : for (; iCol < nXSize; iCol++)
565 : {
566 596 : float d = static_cast<float>(dfK);
567 1708 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
568 : {
569 1112 : d += static_cast<const Tsrc * CPL_RESTRICT>(
570 1112 : papoSources[iSrc])[iOffsetLine + iCol];
571 : }
572 596 : pDest[iCol] = d;
573 : }
574 : }
575 87 : }
576 :
577 : /************************************************************************/
578 : /* OptimizedSumToDouble_SSE2() */
579 : /************************************************************************/
580 :
581 : template <typename Tsrc>
582 107 : static void OptimizedSumToDouble_SSE2(double dfK, void *pOutBuffer,
583 : int nLineSpace, int nXSize, int nYSize,
584 : int nSources,
585 : const void *const *papoSources)
586 : {
587 107 : const XMMReg4Double cst = XMMReg4Double::Set1(dfK);
588 :
589 355 : for (int iLine = 0; iLine < nYSize; ++iLine)
590 : {
591 248 : double *CPL_RESTRICT const pDest = reinterpret_cast<double *>(
592 : static_cast<GByte *>(pOutBuffer) +
593 248 : static_cast<GSpacing>(nLineSpace) * iLine);
594 248 : const size_t iOffsetLine = static_cast<size_t>(iLine) * nXSize;
595 :
596 248 : constexpr int VALUES_PER_REG = 4;
597 248 : constexpr int UNROLLING = 2 * VALUES_PER_REG;
598 248 : int iCol = 0;
599 1976 : for (; iCol < nXSize - (UNROLLING - 1); iCol += UNROLLING)
600 : {
601 1728 : XMMReg4Double d0(cst);
602 1728 : XMMReg4Double d1(cst);
603 5104 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
604 : {
605 3376 : XMMReg4Double t0, t1;
606 3376 : XMMReg4Double::Load8Val(
607 3376 : static_cast<const Tsrc * CPL_RESTRICT>(papoSources[iSrc]) +
608 3376 : iOffsetLine + iCol,
609 : t0, t1);
610 3376 : d0 += t0;
611 3376 : d1 += t1;
612 : }
613 1728 : d0.Store4Val(pDest + iCol + VALUES_PER_REG * 0);
614 1728 : d1.Store4Val(pDest + iCol + VALUES_PER_REG * 1);
615 : }
616 :
617 1028 : for (; iCol < nXSize; iCol++)
618 : {
619 780 : double d = dfK;
620 2180 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
621 : {
622 1400 : d += static_cast<const Tsrc * CPL_RESTRICT>(
623 1400 : papoSources[iSrc])[iOffsetLine + iCol];
624 : }
625 780 : pDest[iCol] = d;
626 : }
627 : }
628 107 : }
629 :
630 : /************************************************************************/
631 : /* OptimizedSumSameType_SSE2() */
632 : /************************************************************************/
633 :
634 : template <typename T, typename Tsigned, typename Tacc, class SSEWrapper>
635 1011 : static void OptimizedSumSameType_SSE2(double dfK, void *pOutBuffer,
636 : int nLineSpace, int nXSize, int nYSize,
637 : int nSources,
638 : const void *const *papoSources)
639 : {
640 : static_assert(std::numeric_limits<T>::is_integer);
641 : static_assert(!std::numeric_limits<T>::is_signed);
642 : static_assert(std::numeric_limits<Tsigned>::is_integer);
643 : static_assert(std::numeric_limits<Tsigned>::is_signed);
644 : static_assert(sizeof(T) == sizeof(Tsigned));
645 1011 : const T nK = static_cast<T>(dfK);
646 : Tsigned nKSigned;
647 1011 : memcpy(&nKSigned, &nK, sizeof(T));
648 1011 : const __m128i valInit = SSEWrapper::Set1(nKSigned);
649 1011 : constexpr int VALUES_PER_REG =
650 : static_cast<int>(sizeof(valInit) / sizeof(T));
651 2030 : for (int iLine = 0; iLine < nYSize; ++iLine)
652 : {
653 1019 : T *CPL_RESTRICT const pDest =
654 : reinterpret_cast<T *>(static_cast<GByte *>(pOutBuffer) +
655 1019 : static_cast<GSpacing>(nLineSpace) * iLine);
656 1019 : const size_t iOffsetLine = static_cast<size_t>(iLine) * nXSize;
657 1019 : int iCol = 0;
658 8043 : for (; iCol < nXSize - (4 * VALUES_PER_REG - 1);
659 : iCol += 4 * VALUES_PER_REG)
660 : {
661 7024 : __m128i reg0 = valInit;
662 7024 : __m128i reg1 = valInit;
663 7024 : __m128i reg2 = valInit;
664 7024 : __m128i reg3 = valInit;
665 21072 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
666 : {
667 14048 : reg0 = SSEWrapper::AddSaturate(
668 : reg0,
669 : _mm_loadu_si128(reinterpret_cast<const __m128i *>(
670 14048 : static_cast<const T * CPL_RESTRICT>(papoSources[iSrc]) +
671 14048 : iOffsetLine + iCol)));
672 14048 : reg1 = SSEWrapper::AddSaturate(
673 : reg1,
674 : _mm_loadu_si128(reinterpret_cast<const __m128i *>(
675 14048 : static_cast<const T * CPL_RESTRICT>(papoSources[iSrc]) +
676 14048 : iOffsetLine + iCol + VALUES_PER_REG)));
677 14048 : reg2 = SSEWrapper::AddSaturate(
678 : reg2,
679 : _mm_loadu_si128(reinterpret_cast<const __m128i *>(
680 14048 : static_cast<const T * CPL_RESTRICT>(papoSources[iSrc]) +
681 14048 : iOffsetLine + iCol + 2 * VALUES_PER_REG)));
682 14048 : reg3 = SSEWrapper::AddSaturate(
683 : reg3,
684 : _mm_loadu_si128(reinterpret_cast<const __m128i *>(
685 14048 : static_cast<const T * CPL_RESTRICT>(papoSources[iSrc]) +
686 14048 : iOffsetLine + iCol + 3 * VALUES_PER_REG)));
687 : }
688 7024 : _mm_storeu_si128(reinterpret_cast<__m128i *>(pDest + iCol), reg0);
689 7024 : _mm_storeu_si128(
690 7024 : reinterpret_cast<__m128i *>(pDest + iCol + VALUES_PER_REG),
691 : reg1);
692 7024 : _mm_storeu_si128(
693 7024 : reinterpret_cast<__m128i *>(pDest + iCol + 2 * VALUES_PER_REG),
694 : reg2);
695 7024 : _mm_storeu_si128(
696 7024 : reinterpret_cast<__m128i *>(pDest + iCol + 3 * VALUES_PER_REG),
697 : reg3);
698 : }
699 53070 : for (; iCol < nXSize; ++iCol)
700 : {
701 52051 : Tacc nAcc = nK;
702 156151 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
703 : {
704 208200 : nAcc = std::min<Tacc>(
705 208200 : nAcc + static_cast<const T * CPL_RESTRICT>(
706 104100 : papoSources[iSrc])[iOffsetLine + iCol],
707 104100 : std::numeric_limits<T>::max());
708 : }
709 52051 : pDest[iCol] = static_cast<T>(nAcc);
710 : }
711 : }
712 1011 : }
713 : #endif // USE_SSE2
714 :
715 : /************************************************************************/
716 : /* OptimizedSumPackedOutput() */
717 : /************************************************************************/
718 :
719 : template <typename Tsrc, typename Tdest>
720 230 : static void OptimizedSumPackedOutput(double dfK, void *pOutBuffer,
721 : int nLineSpace, int nXSize, int nYSize,
722 : int nSources,
723 : const void *const *papoSources)
724 : {
725 : #ifdef USE_SSE2
726 : if constexpr (std::is_same_v<Tdest, float> && !std::is_same_v<Tsrc, double>)
727 : {
728 87 : OptimizedSumToFloat_SSE2<Tsrc>(dfK, pOutBuffer, nLineSpace, nXSize,
729 : nYSize, nSources, papoSources);
730 : }
731 : else if constexpr (std::is_same_v<Tdest, double>)
732 : {
733 107 : OptimizedSumToDouble_SSE2<Tsrc>(dfK, pOutBuffer, nLineSpace, nXSize,
734 : nYSize, nSources, papoSources);
735 : }
736 : else
737 : #endif // USE_SSE2
738 : {
739 36 : const Tdest nCst = static_cast<Tdest>(dfK);
740 144 : for (int iLine = 0; iLine < nYSize; ++iLine)
741 : {
742 108 : Tdest *CPL_RESTRICT const pDest = reinterpret_cast<Tdest *>(
743 : static_cast<GByte *>(pOutBuffer) +
744 108 : static_cast<GSpacing>(nLineSpace) * iLine);
745 108 : const size_t iOffsetLine = static_cast<size_t>(iLine) * nXSize;
746 :
747 : #define LOAD_SRCVAL(iSrc_, j_) \
748 : static_cast<Tdest>(static_cast<const Tsrc * CPL_RESTRICT>( \
749 : papoSources[(iSrc_)])[iOffsetLine + iCol + (j_)])
750 :
751 108 : constexpr int UNROLLING = 4;
752 108 : int iCol = 0;
753 1396 : for (; iCol < nXSize - (UNROLLING - 1); iCol += UNROLLING)
754 : {
755 1288 : Tdest d[4] = {nCst, nCst, nCst, nCst};
756 4064 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
757 : {
758 2776 : d[0] += LOAD_SRCVAL(iSrc, 0);
759 2776 : d[1] += LOAD_SRCVAL(iSrc, 1);
760 2776 : d[2] += LOAD_SRCVAL(iSrc, 2);
761 2776 : d[3] += LOAD_SRCVAL(iSrc, 3);
762 : }
763 1288 : pDest[iCol + 0] = d[0];
764 1288 : pDest[iCol + 1] = d[1];
765 1288 : pDest[iCol + 2] = d[2];
766 1288 : pDest[iCol + 3] = d[3];
767 : }
768 312 : for (; iCol < nXSize; iCol++)
769 : {
770 204 : Tdest d0 = nCst;
771 612 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
772 : {
773 408 : d0 += LOAD_SRCVAL(iSrc, 0);
774 : }
775 204 : pDest[iCol] = d0;
776 : }
777 : #undef LOAD_SRCVAL
778 : }
779 : }
780 230 : }
781 :
782 : /************************************************************************/
783 : /* OptimizedSumPackedOutput() */
784 : /************************************************************************/
785 :
786 : template <typename Tdest>
787 253 : static bool OptimizedSumPackedOutput(GDALDataType eSrcType, double dfK,
788 : void *pOutBuffer, int nLineSpace,
789 : int nXSize, int nYSize, int nSources,
790 : const void *const *papoSources)
791 : {
792 253 : switch (eSrcType)
793 : {
794 39 : case GDT_Byte:
795 39 : OptimizedSumPackedOutput<uint8_t, Tdest>(dfK, pOutBuffer,
796 : nLineSpace, nXSize, nYSize,
797 : nSources, papoSources);
798 39 : return true;
799 :
800 34 : case GDT_UInt16:
801 34 : OptimizedSumPackedOutput<uint16_t, Tdest>(
802 : dfK, pOutBuffer, nLineSpace, nXSize, nYSize, nSources,
803 : papoSources);
804 34 : return true;
805 :
806 34 : case GDT_Int16:
807 34 : OptimizedSumPackedOutput<int16_t, Tdest>(dfK, pOutBuffer,
808 : nLineSpace, nXSize, nYSize,
809 : nSources, papoSources);
810 34 : return true;
811 :
812 34 : case GDT_Int32:
813 34 : OptimizedSumPackedOutput<int32_t, Tdest>(dfK, pOutBuffer,
814 : nLineSpace, nXSize, nYSize,
815 : nSources, papoSources);
816 34 : return true;
817 :
818 36 : case GDT_Float32:
819 36 : OptimizedSumPackedOutput<float, Tdest>(dfK, pOutBuffer, nLineSpace,
820 : nXSize, nYSize, nSources,
821 : papoSources);
822 36 : return true;
823 :
824 36 : case GDT_Float64:
825 36 : OptimizedSumPackedOutput<double, Tdest>(dfK, pOutBuffer, nLineSpace,
826 : nXSize, nYSize, nSources,
827 : papoSources);
828 36 : return true;
829 :
830 40 : default:
831 40 : break;
832 : }
833 40 : return false;
834 : }
835 :
836 : /************************************************************************/
837 : /* OptimizedSumThroughLargerType() */
838 : /************************************************************************/
839 :
840 : namespace
841 : {
842 : template <typename Tsrc, typename Tdest, typename Enable = void>
843 : struct TintermediateS
844 : {
845 : using type = double;
846 : };
847 :
848 : template <typename Tsrc, typename Tdest>
849 : struct TintermediateS<
850 : Tsrc, Tdest,
851 : std::enable_if_t<
852 : (std::is_same_v<Tsrc, uint8_t> || std::is_same_v<Tsrc, int16_t> ||
853 : std::is_same_v<Tsrc, uint16_t>)&&(std::is_same_v<Tdest, uint8_t> ||
854 : std::is_same_v<Tdest, int16_t> ||
855 : std::is_same_v<Tdest, uint16_t>),
856 : bool>>
857 : {
858 : using type = int32_t;
859 : };
860 :
861 : } // namespace
862 :
863 : template <typename Tsrc, typename Tdest>
864 387 : static bool OptimizedSumThroughLargerType(double dfK, void *pOutBuffer,
865 : int nPixelSpace, int nLineSpace,
866 : int nXSize, int nYSize, int nSources,
867 : const void *const *papoSources)
868 : {
869 : using Tintermediate = typename TintermediateS<Tsrc, Tdest>::type;
870 387 : const Tintermediate k = static_cast<Tintermediate>(dfK);
871 :
872 387 : size_t ii = 0;
873 1161 : for (int iLine = 0; iLine < nYSize; ++iLine)
874 : {
875 774 : GByte *CPL_RESTRICT pDest = static_cast<GByte *>(pOutBuffer) +
876 774 : static_cast<GSpacing>(nLineSpace) * iLine;
877 :
878 774 : constexpr int UNROLLING = 4;
879 774 : int iCol = 0;
880 13158 : for (; iCol < nXSize - (UNROLLING - 1);
881 : iCol += UNROLLING, ii += UNROLLING)
882 : {
883 12384 : Tintermediate aSum[4] = {k, k, k, k};
884 :
885 37152 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
886 : {
887 24768 : aSum[0] += static_cast<const Tsrc *>(papoSources[iSrc])[ii + 0];
888 24768 : aSum[1] += static_cast<const Tsrc *>(papoSources[iSrc])[ii + 1];
889 24768 : aSum[2] += static_cast<const Tsrc *>(papoSources[iSrc])[ii + 2];
890 24768 : aSum[3] += static_cast<const Tsrc *>(papoSources[iSrc])[ii + 3];
891 : }
892 :
893 12384 : GDALCopyWord(aSum[0], *reinterpret_cast<Tdest *>(pDest));
894 12384 : pDest += nPixelSpace;
895 12384 : GDALCopyWord(aSum[1], *reinterpret_cast<Tdest *>(pDest));
896 12384 : pDest += nPixelSpace;
897 12384 : GDALCopyWord(aSum[2], *reinterpret_cast<Tdest *>(pDest));
898 12384 : pDest += nPixelSpace;
899 12384 : GDALCopyWord(aSum[3], *reinterpret_cast<Tdest *>(pDest));
900 12384 : pDest += nPixelSpace;
901 : }
902 3096 : for (; iCol < nXSize; ++iCol, ++ii, pDest += nPixelSpace)
903 : {
904 2322 : Tintermediate sum = k;
905 6966 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
906 : {
907 4644 : sum += static_cast<const Tsrc *>(papoSources[iSrc])[ii];
908 : }
909 :
910 2322 : auto pDst = reinterpret_cast<Tdest *>(pDest);
911 2322 : GDALCopyWord(sum, *pDst);
912 : }
913 : }
914 387 : return true;
915 : }
916 :
917 : /************************************************************************/
918 : /* OptimizedSumThroughLargerType() */
919 : /************************************************************************/
920 :
921 : template <typename Tsrc>
922 490 : static bool OptimizedSumThroughLargerType(GDALDataType eBufType, double dfK,
923 : void *pOutBuffer, int nPixelSpace,
924 : int nLineSpace, int nXSize,
925 : int nYSize, int nSources,
926 : const void *const *papoSources)
927 : {
928 490 : switch (eBufType)
929 : {
930 99 : case GDT_Byte:
931 99 : return OptimizedSumThroughLargerType<Tsrc, uint8_t>(
932 : dfK, pOutBuffer, nPixelSpace, nLineSpace, nXSize, nYSize,
933 99 : nSources, papoSources);
934 :
935 99 : case GDT_UInt16:
936 99 : return OptimizedSumThroughLargerType<Tsrc, uint16_t>(
937 : dfK, pOutBuffer, nPixelSpace, nLineSpace, nXSize, nYSize,
938 99 : nSources, papoSources);
939 :
940 103 : case GDT_Int16:
941 103 : return OptimizedSumThroughLargerType<Tsrc, int16_t>(
942 : dfK, pOutBuffer, nPixelSpace, nLineSpace, nXSize, nYSize,
943 103 : nSources, papoSources);
944 :
945 86 : case GDT_Int32:
946 86 : return OptimizedSumThroughLargerType<Tsrc, int32_t>(
947 : dfK, pOutBuffer, nPixelSpace, nLineSpace, nXSize, nYSize,
948 86 : nSources, papoSources);
949 :
950 : // Float32 and Float64 already covered by OptimizedSum() for packed case
951 103 : default:
952 103 : break;
953 : }
954 103 : return false;
955 : }
956 :
957 : /************************************************************************/
958 : /* OptimizedSumThroughLargerType() */
959 : /************************************************************************/
960 :
961 630 : static bool OptimizedSumThroughLargerType(GDALDataType eSrcType,
962 : GDALDataType eBufType, double dfK,
963 : void *pOutBuffer, int nPixelSpace,
964 : int nLineSpace, int nXSize,
965 : int nYSize, int nSources,
966 : const void *const *papoSources)
967 : {
968 630 : switch (eSrcType)
969 : {
970 64 : case GDT_Byte:
971 64 : return OptimizedSumThroughLargerType<uint8_t>(
972 : eBufType, dfK, pOutBuffer, nPixelSpace, nLineSpace, nXSize,
973 64 : nYSize, nSources, papoSources);
974 :
975 81 : case GDT_UInt16:
976 81 : return OptimizedSumThroughLargerType<uint16_t>(
977 : eBufType, dfK, pOutBuffer, nPixelSpace, nLineSpace, nXSize,
978 81 : nYSize, nSources, papoSources);
979 :
980 85 : case GDT_Int16:
981 85 : return OptimizedSumThroughLargerType<int16_t>(
982 : eBufType, dfK, pOutBuffer, nPixelSpace, nLineSpace, nXSize,
983 85 : nYSize, nSources, papoSources);
984 :
985 85 : case GDT_Int32:
986 85 : return OptimizedSumThroughLargerType<int32_t>(
987 : eBufType, dfK, pOutBuffer, nPixelSpace, nLineSpace, nXSize,
988 85 : nYSize, nSources, papoSources);
989 :
990 90 : case GDT_Float32:
991 90 : return OptimizedSumThroughLargerType<float>(
992 : eBufType, dfK, pOutBuffer, nPixelSpace, nLineSpace, nXSize,
993 90 : nYSize, nSources, papoSources);
994 :
995 85 : case GDT_Float64:
996 85 : return OptimizedSumThroughLargerType<double>(
997 : eBufType, dfK, pOutBuffer, nPixelSpace, nLineSpace, nXSize,
998 85 : nYSize, nSources, papoSources);
999 :
1000 140 : default:
1001 140 : break;
1002 : }
1003 :
1004 140 : return false;
1005 : }
1006 :
1007 : /************************************************************************/
1008 : /* SumPixelFunc() */
1009 : /************************************************************************/
1010 :
1011 : static const char pszSumPixelFuncMetadata[] =
1012 : "<PixelFunctionArgumentsList>"
1013 : " <Argument name='k' description='Optional constant term' type='double' "
1014 : "default='0.0' />"
1015 : " <Argument name='propagateNoData' description='Whether the output value "
1016 : "should be NoData as as soon as one source is NoData' type='boolean' "
1017 : "default='false' />"
1018 : " <Argument type='builtin' value='NoData' optional='true' />"
1019 : "</PixelFunctionArgumentsList>";
1020 :
1021 1990 : static CPLErr SumPixelFunc(void **papoSources, int nSources, void *pData,
1022 : int nXSize, int nYSize, GDALDataType eSrcType,
1023 : GDALDataType eBufType, int nPixelSpace,
1024 : int nLineSpace, CSLConstList papszArgs)
1025 : {
1026 : /* ---- Init ---- */
1027 1990 : if (nSources < 1)
1028 : {
1029 1 : CPLError(CE_Failure, CPLE_AppDefined,
1030 : "sum requires at least one source");
1031 1 : return CE_Failure;
1032 : }
1033 :
1034 1989 : double dfK = 0.0;
1035 1989 : if (FetchDoubleArg(papszArgs, "k", &dfK, &dfK) != CE_None)
1036 0 : return CE_Failure;
1037 :
1038 1989 : double dfNoData{0};
1039 1989 : bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1040 1989 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1041 0 : return CE_Failure;
1042 :
1043 1989 : const bool bPropagateNoData = CPLTestBool(
1044 : CSLFetchNameValueDef(papszArgs, "propagateNoData", "false"));
1045 :
1046 1989 : if (dfNoData == 0 && !bPropagateNoData)
1047 1907 : bHasNoData = false;
1048 :
1049 : /* ---- Set pixels ---- */
1050 1989 : if (GDALDataTypeIsComplex(eSrcType))
1051 : {
1052 36 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
1053 :
1054 : /* ---- Set pixels ---- */
1055 36 : size_t ii = 0;
1056 112 : for (int iLine = 0; iLine < nYSize; ++iLine)
1057 : {
1058 4796 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1059 : {
1060 4720 : double adfSum[2] = {dfK, 0.0};
1061 :
1062 14190 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
1063 : {
1064 9470 : const void *const pReal = papoSources[iSrc];
1065 9470 : const void *const pImag =
1066 9470 : static_cast<const GByte *>(pReal) + nOffset;
1067 :
1068 : // Source raster pixels may be obtained with GetSrcVal
1069 : // macro.
1070 9470 : adfSum[0] += GetSrcVal(pReal, eSrcType, ii);
1071 9470 : adfSum[1] += GetSrcVal(pImag, eSrcType, ii);
1072 : }
1073 :
1074 4720 : GDALCopyWords(adfSum, GDT_CFloat64, 0,
1075 : static_cast<GByte *>(pData) +
1076 4720 : static_cast<GSpacing>(nLineSpace) * iLine +
1077 4720 : iCol * nPixelSpace,
1078 : eBufType, nPixelSpace, 1);
1079 : }
1080 : }
1081 : }
1082 : else
1083 : {
1084 : /* ---- Set pixels ---- */
1085 1953 : bool bGeneralCase = true;
1086 1953 : if (dfNoData == 0 && !bPropagateNoData)
1087 : {
1088 : #ifdef USE_SSE2
1089 1126 : if (eBufType == GDT_Byte && nPixelSpace == sizeof(uint8_t) &&
1090 1020 : eSrcType == GDT_Byte &&
1091 1020 : dfK >= std::numeric_limits<uint8_t>::min() &&
1092 4017 : dfK <= std::numeric_limits<uint8_t>::max() &&
1093 1020 : static_cast<int>(dfK) == dfK)
1094 : {
1095 1007 : bGeneralCase = false;
1096 :
1097 : struct SSEWrapper
1098 : {
1099 1007 : inline static __m128i Set1(int8_t x)
1100 : {
1101 2014 : return _mm_set1_epi8(x);
1102 : }
1103 :
1104 56064 : inline static __m128i AddSaturate(__m128i x, __m128i y)
1105 : {
1106 56064 : return _mm_adds_epu8(x, y);
1107 : }
1108 : };
1109 :
1110 : OptimizedSumSameType_SSE2<uint8_t, int8_t, uint32_t,
1111 1007 : SSEWrapper>(dfK, pData, nLineSpace,
1112 : nXSize, nYSize, nSources,
1113 : papoSources);
1114 : }
1115 123 : else if (eBufType == GDT_UInt16 &&
1116 123 : nPixelSpace == sizeof(uint16_t) &&
1117 17 : eSrcType == GDT_UInt16 &&
1118 17 : dfK >= std::numeric_limits<uint16_t>::min() &&
1119 1004 : dfK <= std::numeric_limits<uint16_t>::max() &&
1120 17 : static_cast<int>(dfK) == dfK)
1121 : {
1122 4 : bGeneralCase = false;
1123 :
1124 : struct SSEWrapper
1125 : {
1126 4 : inline static __m128i Set1(int16_t x)
1127 : {
1128 8 : return _mm_set1_epi16(x);
1129 : }
1130 :
1131 128 : inline static __m128i AddSaturate(__m128i x, __m128i y)
1132 : {
1133 128 : return _mm_adds_epu16(x, y);
1134 : }
1135 : };
1136 :
1137 : OptimizedSumSameType_SSE2<uint16_t, int16_t, uint32_t,
1138 4 : SSEWrapper>(dfK, pData, nLineSpace,
1139 : nXSize, nYSize, nSources,
1140 : papoSources);
1141 : }
1142 : else
1143 : #endif
1144 860 : if (eBufType == GDT_Float32 && nPixelSpace == sizeof(float))
1145 : {
1146 126 : bGeneralCase = !OptimizedSumPackedOutput<float>(
1147 : eSrcType, dfK, pData, nLineSpace, nXSize, nYSize, nSources,
1148 : papoSources);
1149 : }
1150 734 : else if (eBufType == GDT_Float64 && nPixelSpace == sizeof(double))
1151 : {
1152 127 : bGeneralCase = !OptimizedSumPackedOutput<double>(
1153 : eSrcType, dfK, pData, nLineSpace, nXSize, nYSize, nSources,
1154 : papoSources);
1155 : }
1156 607 : else if (
1157 607 : dfK >= 0 && dfK <= INT_MAX && eBufType == GDT_Int32 &&
1158 1231 : nPixelSpace == sizeof(int32_t) && eSrcType == GDT_Byte &&
1159 : // Limitation to avoid overflow of int32 if all source values are at the max of their data type
1160 17 : nSources <=
1161 17 : (INT_MAX - dfK) / std::numeric_limits<uint8_t>::max())
1162 : {
1163 17 : bGeneralCase = false;
1164 17 : OptimizedSumPackedOutput<uint8_t, int32_t>(
1165 : dfK, pData, nLineSpace, nXSize, nYSize, nSources,
1166 : papoSources);
1167 : }
1168 :
1169 2501 : if (bGeneralCase && dfK >= 0 && dfK <= INT_MAX &&
1170 630 : nSources <=
1171 630 : (INT_MAX - dfK) / std::numeric_limits<uint16_t>::max())
1172 : {
1173 630 : bGeneralCase = !OptimizedSumThroughLargerType(
1174 : eSrcType, eBufType, dfK, pData, nPixelSpace, nLineSpace,
1175 : nXSize, nYSize, nSources, papoSources);
1176 : }
1177 : }
1178 :
1179 1953 : if (bGeneralCase)
1180 : {
1181 325 : size_t ii = 0;
1182 1287 : for (int iLine = 0; iLine < nYSize; ++iLine)
1183 : {
1184 53629 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1185 : {
1186 52667 : double dfSum = dfK;
1187 157938 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
1188 : {
1189 : const double dfVal =
1190 105280 : GetSrcVal(papoSources[iSrc], eSrcType, ii);
1191 :
1192 105280 : if (bHasNoData && IsNoData(dfVal, dfNoData))
1193 : {
1194 21 : if (bPropagateNoData)
1195 : {
1196 9 : dfSum = dfNoData;
1197 9 : break;
1198 : }
1199 : }
1200 : else
1201 : {
1202 105259 : dfSum += dfVal;
1203 : }
1204 : }
1205 :
1206 52667 : GDALCopyWords(&dfSum, GDT_Float64, 0,
1207 : static_cast<GByte *>(pData) +
1208 52667 : static_cast<GSpacing>(nLineSpace) *
1209 52667 : iLine +
1210 52667 : iCol * nPixelSpace,
1211 : eBufType, nPixelSpace, 1);
1212 : }
1213 : }
1214 : }
1215 : }
1216 :
1217 : /* ---- Return success ---- */
1218 1989 : return CE_None;
1219 : } /* SumPixelFunc */
1220 :
1221 : static const char pszDiffPixelFuncMetadata[] =
1222 : "<PixelFunctionArgumentsList>"
1223 : " <Argument type='builtin' value='NoData' optional='true' />"
1224 : "</PixelFunctionArgumentsList>";
1225 :
1226 12 : static CPLErr DiffPixelFunc(void **papoSources, int nSources, void *pData,
1227 : int nXSize, int nYSize, GDALDataType eSrcType,
1228 : GDALDataType eBufType, int nPixelSpace,
1229 : int nLineSpace, CSLConstList papszArgs)
1230 : {
1231 : /* ---- Init ---- */
1232 12 : if (nSources != 2)
1233 1 : return CE_Failure;
1234 :
1235 11 : double dfNoData{0};
1236 11 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1237 11 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1238 0 : return CE_Failure;
1239 :
1240 11 : if (GDALDataTypeIsComplex(eSrcType))
1241 : {
1242 1 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
1243 1 : const void *const pReal0 = papoSources[0];
1244 1 : const void *const pImag0 =
1245 1 : static_cast<GByte *>(papoSources[0]) + nOffset;
1246 1 : const void *const pReal1 = papoSources[1];
1247 1 : const void *const pImag1 =
1248 1 : static_cast<GByte *>(papoSources[1]) + nOffset;
1249 :
1250 : /* ---- Set pixels ---- */
1251 1 : size_t ii = 0;
1252 7 : for (int iLine = 0; iLine < nYSize; ++iLine)
1253 : {
1254 36 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1255 : {
1256 : // Source raster pixels may be obtained with GetSrcVal macro.
1257 30 : double adfPixVal[2] = {GetSrcVal(pReal0, eSrcType, ii) -
1258 30 : GetSrcVal(pReal1, eSrcType, ii),
1259 90 : GetSrcVal(pImag0, eSrcType, ii) -
1260 30 : GetSrcVal(pImag1, eSrcType, ii)};
1261 :
1262 30 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
1263 : static_cast<GByte *>(pData) +
1264 30 : static_cast<GSpacing>(nLineSpace) * iLine +
1265 30 : iCol * nPixelSpace,
1266 : eBufType, nPixelSpace, 1);
1267 : }
1268 : }
1269 : }
1270 : else
1271 : {
1272 : /* ---- Set pixels ---- */
1273 10 : size_t ii = 0;
1274 25 : for (int iLine = 0; iLine < nYSize; ++iLine)
1275 : {
1276 57 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1277 : {
1278 42 : const double dfA = GetSrcVal(papoSources[0], eSrcType, ii);
1279 42 : const double dfB = GetSrcVal(papoSources[1], eSrcType, ii);
1280 :
1281 : const double dfPixVal =
1282 46 : bHasNoData &&
1283 4 : (IsNoData(dfA, dfNoData) || IsNoData(dfB, dfNoData))
1284 46 : ? dfNoData
1285 42 : : dfA - dfB;
1286 :
1287 42 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1288 : static_cast<GByte *>(pData) +
1289 42 : static_cast<GSpacing>(nLineSpace) * iLine +
1290 42 : iCol * nPixelSpace,
1291 : eBufType, nPixelSpace, 1);
1292 : }
1293 : }
1294 : }
1295 :
1296 : /* ---- Return success ---- */
1297 11 : return CE_None;
1298 : } // DiffPixelFunc
1299 :
1300 : static const char pszMulPixelFuncMetadata[] =
1301 : "<PixelFunctionArgumentsList>"
1302 : " <Argument name='k' description='Optional constant factor' "
1303 : "type='double' default='1.0' />"
1304 : " <Argument name='propagateNoData' description='Whether the output value "
1305 : "should be NoData as as soon as one source is NoData' type='boolean' "
1306 : "default='false' />"
1307 : " <Argument type='builtin' value='NoData' optional='true' />"
1308 : "</PixelFunctionArgumentsList>";
1309 :
1310 68 : static CPLErr MulPixelFunc(void **papoSources, int nSources, void *pData,
1311 : int nXSize, int nYSize, GDALDataType eSrcType,
1312 : GDALDataType eBufType, int nPixelSpace,
1313 : int nLineSpace, CSLConstList papszArgs)
1314 : {
1315 : /* ---- Init ---- */
1316 68 : if (nSources < 2 && CSLFetchNameValue(papszArgs, "k") == nullptr)
1317 : {
1318 1 : CPLError(CE_Failure, CPLE_AppDefined,
1319 : "mul requires at least two sources or a specified constant k");
1320 1 : return CE_Failure;
1321 : }
1322 :
1323 67 : double dfK = 1.0;
1324 67 : if (FetchDoubleArg(papszArgs, "k", &dfK, &dfK) != CE_None)
1325 0 : return CE_Failure;
1326 :
1327 67 : double dfNoData{0};
1328 67 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1329 67 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1330 0 : return CE_Failure;
1331 :
1332 67 : const bool bPropagateNoData = CPLTestBool(
1333 : CSLFetchNameValueDef(papszArgs, "propagateNoData", "false"));
1334 :
1335 : /* ---- Set pixels ---- */
1336 67 : if (GDALDataTypeIsComplex(eSrcType))
1337 : {
1338 1 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
1339 :
1340 : /* ---- Set pixels ---- */
1341 1 : size_t ii = 0;
1342 7 : for (int iLine = 0; iLine < nYSize; ++iLine)
1343 : {
1344 36 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1345 : {
1346 30 : double adfPixVal[2] = {dfK, 0.0};
1347 :
1348 90 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
1349 : {
1350 60 : const void *const pReal = papoSources[iSrc];
1351 60 : const void *const pImag =
1352 60 : static_cast<const GByte *>(pReal) + nOffset;
1353 :
1354 60 : const double dfOldR = adfPixVal[0];
1355 60 : const double dfOldI = adfPixVal[1];
1356 :
1357 : // Source raster pixels may be obtained with GetSrcVal
1358 : // macro.
1359 60 : const double dfNewR = GetSrcVal(pReal, eSrcType, ii);
1360 60 : const double dfNewI = GetSrcVal(pImag, eSrcType, ii);
1361 :
1362 60 : adfPixVal[0] = dfOldR * dfNewR - dfOldI * dfNewI;
1363 60 : adfPixVal[1] = dfOldR * dfNewI + dfOldI * dfNewR;
1364 : }
1365 :
1366 30 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
1367 : static_cast<GByte *>(pData) +
1368 30 : static_cast<GSpacing>(nLineSpace) * iLine +
1369 30 : iCol * nPixelSpace,
1370 : eBufType, nPixelSpace, 1);
1371 : }
1372 : }
1373 : }
1374 : else
1375 : {
1376 : /* ---- Set pixels ---- */
1377 66 : size_t ii = 0;
1378 679 : for (int iLine = 0; iLine < nYSize; ++iLine)
1379 : {
1380 26883 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1381 : {
1382 26270 : double dfPixVal = dfK; // Not complex.
1383 :
1384 54157 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
1385 : {
1386 : const double dfVal =
1387 27893 : GetSrcVal(papoSources[iSrc], eSrcType, ii);
1388 :
1389 27893 : if (bHasNoData && IsNoData(dfVal, dfNoData))
1390 : {
1391 18 : if (bPropagateNoData)
1392 : {
1393 6 : dfPixVal = dfNoData;
1394 6 : break;
1395 : }
1396 : }
1397 : else
1398 : {
1399 27875 : dfPixVal *= dfVal;
1400 : }
1401 : }
1402 :
1403 26270 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1404 : static_cast<GByte *>(pData) +
1405 26270 : static_cast<GSpacing>(nLineSpace) * iLine +
1406 26270 : iCol * nPixelSpace,
1407 : eBufType, nPixelSpace, 1);
1408 : }
1409 : }
1410 : }
1411 :
1412 : /* ---- Return success ---- */
1413 67 : return CE_None;
1414 : } // MulPixelFunc
1415 :
1416 : static const char pszDivPixelFuncMetadata[] =
1417 : "<PixelFunctionArgumentsList>"
1418 : " "
1419 : "<Argument type='builtin' value='NoData' optional='true' />"
1420 : "</PixelFunctionArgumentsList>";
1421 :
1422 14 : static CPLErr DivPixelFunc(void **papoSources, int nSources, void *pData,
1423 : int nXSize, int nYSize, GDALDataType eSrcType,
1424 : GDALDataType eBufType, int nPixelSpace,
1425 : int nLineSpace, CSLConstList papszArgs)
1426 : {
1427 : /* ---- Init ---- */
1428 14 : if (nSources != 2)
1429 0 : return CE_Failure;
1430 :
1431 14 : double dfNoData{0};
1432 14 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1433 14 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1434 0 : return CE_Failure;
1435 :
1436 : /* ---- Set pixels ---- */
1437 14 : if (GDALDataTypeIsComplex(eSrcType))
1438 : {
1439 1 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
1440 1 : const void *const pReal0 = papoSources[0];
1441 1 : const void *const pImag0 =
1442 1 : static_cast<GByte *>(papoSources[0]) + nOffset;
1443 1 : const void *const pReal1 = papoSources[1];
1444 1 : const void *const pImag1 =
1445 1 : static_cast<GByte *>(papoSources[1]) + nOffset;
1446 :
1447 : /* ---- Set pixels ---- */
1448 1 : size_t ii = 0;
1449 7 : for (int iLine = 0; iLine < nYSize; ++iLine)
1450 : {
1451 36 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1452 : {
1453 : // Source raster pixels may be obtained with GetSrcVal macro.
1454 30 : const double dfReal0 = GetSrcVal(pReal0, eSrcType, ii);
1455 30 : const double dfReal1 = GetSrcVal(pReal1, eSrcType, ii);
1456 30 : const double dfImag0 = GetSrcVal(pImag0, eSrcType, ii);
1457 30 : const double dfImag1 = GetSrcVal(pImag1, eSrcType, ii);
1458 30 : const double dfAux = dfReal1 * dfReal1 + dfImag1 * dfImag1;
1459 :
1460 : const double adfPixVal[2] = {
1461 : dfAux == 0
1462 30 : ? std::numeric_limits<double>::infinity()
1463 30 : : dfReal0 * dfReal1 / dfAux + dfImag0 * dfImag1 / dfAux,
1464 0 : dfAux == 0 ? std::numeric_limits<double>::infinity()
1465 30 : : dfReal1 / dfAux * dfImag0 -
1466 30 : dfReal0 * dfImag1 / dfAux};
1467 :
1468 30 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
1469 : static_cast<GByte *>(pData) +
1470 30 : static_cast<GSpacing>(nLineSpace) * iLine +
1471 30 : iCol * nPixelSpace,
1472 : eBufType, nPixelSpace, 1);
1473 : }
1474 : }
1475 : }
1476 : else
1477 : {
1478 : /* ---- Set pixels ---- */
1479 13 : size_t ii = 0;
1480 31 : for (int iLine = 0; iLine < nYSize; ++iLine)
1481 : {
1482 63 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1483 : {
1484 45 : const double dfNum = GetSrcVal(papoSources[0], eSrcType, ii);
1485 45 : const double dfDenom = GetSrcVal(papoSources[1], eSrcType, ii);
1486 :
1487 45 : double dfPixVal = dfNoData;
1488 49 : if (!bHasNoData || (!IsNoData(dfNum, dfNoData) &&
1489 4 : !IsNoData(dfDenom, dfNoData)))
1490 : {
1491 : // coverity[divide_by_zero]
1492 41 : dfPixVal =
1493 : dfDenom == 0
1494 41 : ? std::numeric_limits<double>::infinity()
1495 : : dfNum /
1496 : #ifdef __COVERITY__
1497 : (dfDenom + std::numeric_limits<double>::min())
1498 : #else
1499 : dfDenom
1500 : #endif
1501 : ;
1502 : }
1503 :
1504 45 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1505 : static_cast<GByte *>(pData) +
1506 45 : static_cast<GSpacing>(nLineSpace) * iLine +
1507 45 : iCol * nPixelSpace,
1508 : eBufType, nPixelSpace, 1);
1509 : }
1510 : }
1511 : }
1512 :
1513 : /* ---- Return success ---- */
1514 14 : return CE_None;
1515 : } // DivPixelFunc
1516 :
1517 3 : static CPLErr CMulPixelFunc(void **papoSources, int nSources, void *pData,
1518 : int nXSize, int nYSize, GDALDataType eSrcType,
1519 : GDALDataType eBufType, int nPixelSpace,
1520 : int nLineSpace)
1521 : {
1522 : /* ---- Init ---- */
1523 3 : if (nSources != 2)
1524 1 : return CE_Failure;
1525 :
1526 : /* ---- Set pixels ---- */
1527 2 : if (GDALDataTypeIsComplex(eSrcType))
1528 : {
1529 1 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
1530 1 : const void *const pReal0 = papoSources[0];
1531 1 : const void *const pImag0 =
1532 1 : static_cast<GByte *>(papoSources[0]) + nOffset;
1533 1 : const void *const pReal1 = papoSources[1];
1534 1 : const void *const pImag1 =
1535 1 : static_cast<GByte *>(papoSources[1]) + nOffset;
1536 :
1537 1 : size_t ii = 0;
1538 7 : for (int iLine = 0; iLine < nYSize; ++iLine)
1539 : {
1540 36 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1541 : {
1542 : // Source raster pixels may be obtained with GetSrcVal macro.
1543 30 : const double dfReal0 = GetSrcVal(pReal0, eSrcType, ii);
1544 30 : const double dfReal1 = GetSrcVal(pReal1, eSrcType, ii);
1545 30 : const double dfImag0 = GetSrcVal(pImag0, eSrcType, ii);
1546 30 : const double dfImag1 = GetSrcVal(pImag1, eSrcType, ii);
1547 : const double adfPixVal[2] = {
1548 30 : dfReal0 * dfReal1 + dfImag0 * dfImag1,
1549 30 : dfReal1 * dfImag0 - dfReal0 * dfImag1};
1550 :
1551 30 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
1552 : static_cast<GByte *>(pData) +
1553 30 : static_cast<GSpacing>(nLineSpace) * iLine +
1554 30 : iCol * nPixelSpace,
1555 : eBufType, nPixelSpace, 1);
1556 : }
1557 : }
1558 : }
1559 : else
1560 : {
1561 1 : size_t ii = 0;
1562 21 : for (int iLine = 0; iLine < nYSize; ++iLine)
1563 : {
1564 420 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1565 : {
1566 : // Source raster pixels may be obtained with GetSrcVal macro.
1567 : // Not complex.
1568 400 : const double adfPixVal[2] = {
1569 400 : GetSrcVal(papoSources[0], eSrcType, ii) *
1570 400 : GetSrcVal(papoSources[1], eSrcType, ii),
1571 400 : 0.0};
1572 :
1573 400 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
1574 : static_cast<GByte *>(pData) +
1575 400 : static_cast<GSpacing>(nLineSpace) * iLine +
1576 400 : iCol * nPixelSpace,
1577 : eBufType, nPixelSpace, 1);
1578 : }
1579 : }
1580 : }
1581 :
1582 : /* ---- Return success ---- */
1583 2 : return CE_None;
1584 : } // CMulPixelFunc
1585 :
1586 : static const char pszInvPixelFuncMetadata[] =
1587 : "<PixelFunctionArgumentsList>"
1588 : " <Argument name='k' description='Optional constant factor' "
1589 : "type='double' default='1.0' />"
1590 : " "
1591 : "<Argument type='builtin' value='NoData' optional='true' />"
1592 : "</PixelFunctionArgumentsList>";
1593 :
1594 13 : static CPLErr InvPixelFunc(void **papoSources, int nSources, void *pData,
1595 : int nXSize, int nYSize, GDALDataType eSrcType,
1596 : GDALDataType eBufType, int nPixelSpace,
1597 : int nLineSpace, CSLConstList papszArgs)
1598 : {
1599 : /* ---- Init ---- */
1600 13 : if (nSources != 1)
1601 1 : return CE_Failure;
1602 :
1603 12 : double dfK = 1.0;
1604 12 : if (FetchDoubleArg(papszArgs, "k", &dfK, &dfK) != CE_None)
1605 0 : return CE_Failure;
1606 :
1607 12 : double dfNoData{0};
1608 12 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1609 12 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1610 0 : return CE_Failure;
1611 :
1612 : /* ---- Set pixels ---- */
1613 12 : if (GDALDataTypeIsComplex(eSrcType))
1614 : {
1615 2 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
1616 2 : const void *const pReal = papoSources[0];
1617 2 : const void *const pImag =
1618 2 : static_cast<GByte *>(papoSources[0]) + nOffset;
1619 :
1620 2 : size_t ii = 0;
1621 9 : for (int iLine = 0; iLine < nYSize; ++iLine)
1622 : {
1623 38 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1624 : {
1625 : // Source raster pixels may be obtained with GetSrcVal macro.
1626 31 : const double dfReal = GetSrcVal(pReal, eSrcType, ii);
1627 31 : const double dfImag = GetSrcVal(pImag, eSrcType, ii);
1628 31 : const double dfAux = dfReal * dfReal + dfImag * dfImag;
1629 : const double adfPixVal[2] = {
1630 31 : dfAux == 0 ? std::numeric_limits<double>::infinity()
1631 30 : : dfK * dfReal / dfAux,
1632 1 : dfAux == 0 ? std::numeric_limits<double>::infinity()
1633 31 : : -dfK * dfImag / dfAux};
1634 :
1635 31 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
1636 : static_cast<GByte *>(pData) +
1637 31 : static_cast<GSpacing>(nLineSpace) * iLine +
1638 31 : iCol * nPixelSpace,
1639 : eBufType, nPixelSpace, 1);
1640 : }
1641 : }
1642 : }
1643 : else
1644 : {
1645 : /* ---- Set pixels ---- */
1646 10 : size_t ii = 0;
1647 58 : for (int iLine = 0; iLine < nYSize; ++iLine)
1648 : {
1649 860 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1650 : {
1651 : // Source raster pixels may be obtained with GetSrcVal macro.
1652 : // Not complex.
1653 812 : const double dfVal = GetSrcVal(papoSources[0], eSrcType, ii);
1654 812 : double dfPixVal = dfNoData;
1655 :
1656 812 : if (!bHasNoData || !IsNoData(dfVal, dfNoData))
1657 : {
1658 807 : dfPixVal =
1659 : dfVal == 0
1660 807 : ? std::numeric_limits<double>::infinity()
1661 806 : : dfK /
1662 : #ifdef __COVERITY__
1663 : (dfVal + std::numeric_limits<double>::min())
1664 : #else
1665 : dfVal
1666 : #endif
1667 : ;
1668 : }
1669 :
1670 812 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1671 : static_cast<GByte *>(pData) +
1672 812 : static_cast<GSpacing>(nLineSpace) * iLine +
1673 812 : iCol * nPixelSpace,
1674 : eBufType, nPixelSpace, 1);
1675 : }
1676 : }
1677 : }
1678 :
1679 : /* ---- Return success ---- */
1680 12 : return CE_None;
1681 : } // InvPixelFunc
1682 :
1683 4 : static CPLErr IntensityPixelFunc(void **papoSources, int nSources, void *pData,
1684 : int nXSize, int nYSize, GDALDataType eSrcType,
1685 : GDALDataType eBufType, int nPixelSpace,
1686 : int nLineSpace)
1687 : {
1688 : /* ---- Init ---- */
1689 4 : if (nSources != 1)
1690 1 : return CE_Failure;
1691 :
1692 3 : if (GDALDataTypeIsComplex(eSrcType))
1693 : {
1694 2 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
1695 2 : const void *const pReal = papoSources[0];
1696 2 : const void *const pImag =
1697 2 : static_cast<GByte *>(papoSources[0]) + nOffset;
1698 :
1699 : /* ---- Set pixels ---- */
1700 2 : size_t ii = 0;
1701 14 : for (int iLine = 0; iLine < nYSize; ++iLine)
1702 : {
1703 72 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1704 : {
1705 : // Source raster pixels may be obtained with GetSrcVal macro.
1706 60 : const double dfReal = GetSrcVal(pReal, eSrcType, ii);
1707 60 : const double dfImag = GetSrcVal(pImag, eSrcType, ii);
1708 :
1709 60 : const double dfPixVal = dfReal * dfReal + dfImag * dfImag;
1710 :
1711 60 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1712 : static_cast<GByte *>(pData) +
1713 60 : static_cast<GSpacing>(nLineSpace) * iLine +
1714 60 : iCol * nPixelSpace,
1715 : eBufType, nPixelSpace, 1);
1716 : }
1717 : }
1718 : }
1719 : else
1720 : {
1721 : /* ---- Set pixels ---- */
1722 1 : size_t ii = 0;
1723 21 : for (int iLine = 0; iLine < nYSize; ++iLine)
1724 : {
1725 420 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1726 : {
1727 : // Source raster pixels may be obtained with GetSrcVal macro.
1728 400 : double dfPixVal = GetSrcVal(papoSources[0], eSrcType, ii);
1729 400 : dfPixVal *= dfPixVal;
1730 :
1731 400 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1732 : static_cast<GByte *>(pData) +
1733 400 : static_cast<GSpacing>(nLineSpace) * iLine +
1734 400 : iCol * nPixelSpace,
1735 : eBufType, nPixelSpace, 1);
1736 : }
1737 : }
1738 : }
1739 :
1740 : /* ---- Return success ---- */
1741 3 : return CE_None;
1742 : } // IntensityPixelFunc
1743 :
1744 : static const char pszSqrtPixelFuncMetadata[] =
1745 : "<PixelFunctionArgumentsList>"
1746 : " <Argument type='builtin' value='NoData' optional='true'/>"
1747 : "</PixelFunctionArgumentsList>";
1748 :
1749 7 : static CPLErr SqrtPixelFunc(void **papoSources, int nSources, void *pData,
1750 : int nXSize, int nYSize, GDALDataType eSrcType,
1751 : GDALDataType eBufType, int nPixelSpace,
1752 : int nLineSpace, CSLConstList papszArgs)
1753 : {
1754 : /* ---- Init ---- */
1755 7 : if (nSources != 1)
1756 1 : return CE_Failure;
1757 6 : if (GDALDataTypeIsComplex(eSrcType))
1758 0 : return CE_Failure;
1759 :
1760 6 : double dfNoData{0};
1761 6 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1762 6 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1763 0 : return CE_Failure;
1764 :
1765 : /* ---- Set pixels ---- */
1766 6 : size_t ii = 0;
1767 31 : for (int iLine = 0; iLine < nYSize; ++iLine)
1768 : {
1769 431 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1770 : {
1771 : // Source raster pixels may be obtained with GetSrcVal macro.
1772 406 : double dfPixVal = GetSrcVal(papoSources[0], eSrcType, ii);
1773 :
1774 406 : if (bHasNoData && IsNoData(dfPixVal, dfNoData))
1775 : {
1776 2 : dfPixVal = dfNoData;
1777 : }
1778 : else
1779 : {
1780 404 : dfPixVal = std::sqrt(dfPixVal);
1781 : }
1782 :
1783 406 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1784 : static_cast<GByte *>(pData) +
1785 406 : static_cast<GSpacing>(nLineSpace) * iLine +
1786 406 : iCol * nPixelSpace,
1787 : eBufType, nPixelSpace, 1);
1788 : }
1789 : }
1790 :
1791 : /* ---- Return success ---- */
1792 6 : return CE_None;
1793 : } // SqrtPixelFunc
1794 :
1795 17 : static CPLErr Log10PixelFuncHelper(void **papoSources, int nSources,
1796 : void *pData, int nXSize, int nYSize,
1797 : GDALDataType eSrcType, GDALDataType eBufType,
1798 : int nPixelSpace, int nLineSpace,
1799 : CSLConstList papszArgs, double fact)
1800 : {
1801 : /* ---- Init ---- */
1802 17 : if (nSources != 1)
1803 2 : return CE_Failure;
1804 :
1805 15 : double dfNoData{0};
1806 15 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1807 15 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1808 0 : return CE_Failure;
1809 :
1810 15 : if (GDALDataTypeIsComplex(eSrcType))
1811 : {
1812 : // Complex input datatype.
1813 5 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
1814 5 : const void *const pReal = papoSources[0];
1815 5 : const void *const pImag =
1816 5 : static_cast<GByte *>(papoSources[0]) + nOffset;
1817 :
1818 : /* We should compute fact * log10( sqrt( dfReal * dfReal + dfImag *
1819 : * dfImag ) ) */
1820 : /* Given that log10(sqrt(x)) = 0.5 * log10(x) */
1821 : /* we can remove the sqrt() by multiplying fact by 0.5 */
1822 5 : fact *= 0.5;
1823 :
1824 : /* ---- Set pixels ---- */
1825 5 : size_t ii = 0;
1826 35 : for (int iLine = 0; iLine < nYSize; ++iLine)
1827 : {
1828 180 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1829 : {
1830 : // Source raster pixels may be obtained with GetSrcVal macro.
1831 150 : const double dfReal = GetSrcVal(pReal, eSrcType, ii);
1832 150 : const double dfImag = GetSrcVal(pImag, eSrcType, ii);
1833 :
1834 : const double dfPixVal =
1835 150 : fact * log10(dfReal * dfReal + dfImag * dfImag);
1836 :
1837 150 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1838 : static_cast<GByte *>(pData) +
1839 150 : static_cast<GSpacing>(nLineSpace) * iLine +
1840 150 : iCol * nPixelSpace,
1841 : eBufType, nPixelSpace, 1);
1842 : }
1843 : }
1844 : }
1845 : else
1846 : {
1847 : /* ---- Set pixels ---- */
1848 10 : size_t ii = 0;
1849 96 : for (int iLine = 0; iLine < nYSize; ++iLine)
1850 : {
1851 1694 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1852 : {
1853 : // Source raster pixels may be obtained with GetSrcVal macro.
1854 1608 : const double dfSrcVal = GetSrcVal(papoSources[0], eSrcType, ii);
1855 : const double dfPixVal =
1856 4 : bHasNoData && IsNoData(dfSrcVal, dfNoData)
1857 1612 : ? dfNoData
1858 1604 : : fact * std::log10(std::abs(dfSrcVal));
1859 :
1860 1608 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1861 : static_cast<GByte *>(pData) +
1862 1608 : static_cast<GSpacing>(nLineSpace) * iLine +
1863 1608 : iCol * nPixelSpace,
1864 : eBufType, nPixelSpace, 1);
1865 : }
1866 : }
1867 : }
1868 :
1869 : /* ---- Return success ---- */
1870 15 : return CE_None;
1871 : } // Log10PixelFuncHelper
1872 :
1873 : static const char pszLog10PixelFuncMetadata[] =
1874 : "<PixelFunctionArgumentsList>"
1875 : " <Argument type='builtin' value='NoData' optional='true'/>"
1876 : "</PixelFunctionArgumentsList>";
1877 :
1878 9 : static CPLErr Log10PixelFunc(void **papoSources, int nSources, void *pData,
1879 : int nXSize, int nYSize, GDALDataType eSrcType,
1880 : GDALDataType eBufType, int nPixelSpace,
1881 : int nLineSpace, CSLConstList papszArgs)
1882 : {
1883 9 : return Log10PixelFuncHelper(papoSources, nSources, pData, nXSize, nYSize,
1884 : eSrcType, eBufType, nPixelSpace, nLineSpace,
1885 9 : papszArgs, 1.0);
1886 : } // Log10PixelFunc
1887 :
1888 : static const char pszDBPixelFuncMetadata[] =
1889 : "<PixelFunctionArgumentsList>"
1890 : " <Argument name='fact' description='Factor' type='double' "
1891 : "default='20.0' />"
1892 : " <Argument type='builtin' value='NoData' optional='true' />"
1893 : "</PixelFunctionArgumentsList>";
1894 :
1895 8 : static CPLErr DBPixelFunc(void **papoSources, int nSources, void *pData,
1896 : int nXSize, int nYSize, GDALDataType eSrcType,
1897 : GDALDataType eBufType, int nPixelSpace,
1898 : int nLineSpace, CSLConstList papszArgs)
1899 : {
1900 8 : double dfFact = 20.;
1901 8 : if (FetchDoubleArg(papszArgs, "fact", &dfFact, &dfFact) != CE_None)
1902 0 : return CE_Failure;
1903 :
1904 8 : return Log10PixelFuncHelper(papoSources, nSources, pData, nXSize, nYSize,
1905 : eSrcType, eBufType, nPixelSpace, nLineSpace,
1906 8 : papszArgs, dfFact);
1907 : } // DBPixelFunc
1908 :
1909 12 : static CPLErr ExpPixelFuncHelper(void **papoSources, int nSources, void *pData,
1910 : int nXSize, int nYSize, GDALDataType eSrcType,
1911 : GDALDataType eBufType, int nPixelSpace,
1912 : int nLineSpace, CSLConstList papszArgs,
1913 : double base, double fact)
1914 : {
1915 : /* ---- Init ---- */
1916 12 : if (nSources != 1)
1917 2 : return CE_Failure;
1918 10 : if (GDALDataTypeIsComplex(eSrcType))
1919 0 : return CE_Failure;
1920 :
1921 10 : double dfNoData{0};
1922 10 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1923 10 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1924 0 : return CE_Failure;
1925 :
1926 : /* ---- Set pixels ---- */
1927 10 : size_t ii = 0;
1928 115 : for (int iLine = 0; iLine < nYSize; ++iLine)
1929 : {
1930 2111 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1931 : {
1932 : // Source raster pixels may be obtained with GetSrcVal macro.
1933 2006 : const double dfVal = GetSrcVal(papoSources[0], eSrcType, ii);
1934 2 : const double dfPixVal = bHasNoData && IsNoData(dfVal, dfNoData)
1935 2008 : ? dfNoData
1936 2004 : : pow(base, dfVal * fact);
1937 :
1938 2006 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1939 : static_cast<GByte *>(pData) +
1940 2006 : static_cast<GSpacing>(nLineSpace) * iLine +
1941 2006 : iCol * nPixelSpace,
1942 : eBufType, nPixelSpace, 1);
1943 : }
1944 : }
1945 :
1946 : /* ---- Return success ---- */
1947 10 : return CE_None;
1948 : } // ExpPixelFuncHelper
1949 :
1950 : static const char pszExpPixelFuncMetadata[] =
1951 : "<PixelFunctionArgumentsList>"
1952 : " <Argument name='base' description='Base' type='double' "
1953 : "default='2.7182818284590452353602874713526624' />"
1954 : " <Argument name='fact' description='Factor' type='double' default='1' />"
1955 : " <Argument type='builtin' value='NoData' optional='true' />"
1956 : "</PixelFunctionArgumentsList>";
1957 :
1958 8 : static CPLErr ExpPixelFunc(void **papoSources, int nSources, void *pData,
1959 : int nXSize, int nYSize, GDALDataType eSrcType,
1960 : GDALDataType eBufType, int nPixelSpace,
1961 : int nLineSpace, CSLConstList papszArgs)
1962 : {
1963 8 : double dfBase = 2.7182818284590452353602874713526624;
1964 8 : double dfFact = 1.;
1965 :
1966 8 : if (FetchDoubleArg(papszArgs, "base", &dfBase, &dfBase) != CE_None)
1967 0 : return CE_Failure;
1968 :
1969 8 : if (FetchDoubleArg(papszArgs, "fact", &dfFact, &dfFact) != CE_None)
1970 0 : return CE_Failure;
1971 :
1972 8 : return ExpPixelFuncHelper(papoSources, nSources, pData, nXSize, nYSize,
1973 : eSrcType, eBufType, nPixelSpace, nLineSpace,
1974 8 : papszArgs, dfBase, dfFact);
1975 : } // ExpPixelFunc
1976 :
1977 2 : static CPLErr dB2AmpPixelFunc(void **papoSources, int nSources, void *pData,
1978 : int nXSize, int nYSize, GDALDataType eSrcType,
1979 : GDALDataType eBufType, int nPixelSpace,
1980 : int nLineSpace)
1981 : {
1982 2 : return ExpPixelFuncHelper(papoSources, nSources, pData, nXSize, nYSize,
1983 : eSrcType, eBufType, nPixelSpace, nLineSpace,
1984 2 : nullptr, 10.0, 1. / 20);
1985 : } // dB2AmpPixelFunc
1986 :
1987 2 : static CPLErr dB2PowPixelFunc(void **papoSources, int nSources, void *pData,
1988 : int nXSize, int nYSize, GDALDataType eSrcType,
1989 : GDALDataType eBufType, int nPixelSpace,
1990 : int nLineSpace)
1991 : {
1992 2 : return ExpPixelFuncHelper(papoSources, nSources, pData, nXSize, nYSize,
1993 : eSrcType, eBufType, nPixelSpace, nLineSpace,
1994 2 : nullptr, 10.0, 1. / 10);
1995 : } // dB2PowPixelFunc
1996 :
1997 : static const char pszPowPixelFuncMetadata[] =
1998 : "<PixelFunctionArgumentsList>"
1999 : " <Argument name='power' description='Exponent' type='double' "
2000 : "mandatory='1' />"
2001 : " <Argument type='builtin' value='NoData' optional='true' />"
2002 : "</PixelFunctionArgumentsList>";
2003 :
2004 6 : static CPLErr PowPixelFunc(void **papoSources, int nSources, void *pData,
2005 : int nXSize, int nYSize, GDALDataType eSrcType,
2006 : GDALDataType eBufType, int nPixelSpace,
2007 : int nLineSpace, CSLConstList papszArgs)
2008 : {
2009 : /* ---- Init ---- */
2010 6 : if (nSources != 1)
2011 0 : return CE_Failure;
2012 6 : if (GDALDataTypeIsComplex(eSrcType))
2013 0 : return CE_Failure;
2014 :
2015 : double power;
2016 6 : if (FetchDoubleArg(papszArgs, "power", &power) != CE_None)
2017 0 : return CE_Failure;
2018 :
2019 6 : double dfNoData{0};
2020 6 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
2021 6 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
2022 0 : return CE_Failure;
2023 :
2024 : /* ---- Set pixels ---- */
2025 6 : size_t ii = 0;
2026 31 : for (int iLine = 0; iLine < nYSize; ++iLine)
2027 : {
2028 431 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
2029 : {
2030 406 : const double dfVal = GetSrcVal(papoSources[0], eSrcType, ii);
2031 :
2032 2 : const double dfPixVal = bHasNoData && IsNoData(dfVal, dfNoData)
2033 408 : ? dfNoData
2034 404 : : std::pow(dfVal, power);
2035 :
2036 406 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
2037 : static_cast<GByte *>(pData) +
2038 406 : static_cast<GSpacing>(nLineSpace) * iLine +
2039 406 : iCol * nPixelSpace,
2040 : eBufType, nPixelSpace, 1);
2041 : }
2042 : }
2043 :
2044 : /* ---- Return success ---- */
2045 6 : return CE_None;
2046 : }
2047 :
2048 : // Given nt intervals spaced by dt and beginning at t0, return the index of
2049 : // the lower bound of the interval that should be used to
2050 : // interpolate/extrapolate a value for t.
2051 17 : static std::size_t intervalLeft(double t0, double dt, std::size_t nt, double t)
2052 : {
2053 17 : if (t < t0)
2054 : {
2055 4 : return 0;
2056 : }
2057 :
2058 13 : std::size_t n = static_cast<std::size_t>((t - t0) / dt);
2059 :
2060 13 : if (n >= nt - 1)
2061 : {
2062 3 : return nt - 2;
2063 : }
2064 :
2065 10 : return n;
2066 : }
2067 :
2068 17 : static double InterpolateLinear(double dfX0, double dfX1, double dfY0,
2069 : double dfY1, double dfX)
2070 : {
2071 17 : return dfY0 + (dfX - dfX0) * (dfY1 - dfY0) / (dfX1 - dfX0);
2072 : }
2073 :
2074 13 : static double InterpolateExponential(double dfX0, double dfX1, double dfY0,
2075 : double dfY1, double dfX)
2076 : {
2077 13 : const double r = std::log(dfY1 / dfY0) / (dfX1 - dfX0);
2078 13 : return dfY0 * std::exp(r * (dfX - dfX0));
2079 : }
2080 :
2081 : static const char pszInterpolatePixelFuncMetadata[] =
2082 : "<PixelFunctionArgumentsList>"
2083 : " <Argument name='t0' description='t0' type='double' mandatory='1' />"
2084 : " <Argument name='dt' description='dt' type='double' mandatory='1' />"
2085 : " <Argument name='t' description='t' type='double' mandatory='1' />"
2086 : " <Argument type='builtin' value='NoData' optional='true' />"
2087 : "</PixelFunctionArgumentsList>";
2088 :
2089 : template <decltype(InterpolateLinear) InterpolationFunction>
2090 17 : CPLErr InterpolatePixelFunc(void **papoSources, int nSources, void *pData,
2091 : int nXSize, int nYSize, GDALDataType eSrcType,
2092 : GDALDataType eBufType, int nPixelSpace,
2093 : int nLineSpace, CSLConstList papszArgs)
2094 : {
2095 : /* ---- Init ---- */
2096 17 : if (GDALDataTypeIsComplex(eSrcType))
2097 0 : return CE_Failure;
2098 :
2099 : double dfT0;
2100 17 : if (FetchDoubleArg(papszArgs, "t0", &dfT0) == CE_Failure)
2101 0 : return CE_Failure;
2102 :
2103 : double dfT;
2104 17 : if (FetchDoubleArg(papszArgs, "t", &dfT) == CE_Failure)
2105 0 : return CE_Failure;
2106 :
2107 : double dfDt;
2108 17 : if (FetchDoubleArg(papszArgs, "dt", &dfDt) == CE_Failure)
2109 0 : return CE_Failure;
2110 :
2111 17 : double dfNoData{0};
2112 17 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
2113 17 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
2114 0 : return CE_Failure;
2115 :
2116 17 : if (nSources < 2)
2117 : {
2118 0 : CPLError(CE_Failure, CPLE_AppDefined,
2119 : "At least two sources required for interpolation.");
2120 0 : return CE_Failure;
2121 : }
2122 :
2123 17 : if (dfT == 0 || !std::isfinite(dfT))
2124 : {
2125 0 : CPLError(CE_Failure, CPLE_AppDefined, "dt must be finite and non-zero");
2126 0 : return CE_Failure;
2127 : }
2128 :
2129 17 : const auto i0 = intervalLeft(dfT0, dfDt, nSources, dfT);
2130 17 : const auto i1 = i0 + 1;
2131 17 : const double dfX0 = dfT0 + static_cast<double>(i0) * dfDt;
2132 17 : const double dfX1 = dfT0 + static_cast<double>(i0 + 1) * dfDt;
2133 :
2134 : /* ---- Set pixels ---- */
2135 17 : size_t ii = 0;
2136 41 : for (int iLine = 0; iLine < nYSize; ++iLine)
2137 : {
2138 72 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
2139 : {
2140 48 : const double dfY0 = GetSrcVal(papoSources[i0], eSrcType, ii);
2141 48 : const double dfY1 = GetSrcVal(papoSources[i1], eSrcType, ii);
2142 :
2143 48 : double dfPixVal = dfNoData;
2144 48 : if (dfT == dfX0)
2145 8 : dfPixVal = dfY0;
2146 40 : else if (dfT == dfX1)
2147 0 : dfPixVal = dfY1;
2148 52 : else if (!bHasNoData ||
2149 12 : (!IsNoData(dfY0, dfNoData) && !IsNoData(dfY1, dfNoData)))
2150 30 : dfPixVal = InterpolationFunction(dfX0, dfX1, dfY0, dfY1, dfT);
2151 :
2152 48 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
2153 : static_cast<GByte *>(pData) +
2154 48 : static_cast<GSpacing>(nLineSpace) * iLine +
2155 48 : iCol * nPixelSpace,
2156 : eBufType, nPixelSpace, 1);
2157 : }
2158 : }
2159 :
2160 : /* ---- Return success ---- */
2161 17 : return CE_None;
2162 : }
2163 :
2164 : static const char pszReplaceNoDataPixelFuncMetadata[] =
2165 : "<PixelFunctionArgumentsList>"
2166 : " <Argument type='builtin' value='NoData' />"
2167 : " <Argument name='to' type='double' description='New NoData value to be "
2168 : "replaced' default='nan' />"
2169 : "</PixelFunctionArgumentsList>";
2170 :
2171 2 : static CPLErr ReplaceNoDataPixelFunc(void **papoSources, int nSources,
2172 : void *pData, int nXSize, int nYSize,
2173 : GDALDataType eSrcType,
2174 : GDALDataType eBufType, int nPixelSpace,
2175 : int nLineSpace, CSLConstList papszArgs)
2176 : {
2177 : /* ---- Init ---- */
2178 2 : if (nSources != 1)
2179 0 : return CE_Failure;
2180 2 : if (GDALDataTypeIsComplex(eSrcType))
2181 : {
2182 0 : CPLError(CE_Failure, CPLE_AppDefined,
2183 : "replace_nodata cannot convert complex data types");
2184 0 : return CE_Failure;
2185 : }
2186 :
2187 2 : double dfOldNoData, dfNewNoData = NAN;
2188 2 : if (FetchDoubleArg(papszArgs, "NoData", &dfOldNoData) != CE_None)
2189 0 : return CE_Failure;
2190 2 : if (FetchDoubleArg(papszArgs, "to", &dfNewNoData, &dfNewNoData) != CE_None)
2191 0 : return CE_Failure;
2192 :
2193 2 : if (!GDALDataTypeIsFloating(eBufType) && std::isnan(dfNewNoData))
2194 : {
2195 0 : CPLError(CE_Failure, CPLE_AppDefined,
2196 : "Using nan requires a floating point type output buffer");
2197 0 : return CE_Failure;
2198 : }
2199 :
2200 : /* ---- Set pixels ---- */
2201 2 : size_t ii = 0;
2202 102 : for (int iLine = 0; iLine < nYSize; ++iLine)
2203 : {
2204 5100 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
2205 : {
2206 5000 : double dfPixVal = GetSrcVal(papoSources[0], eSrcType, ii);
2207 5000 : if (dfPixVal == dfOldNoData || std::isnan(dfPixVal))
2208 3200 : dfPixVal = dfNewNoData;
2209 :
2210 5000 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
2211 : static_cast<GByte *>(pData) +
2212 5000 : static_cast<GSpacing>(nLineSpace) * iLine +
2213 5000 : iCol * nPixelSpace,
2214 : eBufType, nPixelSpace, 1);
2215 : }
2216 : }
2217 :
2218 : /* ---- Return success ---- */
2219 2 : return CE_None;
2220 : }
2221 :
2222 : static const char pszScalePixelFuncMetadata[] =
2223 : "<PixelFunctionArgumentsList>"
2224 : " <Argument type='builtin' value='offset' />"
2225 : " <Argument type='builtin' value='scale' />"
2226 : " <Argument type='builtin' value='NoData' optional='true' />"
2227 : "</PixelFunctionArgumentsList>";
2228 :
2229 2 : static CPLErr ScalePixelFunc(void **papoSources, int nSources, void *pData,
2230 : int nXSize, int nYSize, GDALDataType eSrcType,
2231 : GDALDataType eBufType, int nPixelSpace,
2232 : int nLineSpace, CSLConstList papszArgs)
2233 : {
2234 : /* ---- Init ---- */
2235 2 : if (nSources != 1)
2236 0 : return CE_Failure;
2237 2 : if (GDALDataTypeIsComplex(eSrcType))
2238 : {
2239 0 : CPLError(CE_Failure, CPLE_AppDefined,
2240 : "scale cannot by applied to complex data types");
2241 0 : return CE_Failure;
2242 : }
2243 :
2244 2 : double dfNoData{0};
2245 2 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
2246 2 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
2247 0 : return CE_Failure;
2248 :
2249 : double dfScale, dfOffset;
2250 2 : if (FetchDoubleArg(papszArgs, "scale", &dfScale) != CE_None)
2251 0 : return CE_Failure;
2252 2 : if (FetchDoubleArg(papszArgs, "offset", &dfOffset) != CE_None)
2253 0 : return CE_Failure;
2254 :
2255 : /* ---- Set pixels ---- */
2256 2 : size_t ii = 0;
2257 23 : for (int iLine = 0; iLine < nYSize; ++iLine)
2258 : {
2259 423 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
2260 : {
2261 402 : const double dfVal = GetSrcVal(papoSources[0], eSrcType, ii);
2262 :
2263 2 : const double dfPixVal = bHasNoData && IsNoData(dfVal, dfNoData)
2264 404 : ? dfNoData
2265 400 : : dfVal * dfScale + dfOffset;
2266 :
2267 402 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
2268 : static_cast<GByte *>(pData) +
2269 402 : static_cast<GSpacing>(nLineSpace) * iLine +
2270 402 : iCol * nPixelSpace,
2271 : eBufType, nPixelSpace, 1);
2272 : }
2273 : }
2274 :
2275 : /* ---- Return success ---- */
2276 2 : return CE_None;
2277 : }
2278 :
2279 : static const char pszNormDiffPixelFuncMetadata[] =
2280 : "<PixelFunctionArgumentsList>"
2281 : " <Argument type='builtin' value='NoData' optional='true' />"
2282 : "</PixelFunctionArgumentsList>";
2283 :
2284 3 : static CPLErr NormDiffPixelFunc(void **papoSources, int nSources, void *pData,
2285 : int nXSize, int nYSize, GDALDataType eSrcType,
2286 : GDALDataType eBufType, int nPixelSpace,
2287 : int nLineSpace, CSLConstList papszArgs)
2288 : {
2289 : /* ---- Init ---- */
2290 3 : if (nSources != 2)
2291 0 : return CE_Failure;
2292 :
2293 3 : if (GDALDataTypeIsComplex(eSrcType))
2294 : {
2295 0 : CPLError(CE_Failure, CPLE_AppDefined,
2296 : "norm_diff cannot by applied to complex data types");
2297 0 : return CE_Failure;
2298 : }
2299 :
2300 3 : double dfNoData{0};
2301 3 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
2302 3 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
2303 0 : return CE_Failure;
2304 :
2305 : /* ---- Set pixels ---- */
2306 3 : size_t ii = 0;
2307 11 : for (int iLine = 0; iLine < nYSize; ++iLine)
2308 : {
2309 42 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
2310 : {
2311 34 : const double dfLeftVal = GetSrcVal(papoSources[0], eSrcType, ii);
2312 34 : const double dfRightVal = GetSrcVal(papoSources[1], eSrcType, ii);
2313 :
2314 34 : double dfPixVal = dfNoData;
2315 :
2316 35 : if (!bHasNoData || (!IsNoData(dfLeftVal, dfNoData) &&
2317 1 : !IsNoData(dfRightVal, dfNoData)))
2318 : {
2319 30 : const double dfDenom = (dfLeftVal + dfRightVal);
2320 : // coverity[divide_by_zero]
2321 30 : dfPixVal =
2322 : dfDenom == 0
2323 30 : ? std::numeric_limits<double>::infinity()
2324 30 : : (dfLeftVal - dfRightVal) /
2325 : #ifdef __COVERITY__
2326 : (dfDenom + std::numeric_limits<double>::min())
2327 : #else
2328 : dfDenom
2329 : #endif
2330 : ;
2331 : }
2332 :
2333 34 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
2334 : static_cast<GByte *>(pData) +
2335 34 : static_cast<GSpacing>(nLineSpace) * iLine +
2336 34 : iCol * nPixelSpace,
2337 : eBufType, nPixelSpace, 1);
2338 : }
2339 : }
2340 :
2341 : /* ---- Return success ---- */
2342 3 : return CE_None;
2343 : } // NormDiffPixelFunc
2344 :
2345 : /************************************************************************/
2346 : /* pszMinMaxFuncMetadataNodata */
2347 : /************************************************************************/
2348 :
2349 : static const char pszArgMinMaxFuncMetadataNodata[] =
2350 : "<PixelFunctionArgumentsList>"
2351 : " <Argument type='builtin' value='NoData' optional='true' />"
2352 : " <Argument name='propagateNoData' description='Whether the output value "
2353 : "should be NoData as as soon as one source is NoData' type='boolean' "
2354 : "default='false' />"
2355 : "</PixelFunctionArgumentsList>";
2356 :
2357 : static const char pszMinMaxFuncMetadataNodata[] =
2358 : "<PixelFunctionArgumentsList>"
2359 : " <Argument name='k' description='Optional constant term' type='double' "
2360 : "default='nan' />"
2361 : " <Argument type='builtin' value='NoData' optional='true' />"
2362 : " <Argument name='propagateNoData' description='Whether the output value "
2363 : "should be NoData as as soon as one source is NoData' type='boolean' "
2364 : "default='false' />"
2365 : "</PixelFunctionArgumentsList>";
2366 :
2367 : namespace
2368 : {
2369 : struct ReturnIndex;
2370 : struct ReturnValue;
2371 : } // namespace
2372 :
2373 : template <class Comparator, class ReturnType = ReturnValue>
2374 36 : static CPLErr MinOrMaxPixelFunc(double dfK, void **papoSources, int nSources,
2375 : void *pData, int nXSize, int nYSize,
2376 : GDALDataType eSrcType, GDALDataType eBufType,
2377 : int nPixelSpace, int nLineSpace,
2378 : CSLConstList papszArgs)
2379 : {
2380 : /* ---- Init ---- */
2381 36 : if (GDALDataTypeIsComplex(eSrcType))
2382 : {
2383 0 : CPLError(CE_Failure, CPLE_AppDefined,
2384 : "Complex data type not supported for min/max().");
2385 0 : return CE_Failure;
2386 : }
2387 :
2388 36 : double dfNoData = std::numeric_limits<double>::quiet_NaN();
2389 36 : if (FetchDoubleArg(papszArgs, "NoData", &dfNoData, &dfNoData) != CE_None)
2390 0 : return CE_Failure;
2391 36 : const bool bPropagateNoData = CPLTestBool(
2392 : CSLFetchNameValueDef(papszArgs, "propagateNoData", "false"));
2393 :
2394 : /* ---- Set pixels ---- */
2395 36 : size_t ii = 0;
2396 268 : for (int iLine = 0; iLine < nYSize; ++iLine)
2397 : {
2398 10290 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
2399 : {
2400 10058 : double dfRes = std::numeric_limits<double>::quiet_NaN();
2401 10058 : double dfResSrc = std::numeric_limits<double>::quiet_NaN();
2402 :
2403 33584 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
2404 : {
2405 25588 : const double dfVal = GetSrcVal(papoSources[iSrc], eSrcType, ii);
2406 :
2407 25588 : if (std::isnan(dfVal) || dfVal == dfNoData)
2408 : {
2409 14425 : if (bPropagateNoData)
2410 : {
2411 2062 : dfRes = dfNoData;
2412 : if constexpr (std::is_same_v<ReturnType, ReturnIndex>)
2413 : {
2414 4 : dfResSrc = std::numeric_limits<double>::quiet_NaN();
2415 : }
2416 2062 : break;
2417 : }
2418 : }
2419 11163 : else if (Comparator::compare(dfVal, dfRes))
2420 : {
2421 6596 : dfRes = dfVal;
2422 : if constexpr (std::is_same_v<ReturnType, ReturnIndex>)
2423 : {
2424 7 : dfResSrc = iSrc;
2425 : }
2426 : }
2427 : }
2428 :
2429 : if constexpr (std::is_same_v<ReturnType, ReturnIndex>)
2430 : {
2431 : static_cast<void>(dfK); // Placate gcc 9.4
2432 12 : dfRes = std::isnan(dfResSrc) ? dfNoData : dfResSrc + 1;
2433 : }
2434 : else
2435 : {
2436 10046 : if (std::isnan(dfRes))
2437 : {
2438 3211 : dfRes = dfNoData;
2439 : }
2440 :
2441 10046 : if (IsNoData(dfRes, dfNoData))
2442 : {
2443 5269 : if (!bPropagateNoData && !std::isnan(dfK))
2444 : {
2445 6 : dfRes = dfK;
2446 : }
2447 : }
2448 4777 : else if (!std::isnan(dfK) && Comparator::compare(dfK, dfRes))
2449 : {
2450 10 : dfRes = dfK;
2451 : }
2452 : }
2453 :
2454 10058 : GDALCopyWords(&dfRes, GDT_Float64, 0,
2455 : static_cast<GByte *>(pData) +
2456 10058 : static_cast<GSpacing>(nLineSpace) * iLine +
2457 10058 : iCol * nPixelSpace,
2458 : eBufType, nPixelSpace, 1);
2459 : }
2460 : }
2461 :
2462 : /* ---- Return success ---- */
2463 36 : return CE_None;
2464 : } /* MinOrMaxPixelFunc */
2465 :
2466 : #ifdef USE_SSE2
2467 :
2468 : template <class T, class SSEWrapper>
2469 29 : static void OptimizedMinOrMaxSSE2(const void *const *papoSources, int nSources,
2470 : void *pData, int nXSize, int nYSize,
2471 : int nLineSpace)
2472 : {
2473 29 : assert(nSources >= 1);
2474 29 : constexpr int VALUES_PER_REG =
2475 : static_cast<int>(sizeof(typename SSEWrapper::Vec) / sizeof(T));
2476 597 : for (int iLine = 0; iLine < nYSize; ++iLine)
2477 : {
2478 568 : T *CPL_RESTRICT pDest =
2479 : reinterpret_cast<T *>(static_cast<GByte *>(pData) +
2480 568 : static_cast<GSpacing>(nLineSpace) * iLine);
2481 568 : const size_t iOffsetLine = static_cast<size_t>(iLine) * nXSize;
2482 568 : int iCol = 0;
2483 3118 : for (; iCol < nXSize - (2 * VALUES_PER_REG - 1);
2484 : iCol += 2 * VALUES_PER_REG)
2485 : {
2486 2550 : auto reg0 = SSEWrapper::LoadU(
2487 2550 : static_cast<const T * CPL_RESTRICT>(papoSources[0]) +
2488 2550 : iOffsetLine + iCol);
2489 2550 : auto reg1 = SSEWrapper::LoadU(
2490 2550 : static_cast<const T * CPL_RESTRICT>(papoSources[0]) +
2491 2550 : iOffsetLine + iCol + VALUES_PER_REG);
2492 5150 : for (int iSrc = 1; iSrc < nSources; ++iSrc)
2493 : {
2494 2600 : reg0 = SSEWrapper::MinOrMax(
2495 2600 : reg0, SSEWrapper::LoadU(static_cast<const T * CPL_RESTRICT>(
2496 2600 : papoSources[iSrc]) +
2497 2600 : iOffsetLine + iCol));
2498 2600 : reg1 = SSEWrapper::MinOrMax(
2499 : reg1,
2500 : SSEWrapper::LoadU(
2501 2600 : static_cast<const T * CPL_RESTRICT>(papoSources[iSrc]) +
2502 2600 : iOffsetLine + iCol + VALUES_PER_REG));
2503 : }
2504 2550 : SSEWrapper::StoreU(pDest + iCol, reg0);
2505 2550 : SSEWrapper::StoreU(pDest + iCol + VALUES_PER_REG, reg1);
2506 : }
2507 4092 : for (; iCol < nXSize; ++iCol)
2508 : {
2509 3524 : T v = static_cast<const T * CPL_RESTRICT>(
2510 3524 : papoSources[0])[iOffsetLine + iCol];
2511 7954 : for (int iSrc = 1; iSrc < nSources; ++iSrc)
2512 : {
2513 4430 : v = SSEWrapper::MinOrMax(
2514 4430 : v, static_cast<const T * CPL_RESTRICT>(
2515 4430 : papoSources[iSrc])[iOffsetLine + iCol]);
2516 : }
2517 3524 : pDest[iCol] = v;
2518 : }
2519 : }
2520 29 : }
2521 :
2522 : // clang-format off
2523 : namespace
2524 : {
2525 : struct SSEWrapperMinByte
2526 : {
2527 : using T = uint8_t;
2528 : typedef __m128i Vec;
2529 :
2530 400 : static inline Vec LoadU(const T *x) { return _mm_loadu_si128(reinterpret_cast<const Vec*>(x)); }
2531 100 : static inline void StoreU(T *x, Vec y) { _mm_storeu_si128(reinterpret_cast<Vec*>(x), y); }
2532 200 : static inline Vec MinOrMax(Vec x, Vec y) { return _mm_min_epu8(x, y); }
2533 908 : static inline T MinOrMax(T x, T y) { return std::min(x, y); }
2534 : };
2535 :
2536 : struct SSEWrapperMaxByte
2537 : {
2538 : using T = uint8_t;
2539 : typedef __m128i Vec;
2540 :
2541 1000 : static inline Vec LoadU(const T *x) { return _mm_loadu_si128(reinterpret_cast<const Vec*>(x)); }
2542 200 : static inline void StoreU(T *x, Vec y) { _mm_storeu_si128(reinterpret_cast<Vec*>(x), y); }
2543 600 : static inline Vec MinOrMax(Vec x, Vec y) { return _mm_max_epu8(x, y); }
2544 2708 : static inline T MinOrMax(T x, T y) { return std::max(x, y); }
2545 : };
2546 :
2547 : struct SSEWrapperMinUInt16
2548 : {
2549 : using T = uint16_t;
2550 : typedef __m128i Vec;
2551 :
2552 1200 : static inline Vec LoadU(const T *x) { return _mm_loadu_si128(reinterpret_cast<const Vec*>(x)); }
2553 300 : static inline void StoreU(T *x, Vec y) { _mm_storeu_si128(reinterpret_cast<Vec*>(x), y); }
2554 : #if defined(__SSE4_1__) || defined(USE_NEON_OPTIMIZATIONS)
2555 : static inline Vec MinOrMax(Vec x, Vec y) { return _mm_min_epu16(x, y); }
2556 : #else
2557 300 : static inline Vec MinOrMax(Vec x, Vec y) { return
2558 1800 : _mm_add_epi16(
2559 : _mm_min_epi16(
2560 : _mm_add_epi16(x, _mm_set1_epi16(-32768)),
2561 : _mm_add_epi16(y, _mm_set1_epi16(-32768))),
2562 300 : _mm_set1_epi16(-32768)); }
2563 : #endif
2564 100 : static inline T MinOrMax(T x, T y) { return std::min(x, y); }
2565 : };
2566 :
2567 : struct SSEWrapperMaxUInt16
2568 : {
2569 : using T = uint16_t;
2570 : typedef __m128i Vec;
2571 :
2572 1200 : static inline Vec LoadU(const T *x) { return _mm_loadu_si128(reinterpret_cast<const Vec*>(x)); }
2573 300 : static inline void StoreU(T *x, Vec y) { _mm_storeu_si128(reinterpret_cast<Vec*>(x), y); }
2574 : #if defined(__SSE4_1__) || defined(USE_NEON_OPTIMIZATIONS)
2575 : static inline Vec MinOrMax(Vec x, Vec y) { return _mm_max_epu16(x, y); }
2576 : #else
2577 300 : static inline Vec MinOrMax(Vec x, Vec y) { return
2578 1800 : _mm_add_epi16(
2579 : _mm_max_epi16(
2580 : _mm_add_epi16(x, _mm_set1_epi16(-32768)),
2581 : _mm_add_epi16(y, _mm_set1_epi16(-32768))),
2582 300 : _mm_set1_epi16(-32768)); }
2583 : #endif
2584 100 : static inline T MinOrMax(T x, T y) { return std::max(x, y); }
2585 : };
2586 :
2587 : struct SSEWrapperMinInt16
2588 : {
2589 : using T = int16_t;
2590 : typedef __m128i Vec;
2591 :
2592 1200 : static inline Vec LoadU(const T *x) { return _mm_loadu_si128(reinterpret_cast<const Vec*>(x)); }
2593 300 : static inline void StoreU(T *x, Vec y) { _mm_storeu_si128(reinterpret_cast<Vec*>(x), y); }
2594 600 : static inline Vec MinOrMax(Vec x, Vec y) { return _mm_min_epi16(x, y); }
2595 100 : static inline T MinOrMax(T x, T y) { return std::min(x, y); }
2596 : };
2597 :
2598 : struct SSEWrapperMaxInt16
2599 : {
2600 : using T = int16_t;
2601 : typedef __m128i Vec;
2602 :
2603 1200 : static inline Vec LoadU(const T *x) { return _mm_loadu_si128(reinterpret_cast<const Vec*>(x)); }
2604 300 : static inline void StoreU(T *x, Vec y) { _mm_storeu_si128(reinterpret_cast<Vec*>(x), y); }
2605 600 : static inline Vec MinOrMax(Vec x, Vec y) { return _mm_max_epi16(x, y); }
2606 100 : static inline T MinOrMax(T x, T y) { return std::max(x, y); }
2607 : };
2608 :
2609 : struct SSEWrapperMinFloat
2610 : {
2611 : using T = float;
2612 : typedef __m128 Vec;
2613 :
2614 2400 : static inline Vec LoadU(const T *x) { return _mm_loadu_ps(x); }
2615 600 : static inline void StoreU(T *x, Vec y) { _mm_storeu_ps(x, y); }
2616 1200 : static inline Vec MinOrMax(Vec x, Vec y) { return _mm_min_ps(x, y); }
2617 100 : static inline T MinOrMax(T x, T y) { return std::min(x, y); }
2618 : };
2619 :
2620 : struct SSEWrapperMaxFloat
2621 : {
2622 : using T = float;
2623 : typedef __m128 Vec;
2624 :
2625 2400 : static inline Vec LoadU(const T *x) { return _mm_loadu_ps(x); }
2626 600 : static inline void StoreU(T *x, Vec y) { _mm_storeu_ps(x, y); }
2627 1200 : static inline Vec MinOrMax(Vec x, Vec y) { return _mm_max_ps(x, y); }
2628 100 : static inline T MinOrMax(T x, T y) { return std::max(x, y); }
2629 : };
2630 :
2631 : struct SSEWrapperMinDouble
2632 : {
2633 : using T = double;
2634 : typedef __m128d Vec;
2635 :
2636 4800 : static inline Vec LoadU(const T *x) { return _mm_loadu_pd(x); }
2637 1200 : static inline void StoreU(T *x, Vec y) { _mm_storeu_pd(x, y); }
2638 2400 : static inline Vec MinOrMax(Vec x, Vec y) { return _mm_min_pd(x, y); }
2639 107 : static inline T MinOrMax(T x, T y) { return std::min(x, y); }
2640 : };
2641 :
2642 : struct SSEWrapperMaxDouble
2643 : {
2644 : using T = double;
2645 : typedef __m128d Vec;
2646 :
2647 4800 : static inline Vec LoadU(const T *x) { return _mm_loadu_pd(x); }
2648 1200 : static inline void StoreU(T *x, Vec y) { _mm_storeu_pd(x, y); }
2649 2400 : static inline Vec MinOrMax(Vec x, Vec y) { return _mm_max_pd(x, y); }
2650 107 : static inline T MinOrMax(T x, T y) { return std::max(x, y); }
2651 : };
2652 :
2653 : } // namespace
2654 :
2655 : // clang-format on
2656 :
2657 : #endif // USE_SSE2
2658 :
2659 : template <typename ReturnType>
2660 35 : static CPLErr MinPixelFunc(void **papoSources, int nSources, void *pData,
2661 : int nXSize, int nYSize, GDALDataType eSrcType,
2662 : GDALDataType eBufType, int nPixelSpace,
2663 : int nLineSpace, CSLConstList papszArgs)
2664 : {
2665 : struct Comparator
2666 : {
2667 2746 : static bool compare(double x, double resVal)
2668 : {
2669 : // Written this way to deal with resVal being NaN
2670 2746 : return !(x >= resVal);
2671 : }
2672 : };
2673 :
2674 35 : double dfK = std::numeric_limits<double>::quiet_NaN();
2675 : if constexpr (std::is_same_v<ReturnType, ReturnValue>)
2676 : {
2677 32 : if (FetchDoubleArg(papszArgs, "k", &dfK, &dfK) != CE_None)
2678 0 : return CE_Failure;
2679 :
2680 : #ifdef USE_SSE2
2681 32 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
2682 54 : if (std::isnan(dfK) && nSources > 0 && !bHasNoData &&
2683 54 : eSrcType == eBufType &&
2684 14 : nPixelSpace == GDALGetDataTypeSizeBytes(eSrcType))
2685 : {
2686 14 : if (eSrcType == GDT_Byte)
2687 : {
2688 5 : OptimizedMinOrMaxSSE2<uint8_t, SSEWrapperMinByte>(
2689 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
2690 5 : return CE_None;
2691 : }
2692 9 : else if (eSrcType == GDT_UInt16)
2693 : {
2694 1 : OptimizedMinOrMaxSSE2<uint16_t, SSEWrapperMinUInt16>(
2695 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
2696 1 : return CE_None;
2697 : }
2698 8 : else if (eSrcType == GDT_Int16)
2699 : {
2700 1 : OptimizedMinOrMaxSSE2<int16_t, SSEWrapperMinInt16>(
2701 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
2702 1 : return CE_None;
2703 : }
2704 7 : else if (eSrcType == GDT_Float32)
2705 : {
2706 1 : OptimizedMinOrMaxSSE2<float, SSEWrapperMinFloat>(
2707 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
2708 1 : return CE_None;
2709 : }
2710 6 : else if (eSrcType == GDT_Float64)
2711 : {
2712 6 : OptimizedMinOrMaxSSE2<double, SSEWrapperMinDouble>(
2713 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
2714 6 : return CE_None;
2715 : }
2716 : }
2717 : #endif
2718 : }
2719 :
2720 21 : return MinOrMaxPixelFunc<Comparator, ReturnType>(
2721 : dfK, papoSources, nSources, pData, nXSize, nYSize, eSrcType, eBufType,
2722 24 : nPixelSpace, nLineSpace, papszArgs);
2723 : }
2724 :
2725 : template <typename ReturnType>
2726 30 : static CPLErr MaxPixelFunc(void **papoSources, int nSources, void *pData,
2727 : int nXSize, int nYSize, GDALDataType eSrcType,
2728 : GDALDataType eBufType, int nPixelSpace,
2729 : int nLineSpace, CSLConstList papszArgs)
2730 : {
2731 : struct Comparator
2732 : {
2733 8437 : static bool compare(double x, double resVal)
2734 : {
2735 : // Written this way to deal with resVal being NaN
2736 8437 : return !(x <= resVal);
2737 : }
2738 : };
2739 :
2740 30 : double dfK = std::numeric_limits<double>::quiet_NaN();
2741 : if constexpr (std::is_same_v<ReturnType, ReturnValue>)
2742 : {
2743 27 : if (FetchDoubleArg(papszArgs, "k", &dfK, &dfK) != CE_None)
2744 0 : return CE_Failure;
2745 :
2746 : #ifdef USE_SSE2
2747 27 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
2748 46 : if (std::isnan(dfK) && nSources > 0 && !bHasNoData &&
2749 46 : eSrcType == eBufType &&
2750 15 : nPixelSpace == GDALGetDataTypeSizeBytes(eSrcType))
2751 : {
2752 15 : if (eSrcType == GDT_Byte)
2753 : {
2754 6 : OptimizedMinOrMaxSSE2<uint8_t, SSEWrapperMaxByte>(
2755 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
2756 6 : return CE_None;
2757 : }
2758 9 : else if (eSrcType == GDT_UInt16)
2759 : {
2760 1 : OptimizedMinOrMaxSSE2<uint16_t, SSEWrapperMaxUInt16>(
2761 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
2762 1 : return CE_None;
2763 : }
2764 8 : else if (eSrcType == GDT_Int16)
2765 : {
2766 1 : OptimizedMinOrMaxSSE2<int16_t, SSEWrapperMaxInt16>(
2767 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
2768 1 : return CE_None;
2769 : }
2770 7 : else if (eSrcType == GDT_Float32)
2771 : {
2772 1 : OptimizedMinOrMaxSSE2<float, SSEWrapperMaxFloat>(
2773 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
2774 1 : return CE_None;
2775 : }
2776 6 : else if (eSrcType == GDT_Float64)
2777 : {
2778 6 : OptimizedMinOrMaxSSE2<double, SSEWrapperMaxDouble>(
2779 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
2780 6 : return CE_None;
2781 : }
2782 : }
2783 : #endif
2784 : }
2785 :
2786 15 : return MinOrMaxPixelFunc<Comparator, ReturnType>(
2787 : dfK, papoSources, nSources, pData, nXSize, nYSize, eSrcType, eBufType,
2788 18 : nPixelSpace, nLineSpace, papszArgs);
2789 : }
2790 :
2791 : static const char pszExprPixelFuncMetadata[] =
2792 : "<PixelFunctionArgumentsList>"
2793 : " <Argument type='builtin' value='NoData' optional='true' />"
2794 : " <Argument name='propagateNoData' description='Whether the output value "
2795 : "should be NoData as as soon as one source is NoData' type='boolean' "
2796 : "default='false' />"
2797 : " <Argument name='expression' "
2798 : " description='Expression to be evaluated' "
2799 : " type='string'></Argument>"
2800 : " <Argument name='dialect' "
2801 : " description='Expression dialect' "
2802 : " type='string-select'"
2803 : " default='muparser'>"
2804 : " <Value>exprtk</Value>"
2805 : " <Value>muparser</Value>"
2806 : " </Argument>"
2807 : " <Argument type='builtin' value='source_names' />"
2808 : " <Argument type='builtin' value='xoff' />"
2809 : " <Argument type='builtin' value='yoff' />"
2810 : " <Argument type='builtin' value='geotransform' />"
2811 : "</PixelFunctionArgumentsList>";
2812 :
2813 370 : static CPLErr ExprPixelFunc(void **papoSources, int nSources, void *pData,
2814 : int nXSize, int nYSize, GDALDataType eSrcType,
2815 : GDALDataType eBufType, int nPixelSpace,
2816 : int nLineSpace, CSLConstList papszArgs)
2817 : {
2818 : /* ---- Init ---- */
2819 370 : if (GDALDataTypeIsComplex(eSrcType))
2820 : {
2821 0 : CPLError(CE_Failure, CPLE_AppDefined,
2822 : "expression cannot by applied to complex data types");
2823 0 : return CE_Failure;
2824 : }
2825 :
2826 370 : double dfNoData{0};
2827 370 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
2828 370 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
2829 0 : return CE_Failure;
2830 :
2831 370 : const bool bPropagateNoData = CPLTestBool(
2832 : CSLFetchNameValueDef(papszArgs, "propagateNoData", "false"));
2833 :
2834 370 : const char *pszExpression = CSLFetchNameValue(papszArgs, "expression");
2835 370 : if (!pszExpression)
2836 : {
2837 1 : CPLError(CE_Failure, CPLE_AppDefined,
2838 : "Missing 'expression' pixel function argument");
2839 1 : return CE_Failure;
2840 : }
2841 :
2842 369 : const char *pszSourceNames = CSLFetchNameValue(papszArgs, "source_names");
2843 : const CPLStringList aosSourceNames(
2844 738 : CSLTokenizeString2(pszSourceNames, "|", 0));
2845 369 : if (aosSourceNames.size() != nSources)
2846 : {
2847 7 : CPLError(CE_Failure, CPLE_AppDefined,
2848 : "The source_names variable passed to ExprPixelFunc() has %d "
2849 : "values, whereas %d were expected. An invalid variable name "
2850 : "has likely been used",
2851 : aosSourceNames.size(), nSources);
2852 7 : return CE_Failure;
2853 : }
2854 :
2855 724 : std::vector<double> adfValuesForPixel(nSources);
2856 :
2857 362 : const char *pszDialect = CSLFetchNameValue(papszArgs, "dialect");
2858 362 : if (!pszDialect)
2859 : {
2860 230 : pszDialect = "muparser";
2861 : }
2862 :
2863 724 : auto poExpression = gdal::MathExpression::Create(pszExpression, pszDialect);
2864 :
2865 : // cppcheck-suppress knownConditionTrueFalse
2866 362 : if (!poExpression)
2867 : {
2868 0 : return CE_Failure;
2869 : }
2870 :
2871 362 : int nXOff = 0;
2872 362 : int nYOff = 0;
2873 362 : GDALGeoTransform gt;
2874 362 : double dfCenterX = 0;
2875 362 : double dfCenterY = 0;
2876 :
2877 362 : bool includeCenterCoords = false;
2878 362 : if (strstr(pszExpression, "_CENTER_X_") ||
2879 360 : strstr(pszExpression, "_CENTER_Y_"))
2880 : {
2881 2 : includeCenterCoords = true;
2882 :
2883 2 : const char *pszXOff = CSLFetchNameValue(papszArgs, "xoff");
2884 2 : nXOff = std::atoi(pszXOff);
2885 :
2886 2 : const char *pszYOff = CSLFetchNameValue(papszArgs, "yoff");
2887 2 : nYOff = std::atoi(pszYOff);
2888 :
2889 2 : const char *pszGT = CSLFetchNameValue(papszArgs, "geotransform");
2890 2 : if (pszGT == nullptr)
2891 : {
2892 1 : CPLError(CE_Failure, CPLE_AppDefined,
2893 : "To use _CENTER_X_ or _CENTER_Y_ in an expression, "
2894 : "VRTDataset must have a <GeoTransform> element.");
2895 1 : return CE_Failure;
2896 : }
2897 :
2898 : CPLStringList aosGeoTransform(
2899 1 : CSLTokenizeString2(pszGT, ",", CSLT_HONOURSTRINGS));
2900 1 : if (aosGeoTransform.size() != 6)
2901 : {
2902 0 : CPLError(CE_Failure, CPLE_AppDefined,
2903 : "Invalid GeoTransform argument");
2904 0 : return CE_Failure;
2905 : }
2906 7 : for (int i = 0; i < 6; i++)
2907 : {
2908 6 : gt[i] = CPLAtof(aosGeoTransform[i]);
2909 : }
2910 : }
2911 :
2912 : {
2913 361 : int iSource = 0;
2914 966 : for (const auto &osName : aosSourceNames)
2915 : {
2916 1210 : poExpression->RegisterVariable(osName,
2917 605 : &adfValuesForPixel[iSource++]);
2918 : }
2919 : }
2920 :
2921 361 : if (includeCenterCoords)
2922 : {
2923 1 : poExpression->RegisterVariable("_CENTER_X_", &dfCenterX);
2924 1 : poExpression->RegisterVariable("_CENTER_Y_", &dfCenterY);
2925 : }
2926 :
2927 361 : if (bHasNoData)
2928 : {
2929 8 : poExpression->RegisterVariable("NODATA", &dfNoData);
2930 : }
2931 :
2932 361 : if (strstr(pszExpression, "BANDS"))
2933 : {
2934 2 : poExpression->RegisterVector("BANDS", &adfValuesForPixel);
2935 : }
2936 :
2937 : std::unique_ptr<double, VSIFreeReleaser> padfResults(
2938 722 : static_cast<double *>(VSI_MALLOC2_VERBOSE(nXSize, sizeof(double))));
2939 361 : if (!padfResults)
2940 0 : return CE_Failure;
2941 :
2942 : /* ---- Set pixels ---- */
2943 361 : size_t ii = 0;
2944 5097 : for (int iLine = 0; iLine < nYSize; ++iLine)
2945 : {
2946 12996300 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
2947 : {
2948 12991500 : double &dfResult = padfResults.get()[iCol];
2949 12991500 : bool resultIsNoData = false;
2950 :
2951 38978400 : for (int iSrc = 0; iSrc < nSources; iSrc++)
2952 : {
2953 : // cppcheck-suppress unreadVariable
2954 25986900 : double dfVal = GetSrcVal(papoSources[iSrc], eSrcType, ii);
2955 :
2956 25986900 : if (bHasNoData && bPropagateNoData && IsNoData(dfVal, dfNoData))
2957 : {
2958 1 : resultIsNoData = true;
2959 : }
2960 :
2961 25986900 : adfValuesForPixel[iSrc] = dfVal;
2962 : }
2963 :
2964 12991500 : if (includeCenterCoords)
2965 : {
2966 : // Add 0.5 to pixel / line to move from pixel corner to cell center
2967 400 : gt.Apply(static_cast<double>(iCol + nXOff) + 0.5,
2968 400 : static_cast<double>(iLine + nYOff) + 0.5, &dfCenterX,
2969 : &dfCenterY);
2970 : }
2971 :
2972 12991500 : if (resultIsNoData)
2973 : {
2974 1 : dfResult = dfNoData;
2975 : }
2976 : else
2977 : {
2978 12991500 : if (auto eErr = poExpression->Evaluate(); eErr != CE_None)
2979 : {
2980 5 : return CE_Failure;
2981 : }
2982 :
2983 12991500 : dfResult = poExpression->Results()[0];
2984 : }
2985 : }
2986 :
2987 4736 : GDALCopyWords(padfResults.get(), GDT_Float64, sizeof(double),
2988 : static_cast<GByte *>(pData) +
2989 4736 : static_cast<GSpacing>(nLineSpace) * iLine,
2990 : eBufType, nPixelSpace, nXSize);
2991 : }
2992 :
2993 : /* ---- Return success ---- */
2994 356 : return CE_None;
2995 : } // ExprPixelFunc
2996 :
2997 : static const char pszReclassifyPixelFuncMetadata[] =
2998 : "<PixelFunctionArgumentsList>"
2999 : " <Argument name='mapping' "
3000 : " description='Lookup table for mapping, in format "
3001 : "from=to,from=to' "
3002 : " type='string'></Argument>"
3003 : " <Argument type='builtin' value='NoData' optional='true' />"
3004 : "</PixelFunctionArgumentsList>";
3005 :
3006 33 : static CPLErr ReclassifyPixelFunc(void **papoSources, int nSources, void *pData,
3007 : int nXSize, int nYSize, GDALDataType eSrcType,
3008 : GDALDataType eBufType, int nPixelSpace,
3009 : int nLineSpace, CSLConstList papszArgs)
3010 : {
3011 33 : if (GDALDataTypeIsComplex(eSrcType))
3012 : {
3013 0 : CPLError(CE_Failure, CPLE_AppDefined,
3014 : "reclassify cannot by applied to complex data types");
3015 0 : return CE_Failure;
3016 : }
3017 :
3018 33 : if (nSources != 1)
3019 : {
3020 0 : CPLError(CE_Failure, CPLE_AppDefined,
3021 : "reclassify only be applied to a single source at a time");
3022 0 : return CE_Failure;
3023 : }
3024 33 : std::optional<double> noDataValue{};
3025 :
3026 33 : const char *pszNoData = CSLFetchNameValue(papszArgs, "NoData");
3027 33 : if (pszNoData != nullptr)
3028 : {
3029 10 : noDataValue = CPLAtof(pszNoData);
3030 : }
3031 :
3032 33 : const char *pszMappings = CSLFetchNameValue(papszArgs, "mapping");
3033 33 : if (pszMappings == nullptr)
3034 : {
3035 0 : CPLError(CE_Failure, CPLE_AppDefined,
3036 : "reclassify must be called with 'mapping' argument");
3037 0 : return CE_Failure;
3038 : }
3039 :
3040 66 : gdal::Reclassifier oReclassifier;
3041 33 : if (auto eErr = oReclassifier.Init(pszMappings, noDataValue, eBufType);
3042 : eErr != CE_None)
3043 : {
3044 14 : return eErr;
3045 : }
3046 :
3047 : std::unique_ptr<double, VSIFreeReleaser> padfResults(
3048 38 : static_cast<double *>(VSI_MALLOC2_VERBOSE(nXSize, sizeof(double))));
3049 19 : if (!padfResults)
3050 0 : return CE_Failure;
3051 :
3052 19 : size_t ii = 0;
3053 19 : bool bSuccess = false;
3054 435 : for (int iLine = 0; iLine < nYSize; ++iLine)
3055 : {
3056 20805 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
3057 : {
3058 20389 : double srcVal = GetSrcVal(papoSources[0], eSrcType, ii);
3059 40778 : padfResults.get()[iCol] =
3060 20389 : oReclassifier.Reclassify(srcVal, bSuccess);
3061 20389 : if (!bSuccess)
3062 : {
3063 2 : CPLError(CE_Failure, CPLE_AppDefined,
3064 : "Encountered value %g with no specified mapping",
3065 : srcVal);
3066 2 : return CE_Failure;
3067 : }
3068 : }
3069 :
3070 416 : GDALCopyWords(padfResults.get(), GDT_Float64, sizeof(double),
3071 : static_cast<GByte *>(pData) +
3072 416 : static_cast<GSpacing>(nLineSpace) * iLine,
3073 : eBufType, nPixelSpace, nXSize);
3074 : }
3075 :
3076 17 : return CE_None;
3077 : } // ReclassifyPixelFunc
3078 :
3079 : struct MeanKernel
3080 : {
3081 : static constexpr const char *pszName = "mean";
3082 :
3083 : double dfMean = 0;
3084 : int nValidSources = 0;
3085 :
3086 3466 : void Reset()
3087 : {
3088 3466 : dfMean = 0;
3089 3466 : nValidSources = 0;
3090 3466 : }
3091 :
3092 93 : static CPLErr ProcessArguments(CSLConstList)
3093 : {
3094 93 : return CE_None;
3095 : }
3096 :
3097 9271 : void ProcessPixel(double dfVal)
3098 : {
3099 9271 : ++nValidSources;
3100 :
3101 9271 : if (CPL_UNLIKELY(std::isinf(dfVal)))
3102 : {
3103 620 : if (nValidSources == 1)
3104 : {
3105 310 : dfMean = dfVal;
3106 : }
3107 310 : else if (dfVal == -dfMean)
3108 : {
3109 62 : dfMean = std::numeric_limits<double>::quiet_NaN();
3110 : }
3111 : }
3112 8651 : else if (CPL_UNLIKELY(std::isinf(dfMean)))
3113 : {
3114 186 : if (!std::isfinite(dfVal))
3115 : {
3116 62 : dfMean = std::numeric_limits<double>::quiet_NaN();
3117 : }
3118 : }
3119 : else
3120 : {
3121 8465 : const double delta = dfVal - dfMean;
3122 8465 : if (CPL_UNLIKELY(std::isinf(delta)))
3123 0 : dfMean += dfVal / nValidSources - dfMean / nValidSources;
3124 : else
3125 8465 : dfMean += delta / nValidSources;
3126 : }
3127 9271 : }
3128 :
3129 3464 : bool HasValue() const
3130 : {
3131 3464 : return nValidSources > 0;
3132 : }
3133 :
3134 3461 : double GetValue() const
3135 : {
3136 3461 : return dfMean;
3137 : }
3138 : };
3139 :
3140 : struct GeoMeanKernel
3141 : {
3142 : static constexpr const char *pszName = "geometric_mean";
3143 :
3144 : double dfProduct = 1;
3145 : int nValidSources = 0;
3146 :
3147 6 : void Reset()
3148 : {
3149 6 : dfProduct = 1;
3150 6 : nValidSources = 0;
3151 6 : }
3152 :
3153 3 : static CPLErr ProcessArguments(CSLConstList)
3154 : {
3155 3 : return CE_None;
3156 : }
3157 :
3158 3 : void ProcessPixel(double dfVal)
3159 : {
3160 3 : dfProduct *= dfVal;
3161 3 : nValidSources++;
3162 3 : }
3163 :
3164 4 : bool HasValue() const
3165 : {
3166 4 : return nValidSources > 0;
3167 : }
3168 :
3169 1 : double GetValue() const
3170 : {
3171 1 : return std::pow(dfProduct, 1.0 / nValidSources);
3172 : }
3173 : };
3174 :
3175 : struct HarmonicMeanKernel
3176 : {
3177 : static constexpr const char *pszName = "harmonic_mean";
3178 :
3179 : double dfDenom = 0;
3180 : int nValidSources = 0;
3181 : bool bValueIsZero = false;
3182 : bool bPropagateZero = false;
3183 :
3184 10 : void Reset()
3185 : {
3186 10 : dfDenom = 0;
3187 10 : nValidSources = 0;
3188 10 : bValueIsZero = false;
3189 10 : }
3190 :
3191 7 : void ProcessPixel(double dfVal)
3192 : {
3193 7 : if (dfVal == 0)
3194 : {
3195 2 : bValueIsZero = true;
3196 : }
3197 : else
3198 : {
3199 5 : dfDenom += 1 / dfVal;
3200 : }
3201 7 : nValidSources++;
3202 7 : }
3203 :
3204 5 : CPLErr ProcessArguments(CSLConstList papszArgs)
3205 : {
3206 5 : bPropagateZero =
3207 5 : CPLTestBool(CSLFetchNameValueDef(papszArgs, "propagateZero", "0"));
3208 5 : return CE_None;
3209 : }
3210 :
3211 8 : bool HasValue() const
3212 : {
3213 8 : return dfDenom > 0 && (bPropagateZero || !bValueIsZero);
3214 : }
3215 :
3216 2 : double GetValue() const
3217 : {
3218 2 : if (bPropagateZero && bValueIsZero)
3219 : {
3220 1 : return 0;
3221 : }
3222 1 : return static_cast<double>(nValidSources) / dfDenom;
3223 : }
3224 : };
3225 :
3226 : struct MedianKernel
3227 : {
3228 : static constexpr const char *pszName = "median";
3229 :
3230 : mutable std::vector<double> values{};
3231 :
3232 9 : void Reset()
3233 : {
3234 9 : values.clear();
3235 9 : }
3236 :
3237 5 : static CPLErr ProcessArguments(CSLConstList)
3238 : {
3239 5 : return CE_None;
3240 : }
3241 :
3242 9 : void ProcessPixel(double dfVal)
3243 : {
3244 9 : if (!std::isnan(dfVal))
3245 : {
3246 9 : values.push_back(dfVal);
3247 : }
3248 9 : }
3249 :
3250 7 : bool HasValue() const
3251 : {
3252 7 : return !values.empty();
3253 : }
3254 :
3255 3 : double GetValue() const
3256 : {
3257 3 : std::sort(values.begin(), values.end());
3258 3 : if (values.size() % 2 == 0)
3259 : {
3260 : return 0.5 *
3261 1 : (values[values.size() / 2 - 1] + values[values.size() / 2]);
3262 : }
3263 :
3264 2 : return values[values.size() / 2];
3265 : }
3266 : };
3267 :
3268 : struct ModeKernel
3269 : {
3270 : static constexpr const char *pszName = "mode";
3271 :
3272 : std::map<double, size_t> counts{};
3273 : std::size_t nanCount{0};
3274 : double dfMax = std::numeric_limits<double>::quiet_NaN();
3275 : decltype(counts.begin()) oMax = counts.end();
3276 :
3277 7 : void Reset()
3278 : {
3279 7 : nanCount = 0;
3280 7 : counts.clear();
3281 7 : oMax = counts.end();
3282 7 : }
3283 :
3284 4 : static CPLErr ProcessArguments(CSLConstList)
3285 : {
3286 4 : return CE_None;
3287 : }
3288 :
3289 11 : void ProcessPixel(double dfVal)
3290 : {
3291 11 : if (std::isnan(dfVal))
3292 : {
3293 2 : nanCount += 1;
3294 2 : return;
3295 : }
3296 :
3297 : // if dfVal is NaN, try_emplace will return an entry for a different key!
3298 9 : auto [it, inserted] = counts.try_emplace(dfVal, 0);
3299 :
3300 9 : it->second += 1;
3301 :
3302 9 : if (oMax == counts.end() || it->second > oMax->second)
3303 : {
3304 5 : oMax = it;
3305 : }
3306 : }
3307 :
3308 5 : bool HasValue() const
3309 : {
3310 5 : return nanCount > 0 || oMax != counts.end();
3311 : }
3312 :
3313 3 : double GetValue() const
3314 : {
3315 3 : double ret = std::numeric_limits<double>::quiet_NaN();
3316 3 : if (oMax != counts.end())
3317 : {
3318 3 : const size_t nCount = oMax->second;
3319 3 : if (nCount > nanCount)
3320 2 : ret = oMax->first;
3321 : }
3322 3 : return ret;
3323 : }
3324 : };
3325 :
3326 : static const char pszBasicPixelFuncMetadata[] =
3327 : "<PixelFunctionArgumentsList>"
3328 : " <Argument type='builtin' value='NoData' optional='true' />"
3329 : " <Argument name='propagateNoData' description='Whether the output value "
3330 : "should be NoData as as soon as one source is NoData' type='boolean' "
3331 : "default='false' />"
3332 : "</PixelFunctionArgumentsList>";
3333 :
3334 : #if defined(USE_SSE2) && !defined(USE_NEON_OPTIMIZATIONS)
3335 636 : inline __m128i packus_epi32(__m128i low, __m128i high)
3336 : {
3337 : #if __SSE4_1__
3338 : return _mm_packus_epi32(low, high); // Pack uint32 to uint16
3339 : #else
3340 1272 : low = _mm_add_epi32(low, _mm_set1_epi32(-32768));
3341 1272 : high = _mm_add_epi32(high, _mm_set1_epi32(-32768));
3342 1908 : return _mm_sub_epi16(_mm_packs_epi32(low, high), _mm_set1_epi16(-32768));
3343 : #endif
3344 : }
3345 : #endif
3346 :
3347 : #ifdef USE_SSE2
3348 :
3349 : template <class T, class SSEWrapper>
3350 41 : static void OptimizedMeanFloatSSE2(const void *const *papoSources, int nSources,
3351 : void *pData, int nXSize, int nYSize,
3352 : int nLineSpace)
3353 : {
3354 41 : assert(nSources >= 1);
3355 41 : constexpr int VALUES_PER_REG =
3356 : static_cast<int>(sizeof(typename SSEWrapper::Vec) / sizeof(T));
3357 41 : const T invSources = static_cast<T>(1.0) / static_cast<T>(nSources);
3358 41 : const auto invSourcesSSE = SSEWrapper::Set1(invSources);
3359 41 : const auto signMaskSSE = SSEWrapper::Set1(static_cast<T>(-0.0));
3360 41 : const auto infSSE = SSEWrapper::Set1(std::numeric_limits<T>::infinity());
3361 180 : for (int iLine = 0; iLine < nYSize; ++iLine)
3362 : {
3363 139 : T *CPL_RESTRICT pDest =
3364 : reinterpret_cast<T *>(static_cast<GByte *>(pData) +
3365 139 : static_cast<GSpacing>(nLineSpace) * iLine);
3366 139 : const size_t iOffsetLine = static_cast<size_t>(iLine) * nXSize;
3367 139 : int iCol = 0;
3368 1149 : for (; iCol < nXSize - (2 * VALUES_PER_REG - 1);
3369 : iCol += 2 * VALUES_PER_REG)
3370 : {
3371 1025 : auto reg0 = SSEWrapper::LoadU(
3372 1025 : static_cast<const T * CPL_RESTRICT>(papoSources[0]) +
3373 1025 : iOffsetLine + iCol);
3374 1025 : auto reg1 = SSEWrapper::LoadU(
3375 1025 : static_cast<const T * CPL_RESTRICT>(papoSources[0]) +
3376 1025 : iOffsetLine + iCol + VALUES_PER_REG);
3377 3000 : for (int iSrc = 1; iSrc < nSources; ++iSrc)
3378 : {
3379 1975 : const auto inputVal0 = SSEWrapper::LoadU(
3380 1975 : static_cast<const T * CPL_RESTRICT>(papoSources[iSrc]) +
3381 1975 : iOffsetLine + iCol);
3382 1975 : const auto inputVal1 = SSEWrapper::LoadU(
3383 1975 : static_cast<const T * CPL_RESTRICT>(papoSources[iSrc]) +
3384 1975 : iOffsetLine + iCol + VALUES_PER_REG);
3385 1975 : reg0 = SSEWrapper::Add(reg0, inputVal0);
3386 1975 : reg1 = SSEWrapper::Add(reg1, inputVal1);
3387 : }
3388 1025 : reg0 = SSEWrapper::Mul(reg0, invSourcesSSE);
3389 1025 : reg1 = SSEWrapper::Mul(reg1, invSourcesSSE);
3390 :
3391 : // Detect infinity that could happen when summing huge
3392 : // values
3393 1025 : if (SSEWrapper::MoveMask(SSEWrapper::Or(
3394 : SSEWrapper::CmpEq(SSEWrapper::AndNot(signMaskSSE, reg0),
3395 : infSSE),
3396 : SSEWrapper::CmpEq(SSEWrapper::AndNot(signMaskSSE, reg1),
3397 : infSSE))))
3398 : {
3399 15 : break;
3400 : }
3401 :
3402 1010 : SSEWrapper::StoreU(pDest + iCol, reg0);
3403 1010 : SSEWrapper::StoreU(pDest + iCol + VALUES_PER_REG, reg1);
3404 : }
3405 :
3406 : // Use numerically stable mean computation
3407 916 : for (; iCol < nXSize; ++iCol)
3408 : {
3409 777 : T mean = static_cast<const T * CPL_RESTRICT>(
3410 777 : papoSources[0])[iOffsetLine + iCol];
3411 777 : if (nSources >= 2)
3412 : {
3413 767 : T new_val = static_cast<const T * CPL_RESTRICT>(
3414 767 : papoSources[1])[iOffsetLine + iCol];
3415 767 : if (CPL_UNLIKELY(std::isinf(new_val)))
3416 : {
3417 268 : if (new_val == -mean)
3418 : {
3419 10 : pDest[iCol] = std::numeric_limits<T>::quiet_NaN();
3420 10 : continue;
3421 : }
3422 : }
3423 499 : else if (CPL_UNLIKELY(std::isinf(mean)))
3424 : {
3425 144 : if (!std::isfinite(new_val))
3426 : {
3427 10 : pDest[iCol] = std::numeric_limits<T>::quiet_NaN();
3428 10 : continue;
3429 : }
3430 : }
3431 : else
3432 : {
3433 355 : const T delta = new_val - mean;
3434 355 : if (CPL_UNLIKELY(std::isinf(delta)))
3435 10 : mean += new_val * static_cast<T>(0.5) -
3436 10 : mean * static_cast<T>(0.5);
3437 : else
3438 345 : mean += delta * static_cast<T>(0.5);
3439 : }
3440 :
3441 1279 : for (int iSrc = 2; iSrc < nSources; ++iSrc)
3442 : {
3443 552 : new_val = static_cast<const T * CPL_RESTRICT>(
3444 552 : papoSources[iSrc])[iOffsetLine + iCol];
3445 552 : if (CPL_UNLIKELY(std::isinf(new_val)))
3446 : {
3447 196 : if (new_val == -mean)
3448 : {
3449 10 : mean = std::numeric_limits<T>::quiet_NaN();
3450 10 : break;
3451 : }
3452 : }
3453 356 : else if (CPL_UNLIKELY(std::isinf(mean)))
3454 : {
3455 72 : if (!std::isfinite(new_val))
3456 : {
3457 10 : mean = std::numeric_limits<T>::quiet_NaN();
3458 10 : break;
3459 : }
3460 : }
3461 : else
3462 : {
3463 284 : const T delta = new_val - mean;
3464 284 : if (CPL_UNLIKELY(std::isinf(delta)))
3465 62 : mean += new_val / static_cast<T>(iSrc + 1) -
3466 62 : mean / static_cast<T>(iSrc + 1);
3467 : else
3468 222 : mean += delta / static_cast<T>(iSrc + 1);
3469 : }
3470 : }
3471 : }
3472 757 : pDest[iCol] = mean;
3473 : }
3474 : }
3475 41 : }
3476 :
3477 : // clang-format off
3478 : namespace
3479 : {
3480 : #ifdef __AVX2__
3481 : struct SSEWrapperFloat
3482 : {
3483 : typedef __m256 Vec;
3484 :
3485 : static inline Vec Set1(float x) { return _mm256_set1_ps(x); }
3486 : static inline Vec LoadU(const float *x) { return _mm256_loadu_ps(x); }
3487 : static inline void StoreU(float *x, Vec y) { _mm256_storeu_ps(x, y); }
3488 : static inline Vec Add(Vec x, Vec y) { return _mm256_add_ps(x, y); }
3489 : static inline Vec Mul(Vec x, Vec y) { return _mm256_mul_ps(x, y); }
3490 : static inline Vec Or(Vec x, Vec y) { return _mm256_or_ps(x, y); }
3491 : static inline Vec AndNot(Vec x, Vec y) { return _mm256_andnot_ps(x, y); }
3492 : static inline Vec CmpEq(Vec x, Vec y) { return _mm256_cmp_ps(x, y, _CMP_EQ_OQ); }
3493 : static inline int MoveMask(Vec x) { return _mm256_movemask_ps(x); }
3494 : };
3495 :
3496 : struct SSEWrapperDouble
3497 : {
3498 : typedef __m256d Vec;
3499 :
3500 : static inline Vec Set1(double x) { return _mm256_set1_pd(x); }
3501 : static inline Vec LoadU(const double *x) { return _mm256_loadu_pd(x); }
3502 : static inline void StoreU(double *x, Vec y) { _mm256_storeu_pd(x, y); }
3503 : static inline Vec Add(Vec x, Vec y) { return _mm256_add_pd(x, y); }
3504 : static inline Vec Mul(Vec x, Vec y) { return _mm256_mul_pd(x, y); }
3505 : static inline Vec Or(Vec x, Vec y) { return _mm256_or_pd(x, y); }
3506 : static inline Vec AndNot(Vec x, Vec y) { return _mm256_andnot_pd(x, y); }
3507 : static inline Vec CmpEq(Vec x, Vec y) { return _mm256_cmp_pd(x, y, _CMP_EQ_OQ); }
3508 : static inline int MoveMask(Vec x) { return _mm256_movemask_pd(x); }
3509 : };
3510 :
3511 : #else
3512 :
3513 : struct SSEWrapperFloat
3514 : {
3515 : typedef __m128 Vec;
3516 :
3517 114 : static inline Vec Set1(float x) { return _mm_set1_ps(x); }
3518 3984 : static inline Vec LoadU(const float *x) { return _mm_loadu_ps(x); }
3519 666 : static inline void StoreU(float *x, Vec y) { _mm_storeu_ps(x, y); }
3520 2624 : static inline Vec Add(Vec x, Vec y) { return _mm_add_ps(x, y); }
3521 1360 : static inline Vec Mul(Vec x, Vec y) { return _mm_mul_ps(x, y); }
3522 680 : static inline Vec Or(Vec x, Vec y) { return _mm_or_ps(x, y); }
3523 1360 : static inline Vec AndNot(Vec x, Vec y) { return _mm_andnot_ps(x, y); }
3524 1360 : static inline Vec CmpEq(Vec x, Vec y) { return _mm_cmpeq_ps(x, y); }
3525 680 : static inline int MoveMask(Vec x) { return _mm_movemask_ps(x); }
3526 : };
3527 :
3528 : struct SSEWrapperDouble
3529 : {
3530 : typedef __m128d Vec;
3531 :
3532 132 : static inline Vec Set1(double x) { return _mm_set1_pd(x); }
3533 8016 : static inline Vec LoadU(const double *x) { return _mm_loadu_pd(x); }
3534 1354 : static inline void StoreU(double *x, Vec y) { _mm_storeu_pd(x, y); }
3535 5276 : static inline Vec Add(Vec x, Vec y) { return _mm_add_pd(x, y); }
3536 2740 : static inline Vec Mul(Vec x, Vec y) { return _mm_mul_pd(x, y); }
3537 1370 : static inline Vec Or(Vec x, Vec y) { return _mm_or_pd(x, y); }
3538 2740 : static inline Vec AndNot(Vec x, Vec y) { return _mm_andnot_pd(x, y); }
3539 2740 : static inline Vec CmpEq(Vec x, Vec y) { return _mm_cmpeq_pd(x, y); }
3540 1370 : static inline int MoveMask(Vec x) { return _mm_movemask_pd(x); }
3541 : };
3542 : #endif
3543 : } // namespace
3544 :
3545 : // clang-format on
3546 :
3547 : #endif // USE_SSE2
3548 :
3549 : template <typename Kernel>
3550 110 : static CPLErr BasicPixelFunc(void **papoSources, int nSources, void *pData,
3551 : int nXSize, int nYSize, GDALDataType eSrcType,
3552 : GDALDataType eBufType, int nPixelSpace,
3553 : int nLineSpace, CSLConstList papszArgs)
3554 : {
3555 : /* ---- Init ---- */
3556 119 : Kernel oKernel;
3557 :
3558 110 : if (GDALDataTypeIsComplex(eSrcType))
3559 : {
3560 0 : CPLError(CE_Failure, CPLE_AppDefined,
3561 : "Complex data types not supported by %s", oKernel.pszName);
3562 0 : return CE_Failure;
3563 : }
3564 :
3565 110 : double dfNoData{0};
3566 110 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
3567 110 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
3568 0 : return CE_Failure;
3569 :
3570 110 : const bool bPropagateNoData = CPLTestBool(
3571 : CSLFetchNameValueDef(papszArgs, "propagateNoData", "false"));
3572 :
3573 110 : if (oKernel.ProcessArguments(papszArgs) == CE_Failure)
3574 : {
3575 0 : return CE_Failure;
3576 : }
3577 :
3578 : #if defined(USE_SSE2) && !defined(USE_NEON_OPTIMIZATIONS)
3579 : if constexpr (std::is_same_v<Kernel, MeanKernel>)
3580 : {
3581 90 : if (!bHasNoData && eSrcType == GDT_Byte && eBufType == GDT_Byte &&
3582 183 : nPixelSpace == 1 &&
3583 : // We use signed int16 to accumulate
3584 15 : nSources <= std::numeric_limits<int16_t>::max() /
3585 15 : std::numeric_limits<uint8_t>::max())
3586 : {
3587 : using T = uint8_t;
3588 14 : constexpr int VALUES_PER_REG = 16;
3589 14 : if (nSources == 2)
3590 : {
3591 211 : for (int iLine = 0; iLine < nYSize; ++iLine)
3592 : {
3593 205 : T *CPL_RESTRICT pDest = reinterpret_cast<T *>(
3594 : static_cast<GByte *>(pData) +
3595 205 : static_cast<GSpacing>(nLineSpace) * iLine);
3596 205 : const size_t iOffsetLine =
3597 205 : static_cast<size_t>(iLine) * nXSize;
3598 205 : int iCol = 0;
3599 5211 : for (; iCol < nXSize - (VALUES_PER_REG - 1);
3600 : iCol += VALUES_PER_REG)
3601 : {
3602 : const __m128i inputVal0 =
3603 10012 : _mm_loadu_si128(reinterpret_cast<const __m128i *>(
3604 : static_cast<const T * CPL_RESTRICT>(
3605 : papoSources[0]) +
3606 5006 : iOffsetLine + iCol));
3607 : const __m128i inputVal1 =
3608 5006 : _mm_loadu_si128(reinterpret_cast<const __m128i *>(
3609 5006 : static_cast<const T * CPL_RESTRICT>(
3610 : papoSources[1]) +
3611 5006 : iOffsetLine + iCol));
3612 5006 : _mm_storeu_si128(
3613 5006 : reinterpret_cast<__m128i *>(pDest + iCol),
3614 : _mm_avg_epu8(inputVal0, inputVal1));
3615 : }
3616 239 : for (; iCol < nXSize; ++iCol)
3617 : {
3618 34 : uint32_t acc = 1 +
3619 34 : static_cast<const T * CPL_RESTRICT>(
3620 34 : papoSources[0])[iOffsetLine + iCol] +
3621 34 : static_cast<const T * CPL_RESTRICT>(
3622 34 : papoSources[1])[iOffsetLine + iCol];
3623 34 : pDest[iCol] = static_cast<T>(acc / 2);
3624 : }
3625 : }
3626 : }
3627 : else
3628 : {
3629 8 : libdivide::divider<uint16_t> fast_d(
3630 : static_cast<uint16_t>(nSources));
3631 : const auto halfConstant =
3632 8 : _mm_set1_epi16(static_cast<int16_t>(nSources / 2));
3633 215 : for (int iLine = 0; iLine < nYSize; ++iLine)
3634 : {
3635 207 : T *CPL_RESTRICT pDest =
3636 : static_cast<GByte *>(pData) +
3637 207 : static_cast<GSpacing>(nLineSpace) * iLine;
3638 207 : const size_t iOffsetLine =
3639 207 : static_cast<size_t>(iLine) * nXSize;
3640 207 : int iCol = 0;
3641 5222 : for (; iCol < nXSize - (VALUES_PER_REG - 1);
3642 : iCol += VALUES_PER_REG)
3643 : {
3644 5015 : __m128i reg0 = halfConstant;
3645 5015 : __m128i reg1 = halfConstant;
3646 20435 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
3647 : {
3648 15420 : const __m128i inputVal = _mm_loadu_si128(
3649 : reinterpret_cast<const __m128i *>(
3650 15420 : static_cast<const T * CPL_RESTRICT>(
3651 15420 : papoSources[iSrc]) +
3652 15420 : iOffsetLine + iCol));
3653 46260 : reg0 = _mm_add_epi16(
3654 : reg0, _mm_unpacklo_epi8(inputVal,
3655 : _mm_setzero_si128()));
3656 46260 : reg1 = _mm_add_epi16(
3657 : reg1, _mm_unpackhi_epi8(inputVal,
3658 : _mm_setzero_si128()));
3659 : }
3660 : reg0 /= fast_d;
3661 : reg1 /= fast_d;
3662 5015 : _mm_storeu_si128(
3663 5015 : reinterpret_cast<__m128i *>(pDest + iCol),
3664 : _mm_packus_epi16(reg0, reg1));
3665 : }
3666 284 : for (; iCol < nXSize; ++iCol)
3667 : {
3668 77 : uint32_t acc = nSources / 2;
3669 2181 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
3670 : {
3671 2104 : acc += static_cast<const T * CPL_RESTRICT>(
3672 2104 : papoSources[iSrc])[iOffsetLine + iCol];
3673 : }
3674 77 : pDest[iCol] = static_cast<T>(acc / nSources);
3675 : }
3676 : }
3677 : }
3678 14 : return CE_None;
3679 : }
3680 :
3681 76 : if (!bHasNoData && eSrcType == GDT_Byte && eBufType == GDT_Byte &&
3682 155 : nPixelSpace == 1 &&
3683 : // We use signed int32 to accumulate
3684 1 : nSources <= std::numeric_limits<int32_t>::max() /
3685 1 : std::numeric_limits<uint8_t>::max())
3686 : {
3687 : using T = uint8_t;
3688 1 : constexpr int VALUES_PER_REG = 16;
3689 1 : libdivide::divider<uint32_t> fast_d(nSources);
3690 1 : const auto halfConstant = _mm_set1_epi32(nSources / 2);
3691 2 : for (int iLine = 0; iLine < nYSize; ++iLine)
3692 : {
3693 1 : T *CPL_RESTRICT pDest =
3694 : static_cast<GByte *>(pData) +
3695 1 : static_cast<GSpacing>(nLineSpace) * iLine;
3696 1 : const size_t iOffsetLine = static_cast<size_t>(iLine) * nXSize;
3697 1 : int iCol = 0;
3698 4 : for (; iCol < nXSize - (VALUES_PER_REG - 1);
3699 : iCol += VALUES_PER_REG)
3700 : {
3701 3 : __m128i reg0 = halfConstant;
3702 3 : __m128i reg1 = halfConstant;
3703 3 : __m128i reg2 = halfConstant;
3704 3 : __m128i reg3 = halfConstant;
3705 98307 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
3706 : {
3707 : const __m128i inputVal =
3708 98304 : _mm_loadu_si128(reinterpret_cast<const __m128i *>(
3709 98304 : static_cast<const T * CPL_RESTRICT>(
3710 98304 : papoSources[iSrc]) +
3711 98304 : iOffsetLine + iCol));
3712 : const __m128i low =
3713 196608 : _mm_unpacklo_epi8(inputVal, _mm_setzero_si128());
3714 : const __m128i high =
3715 196608 : _mm_unpackhi_epi8(inputVal, _mm_setzero_si128());
3716 294912 : reg0 = _mm_add_epi32(
3717 : reg0, _mm_unpacklo_epi16(low, _mm_setzero_si128()));
3718 294912 : reg1 = _mm_add_epi32(
3719 : reg1, _mm_unpackhi_epi16(low, _mm_setzero_si128()));
3720 294912 : reg2 = _mm_add_epi32(
3721 : reg2,
3722 : _mm_unpacklo_epi16(high, _mm_setzero_si128()));
3723 294912 : reg3 = _mm_add_epi32(
3724 : reg3,
3725 : _mm_unpackhi_epi16(high, _mm_setzero_si128()));
3726 : }
3727 : reg0 /= fast_d;
3728 : reg1 /= fast_d;
3729 : reg2 /= fast_d;
3730 : reg3 /= fast_d;
3731 3 : _mm_storeu_si128(
3732 3 : reinterpret_cast<__m128i *>(pDest + iCol),
3733 : _mm_packus_epi16(packus_epi32(reg0, reg1),
3734 : packus_epi32(reg2, reg3)));
3735 : }
3736 16 : for (; iCol < nXSize; ++iCol)
3737 : {
3738 15 : uint32_t acc = nSources / 2;
3739 491535 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
3740 : {
3741 491520 : acc += static_cast<const T * CPL_RESTRICT>(
3742 491520 : papoSources[iSrc])[iOffsetLine + iCol];
3743 : }
3744 15 : pDest[iCol] = static_cast<T>(acc / nSources);
3745 : }
3746 : }
3747 1 : return CE_None;
3748 : }
3749 :
3750 75 : if (!bHasNoData && eSrcType == GDT_UInt16 && eBufType == GDT_UInt16 &&
3751 153 : nPixelSpace == 2 &&
3752 5 : nSources <= std::numeric_limits<int32_t>::max() /
3753 5 : std::numeric_limits<uint16_t>::max())
3754 : {
3755 5 : libdivide::divider<uint32_t> fast_d(nSources);
3756 : using T = uint16_t;
3757 5 : const auto halfConstant = _mm_set1_epi32(nSources / 2);
3758 5 : constexpr int VALUES_PER_REG = 8;
3759 59 : for (int iLine = 0; iLine < nYSize; ++iLine)
3760 : {
3761 54 : T *CPL_RESTRICT pDest = reinterpret_cast<T *>(
3762 : static_cast<GByte *>(pData) +
3763 54 : static_cast<GSpacing>(nLineSpace) * iLine);
3764 54 : const size_t iOffsetLine = static_cast<size_t>(iLine) * nXSize;
3765 54 : int iCol = 0;
3766 366 : for (; iCol < nXSize - (VALUES_PER_REG - 1);
3767 : iCol += VALUES_PER_REG)
3768 : {
3769 312 : __m128i reg0 = halfConstant;
3770 312 : __m128i reg1 = halfConstant;
3771 99534 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
3772 : {
3773 : const __m128i inputVal =
3774 99222 : _mm_loadu_si128(reinterpret_cast<const __m128i *>(
3775 99222 : static_cast<const T * CPL_RESTRICT>(
3776 99222 : papoSources[iSrc]) +
3777 99222 : iOffsetLine + iCol));
3778 297666 : reg0 = _mm_add_epi32(
3779 : reg0,
3780 : _mm_unpacklo_epi16(inputVal, _mm_setzero_si128()));
3781 297666 : reg1 = _mm_add_epi32(
3782 : reg1,
3783 : _mm_unpackhi_epi16(inputVal, _mm_setzero_si128()));
3784 : }
3785 : reg0 /= fast_d;
3786 : reg1 /= fast_d;
3787 312 : _mm_storeu_si128(reinterpret_cast<__m128i *>(pDest + iCol),
3788 : packus_epi32(reg0, reg1));
3789 : }
3790 182 : for (; iCol < nXSize; ++iCol)
3791 : {
3792 128 : uint32_t acc = nSources / 2;
3793 229846 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
3794 : {
3795 229718 : acc += static_cast<const T * CPL_RESTRICT>(
3796 229718 : papoSources[iSrc])[iOffsetLine + iCol];
3797 : }
3798 128 : pDest[iCol] = static_cast<T>(acc / nSources);
3799 : }
3800 : }
3801 5 : return CE_None;
3802 : }
3803 :
3804 70 : if (!bHasNoData && eSrcType == GDT_Int16 && eBufType == GDT_Int16 &&
3805 143 : nPixelSpace == 2 &&
3806 7 : nSources <= std::numeric_limits<int32_t>::max() /
3807 7 : std::numeric_limits<uint16_t>::max())
3808 : {
3809 7 : libdivide::divider<uint32_t> fast_d(nSources);
3810 : using T = int16_t;
3811 7 : const auto halfConstant = _mm_set1_epi32(nSources / 2);
3812 7 : const auto shift = _mm_set1_epi16(std::numeric_limits<T>::min());
3813 7 : constexpr int VALUES_PER_REG = 8;
3814 63 : for (int iLine = 0; iLine < nYSize; ++iLine)
3815 : {
3816 56 : T *CPL_RESTRICT pDest = reinterpret_cast<T *>(
3817 : static_cast<GByte *>(pData) +
3818 56 : static_cast<GSpacing>(nLineSpace) * iLine);
3819 56 : const size_t iOffsetLine = static_cast<size_t>(iLine) * nXSize;
3820 56 : int iCol = 0;
3821 374 : for (; iCol < nXSize - (VALUES_PER_REG - 1);
3822 : iCol += VALUES_PER_REG)
3823 : {
3824 318 : __m128i reg0 = halfConstant;
3825 318 : __m128i reg1 = halfConstant;
3826 99555 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
3827 : {
3828 : // Shift input values by 32768 to get unsigned values
3829 198474 : const __m128i inputVal = _mm_add_epi16(
3830 : _mm_loadu_si128(reinterpret_cast<const __m128i *>(
3831 99237 : static_cast<const T * CPL_RESTRICT>(
3832 99237 : papoSources[iSrc]) +
3833 99237 : iOffsetLine + iCol)),
3834 : shift);
3835 297711 : reg0 = _mm_add_epi32(
3836 : reg0,
3837 : _mm_unpacklo_epi16(inputVal, _mm_setzero_si128()));
3838 297711 : reg1 = _mm_add_epi32(
3839 : reg1,
3840 : _mm_unpackhi_epi16(inputVal, _mm_setzero_si128()));
3841 : }
3842 : reg0 /= fast_d;
3843 : reg1 /= fast_d;
3844 318 : _mm_storeu_si128(
3845 318 : reinterpret_cast<__m128i *>(pDest + iCol),
3846 : _mm_add_epi16(packus_epi32(reg0, reg1), shift));
3847 : }
3848 198 : for (; iCol < nXSize; ++iCol)
3849 : {
3850 142 : int32_t acc = (-std::numeric_limits<T>::min()) * nSources +
3851 142 : nSources / 2;
3852 229895 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
3853 : {
3854 229753 : acc += static_cast<const T * CPL_RESTRICT>(
3855 229753 : papoSources[iSrc])[iOffsetLine + iCol];
3856 : }
3857 284 : pDest[iCol] = static_cast<T>(acc / nSources +
3858 142 : std::numeric_limits<T>::min());
3859 : }
3860 : }
3861 7 : return CE_None;
3862 : }
3863 : }
3864 : #endif // defined(USE_SSE2) && !defined(USE_NEON_OPTIMIZATIONS)
3865 :
3866 : #if defined(USE_SSE2)
3867 : if constexpr (std::is_same_v<Kernel, MeanKernel>)
3868 : {
3869 66 : if (!bHasNoData && eSrcType == GDT_Float32 && eBufType == GDT_Float32 &&
3870 19 : nPixelSpace == 4 && nSources > 0)
3871 : {
3872 19 : OptimizedMeanFloatSSE2<float, SSEWrapperFloat>(
3873 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
3874 19 : return CE_None;
3875 : }
3876 :
3877 47 : if (!bHasNoData && eSrcType == GDT_Float64 && eBufType == GDT_Float64 &&
3878 22 : nPixelSpace == 8 && nSources > 0)
3879 : {
3880 22 : OptimizedMeanFloatSSE2<double, SSEWrapperDouble>(
3881 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
3882 22 : return CE_None;
3883 : }
3884 : }
3885 : #endif // USE_SSE2
3886 :
3887 : /* ---- Set pixels ---- */
3888 42 : size_t ii = 0;
3889 152 : for (int iLine = 0; iLine < nYSize; ++iLine)
3890 : {
3891 3608 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
3892 : {
3893 3498 : oKernel.Reset();
3894 3498 : bool bWriteNoData = false;
3895 :
3896 12863 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
3897 : {
3898 9375 : const double dfVal = GetSrcVal(papoSources[iSrc], eSrcType, ii);
3899 :
3900 9375 : if (bHasNoData && IsNoData(dfVal, dfNoData))
3901 : {
3902 74 : if (bPropagateNoData)
3903 : {
3904 10 : bWriteNoData = true;
3905 10 : break;
3906 : }
3907 : }
3908 : else
3909 : {
3910 9301 : oKernel.ProcessPixel(dfVal);
3911 : }
3912 : }
3913 :
3914 3498 : double dfPixVal{dfNoData};
3915 3498 : if (!bWriteNoData && oKernel.HasValue())
3916 : {
3917 3470 : dfPixVal = oKernel.GetValue();
3918 : }
3919 :
3920 3498 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
3921 : static_cast<GByte *>(pData) +
3922 3498 : static_cast<GSpacing>(nLineSpace) * iLine +
3923 3498 : iCol * nPixelSpace,
3924 : eBufType, nPixelSpace, 1);
3925 : }
3926 : }
3927 :
3928 : /* ---- Return success ---- */
3929 42 : return CE_None;
3930 : } // BasicPixelFunc
3931 :
3932 : /************************************************************************/
3933 : /* GDALRegisterDefaultPixelFunc() */
3934 : /************************************************************************/
3935 :
3936 : /**
3937 : * This adds a default set of pixel functions to the global list of
3938 : * available pixel functions for derived bands:
3939 : *
3940 : * - "real": extract real part from a single raster band (just a copy if the
3941 : * input is non-complex)
3942 : * - "imag": extract imaginary part from a single raster band (0 for
3943 : * non-complex)
3944 : * - "complex": make a complex band merging two bands used as real and
3945 : * imag values
3946 : * - "polar": make a complex band using input bands for amplitude and
3947 : * phase values (b1 * exp( j * b2 ))
3948 : * - "mod": extract module from a single raster band (real or complex)
3949 : * - "phase": extract phase from a single raster band [-PI,PI] (0 or PI for
3950 : non-complex)
3951 : * - "conj": computes the complex conjugate of a single raster band (just a
3952 : * copy if the input is non-complex)
3953 : * - "sum": sum 2 or more raster bands
3954 : * - "diff": computes the difference between 2 raster bands (b1 - b2)
3955 : * - "mul": multiply 2 or more raster bands
3956 : * - "div": divide one raster band by another (b1 / b2).
3957 : * - "min": minimum value of 2 or more raster bands
3958 : * - "max": maximum value of 2 or more raster bands
3959 : * - "norm_diff": computes the normalized difference between two raster bands:
3960 : * ``(b1 - b2)/(b1 + b2)``.
3961 : * - "cmul": multiply the first band for the complex conjugate of the second
3962 : * - "inv": inverse (1./x).
3963 : * - "intensity": computes the intensity Re(x*conj(x)) of a single raster band
3964 : * (real or complex)
3965 : * - "sqrt": perform the square root of a single raster band (real only)
3966 : * - "log10": compute the logarithm (base 10) of the abs of a single raster
3967 : * band (real or complex): log10( abs( x ) )
3968 : * - "dB": perform conversion to dB of the abs of a single raster
3969 : * band (real or complex): 20. * log10( abs( x ) ).
3970 : * Note: the optional fact parameter can be set to 10. to get the
3971 : * alternative formula: 10. * log10( abs( x ) )
3972 : * - "exp": computes the exponential of each element in the input band ``x``
3973 : * (of real values): ``e ^ x``.
3974 : * The function also accepts two optional parameters: ``base`` and
3975 : ``fact``
3976 : * that allow to compute the generalized formula: ``base ^ ( fact *
3977 : x)``.
3978 : * Note: this function is the recommended one to perform conversion
3979 : * form logarithmic scale (dB): `` 10. ^ (x / 20.)``, in this case
3980 : * ``base = 10.`` and ``fact = 1./20``
3981 : * - "dB2amp": perform scale conversion from logarithmic to linear
3982 : * (amplitude) (i.e. 10 ^ ( x / 20 ) ) of a single raster
3983 : * band (real only).
3984 : * Deprecated in GDAL v3.5. Please use the ``exp`` pixel function
3985 : with
3986 : * ``base = 10.`` and ``fact = 0.05`` i.e. ``1./20``
3987 : * - "dB2pow": perform scale conversion from logarithmic to linear
3988 : * (power) (i.e. 10 ^ ( x / 10 ) ) of a single raster
3989 : * band (real only)
3990 : * Deprecated in GDAL v3.5. Please use the ``exp`` pixel function
3991 : with
3992 : * ``base = 10.`` and ``fact = 0.1`` i.e. ``1./10``
3993 : * - "pow": raise a single raster band to a constant power
3994 : * - "interpolate_linear": interpolate values between two raster bands
3995 : * using linear interpolation
3996 : * - "interpolate_exp": interpolate values between two raster bands using
3997 : * exponential interpolation
3998 : * - "scale": Apply the RasterBand metadata values of "offset" and "scale"
3999 : * - "reclassify": Reclassify values matching ranges in a table
4000 : * - "nan": Convert incoming NoData values to IEEE 754 nan
4001 : *
4002 : * @see GDALAddDerivedBandPixelFunc
4003 : *
4004 : * @return CE_None
4005 : */
4006 1478 : CPLErr GDALRegisterDefaultPixelFunc()
4007 : {
4008 1478 : GDALAddDerivedBandPixelFunc("real", RealPixelFunc);
4009 1478 : GDALAddDerivedBandPixelFunc("imag", ImagPixelFunc);
4010 1478 : GDALAddDerivedBandPixelFunc("complex", ComplexPixelFunc);
4011 1478 : GDALAddDerivedBandPixelFuncWithArgs("polar", PolarPixelFunc,
4012 : pszPolarPixelFuncMetadata);
4013 1478 : GDALAddDerivedBandPixelFunc("mod", ModulePixelFunc);
4014 1478 : GDALAddDerivedBandPixelFunc("phase", PhasePixelFunc);
4015 1478 : GDALAddDerivedBandPixelFunc("conj", ConjPixelFunc);
4016 1478 : GDALAddDerivedBandPixelFuncWithArgs("sum", SumPixelFunc,
4017 : pszSumPixelFuncMetadata);
4018 1478 : GDALAddDerivedBandPixelFuncWithArgs("diff", DiffPixelFunc,
4019 : pszDiffPixelFuncMetadata);
4020 1478 : GDALAddDerivedBandPixelFuncWithArgs("mul", MulPixelFunc,
4021 : pszMulPixelFuncMetadata);
4022 1478 : GDALAddDerivedBandPixelFuncWithArgs("div", DivPixelFunc,
4023 : pszDivPixelFuncMetadata);
4024 1478 : GDALAddDerivedBandPixelFunc("cmul", CMulPixelFunc);
4025 1478 : GDALAddDerivedBandPixelFuncWithArgs("inv", InvPixelFunc,
4026 : pszInvPixelFuncMetadata);
4027 1478 : GDALAddDerivedBandPixelFunc("intensity", IntensityPixelFunc);
4028 1478 : GDALAddDerivedBandPixelFuncWithArgs("sqrt", SqrtPixelFunc,
4029 : pszSqrtPixelFuncMetadata);
4030 1478 : GDALAddDerivedBandPixelFuncWithArgs("log10", Log10PixelFunc,
4031 : pszLog10PixelFuncMetadata);
4032 1478 : GDALAddDerivedBandPixelFuncWithArgs("dB", DBPixelFunc,
4033 : pszDBPixelFuncMetadata);
4034 1478 : GDALAddDerivedBandPixelFuncWithArgs("exp", ExpPixelFunc,
4035 : pszExpPixelFuncMetadata);
4036 1478 : GDALAddDerivedBandPixelFunc("dB2amp",
4037 : dB2AmpPixelFunc); // deprecated in v3.5
4038 1478 : GDALAddDerivedBandPixelFunc("dB2pow",
4039 : dB2PowPixelFunc); // deprecated in v3.5
4040 1478 : GDALAddDerivedBandPixelFuncWithArgs("pow", PowPixelFunc,
4041 : pszPowPixelFuncMetadata);
4042 1478 : GDALAddDerivedBandPixelFuncWithArgs("interpolate_linear",
4043 : InterpolatePixelFunc<InterpolateLinear>,
4044 : pszInterpolatePixelFuncMetadata);
4045 1478 : GDALAddDerivedBandPixelFuncWithArgs(
4046 : "interpolate_exp", InterpolatePixelFunc<InterpolateExponential>,
4047 : pszInterpolatePixelFuncMetadata);
4048 1478 : GDALAddDerivedBandPixelFuncWithArgs("replace_nodata",
4049 : ReplaceNoDataPixelFunc,
4050 : pszReplaceNoDataPixelFuncMetadata);
4051 1478 : GDALAddDerivedBandPixelFuncWithArgs("scale", ScalePixelFunc,
4052 : pszScalePixelFuncMetadata);
4053 1478 : GDALAddDerivedBandPixelFuncWithArgs("norm_diff", NormDiffPixelFunc,
4054 : pszNormDiffPixelFuncMetadata);
4055 1478 : GDALAddDerivedBandPixelFuncWithArgs("min", MinPixelFunc<ReturnValue>,
4056 : pszMinMaxFuncMetadataNodata);
4057 1478 : GDALAddDerivedBandPixelFuncWithArgs("argmin", MinPixelFunc<ReturnIndex>,
4058 : pszArgMinMaxFuncMetadataNodata);
4059 1478 : GDALAddDerivedBandPixelFuncWithArgs("max", MaxPixelFunc<ReturnValue>,
4060 : pszMinMaxFuncMetadataNodata);
4061 1478 : GDALAddDerivedBandPixelFuncWithArgs("argmax", MaxPixelFunc<ReturnIndex>,
4062 : pszArgMinMaxFuncMetadataNodata);
4063 1478 : GDALAddDerivedBandPixelFuncWithArgs("expression", ExprPixelFunc,
4064 : pszExprPixelFuncMetadata);
4065 1478 : GDALAddDerivedBandPixelFuncWithArgs("reclassify", ReclassifyPixelFunc,
4066 : pszReclassifyPixelFuncMetadata);
4067 1478 : GDALAddDerivedBandPixelFuncWithArgs("mean", BasicPixelFunc<MeanKernel>,
4068 : pszBasicPixelFuncMetadata);
4069 1478 : GDALAddDerivedBandPixelFuncWithArgs("geometric_mean",
4070 : BasicPixelFunc<GeoMeanKernel>,
4071 : pszBasicPixelFuncMetadata);
4072 1478 : GDALAddDerivedBandPixelFuncWithArgs("harmonic_mean",
4073 : BasicPixelFunc<HarmonicMeanKernel>,
4074 : pszBasicPixelFuncMetadata);
4075 1478 : GDALAddDerivedBandPixelFuncWithArgs("median", BasicPixelFunc<MedianKernel>,
4076 : pszBasicPixelFuncMetadata);
4077 1478 : GDALAddDerivedBandPixelFuncWithArgs("mode", BasicPixelFunc<ModeKernel>,
4078 : pszBasicPixelFuncMetadata);
4079 1478 : return CE_None;
4080 : }
|
