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