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