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 26199500 : inline double GetSrcVal(const void *pSource, GDALDataType eSrcType, T ii)
57 : {
58 26199500 : switch (eSrcType)
59 : {
60 0 : case GDT_Unknown:
61 0 : return 0;
62 62381 : case GDT_UInt8:
63 62381 : 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 9816 : case GDT_Float32:
84 9816 : return static_cast<const float *>(pSource)[ii];
85 26022100 : case GDT_Float64:
86 26022100 : 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 10318 : static bool IsNoData(double dfVal, double dfNoData)
104 : {
105 10318 : return dfVal == dfNoData || (std::isnan(dfVal) && std::isnan(dfNoData));
106 : }
107 :
108 2201 : static CPLErr FetchDoubleArg(CSLConstList papszArgs, const char *pszName,
109 : double *pdfX, double *pdfDefault = nullptr)
110 : {
111 2201 : const char *pszVal = CSLFetchNameValue(papszArgs, pszName);
112 :
113 2201 : if (pszVal == nullptr)
114 : {
115 1092 : 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 1092 : *pdfX = *pdfDefault;
124 1092 : return CE_None;
125 : }
126 : }
127 :
128 1109 : char *pszEnd = nullptr;
129 1109 : *pdfX = std::strtod(pszVal, &pszEnd);
130 1109 : 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 1109 : 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 8 : 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 8 : if (nSources != 1)
397 1 : return CE_Failure;
398 :
399 7 : 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 5 : double dfNoData{0};
428 5 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
429 8 : if (bHasNoData &&
430 3 : FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
431 0 : return CE_Failure;
432 :
433 : /* ---- Set pixels ---- */
434 5 : size_t ii = 0;
435 29 : for (int iLine = 0; iLine < nYSize; ++iLine)
436 : {
437 431 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
438 : {
439 407 : 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 413 : ? dfNoData
444 407 : : std::fabs(dfVal);
445 :
446 407 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
447 : static_cast<GByte *>(pData) +
448 407 : static_cast<GSpacing>(nLineSpace) * iLine +
449 407 : iCol * nPixelSpace,
450 : eBufType, nPixelSpace, 1);
451 : }
452 : }
453 : }
454 :
455 : /* ---- Return success ---- */
456 7 : 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 1927 : 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 1927 : 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 1926 : double dfK = 0.0;
1163 1926 : if (FetchDoubleArg(papszArgs, "k", &dfK, &dfK) != CE_None)
1164 0 : return CE_Failure;
1165 :
1166 1926 : double dfNoData{0};
1167 1926 : bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1168 1926 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1169 0 : return CE_Failure;
1170 :
1171 1926 : const bool bPropagateNoData = CPLTestBool(
1172 : CSLFetchNameValueDef(papszArgs, "propagateNoData", "false"));
1173 :
1174 1926 : if (dfNoData == 0 && !bPropagateNoData)
1175 1908 : bHasNoData = false;
1176 :
1177 : /* ---- Set pixels ---- */
1178 1926 : 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 1890 : bool bGeneralCase = true;
1214 1890 : 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 1890 : if (bGeneralCase)
1308 : {
1309 247 : size_t ii = 0;
1310 725 : for (int iLine = 0; iLine < nYSize; ++iLine)
1311 : {
1312 31206 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1313 : {
1314 30728 : double dfSum = dfK;
1315 92169 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
1316 : {
1317 : const double dfVal =
1318 61451 : GetSrcVal(papoSources[iSrc], eSrcType, ii);
1319 :
1320 61451 : if (bHasNoData && IsNoData(dfVal, dfNoData))
1321 : {
1322 23 : if (bPropagateNoData)
1323 : {
1324 10 : dfSum = dfNoData;
1325 10 : break;
1326 : }
1327 : }
1328 : else
1329 : {
1330 61428 : dfSum += dfVal;
1331 : }
1332 : }
1333 :
1334 30728 : GDALCopyWords(&dfSum, GDT_Float64, 0,
1335 : static_cast<GByte *>(pData) +
1336 30728 : static_cast<GSpacing>(nLineSpace) *
1337 30728 : iLine +
1338 30728 : iCol * nPixelSpace,
1339 : eBufType, nPixelSpace, 1);
1340 : }
1341 : }
1342 : }
1343 : }
1344 :
1345 : /* ---- Return success ---- */
1346 1926 : 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 6 : 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 6 : if (nSources != 2)
1361 1 : return CE_Failure;
1362 :
1363 5 : double dfNoData{0};
1364 5 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1365 5 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1366 0 : return CE_Failure;
1367 :
1368 5 : 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 4 : size_t ii = 0;
1402 13 : for (int iLine = 0; iLine < nYSize; ++iLine)
1403 : {
1404 45 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1405 : {
1406 36 : const double dfA = GetSrcVal(papoSources[0], eSrcType, ii);
1407 36 : const double dfB = GetSrcVal(papoSources[1], eSrcType, ii);
1408 :
1409 : const double dfPixVal =
1410 40 : bHasNoData &&
1411 4 : (IsNoData(dfA, dfNoData) || IsNoData(dfB, dfNoData))
1412 40 : ? dfNoData
1413 36 : : dfA - dfB;
1414 :
1415 36 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1416 : static_cast<GByte *>(pData) +
1417 36 : static_cast<GSpacing>(nLineSpace) * iLine +
1418 36 : iCol * nPixelSpace,
1419 : eBufType, nPixelSpace, 1);
1420 : }
1421 : }
1422 : }
1423 :
1424 : /* ---- Return success ---- */
1425 5 : 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 15 : 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 15 : 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 14 : double dfK = 1.0;
1452 14 : if (FetchDoubleArg(papszArgs, "k", &dfK, &dfK) != CE_None)
1453 0 : return CE_Failure;
1454 :
1455 14 : double dfNoData{0};
1456 14 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1457 14 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1458 0 : return CE_Failure;
1459 :
1460 14 : const bool bPropagateNoData = CPLTestBool(
1461 : CSLFetchNameValueDef(papszArgs, "propagateNoData", "false"));
1462 :
1463 : /* ---- Set pixels ---- */
1464 14 : 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 13 : size_t ii = 0;
1506 83 : for (int iLine = 0; iLine < nYSize; ++iLine)
1507 : {
1508 1290 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1509 : {
1510 1220 : double dfPixVal = dfK; // Not complex.
1511 :
1512 4051 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
1513 : {
1514 : const double dfVal =
1515 2837 : GetSrcVal(papoSources[iSrc], eSrcType, ii);
1516 :
1517 2837 : if (bHasNoData && IsNoData(dfVal, dfNoData))
1518 : {
1519 18 : if (bPropagateNoData)
1520 : {
1521 6 : dfPixVal = dfNoData;
1522 6 : break;
1523 : }
1524 : }
1525 : else
1526 : {
1527 2819 : dfPixVal *= dfVal;
1528 : }
1529 : }
1530 :
1531 1220 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1532 : static_cast<GByte *>(pData) +
1533 1220 : static_cast<GSpacing>(nLineSpace) * iLine +
1534 1220 : iCol * nPixelSpace,
1535 : eBufType, nPixelSpace, 1);
1536 : }
1537 : }
1538 : }
1539 :
1540 : /* ---- Return success ---- */
1541 14 : 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 8 : 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 8 : if (nSources != 2)
1557 0 : return CE_Failure;
1558 :
1559 8 : double dfNoData{0};
1560 8 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1561 8 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1562 0 : return CE_Failure;
1563 :
1564 : /* ---- Set pixels ---- */
1565 8 : 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 7 : size_t ii = 0;
1608 19 : for (int iLine = 0; iLine < nYSize; ++iLine)
1609 : {
1610 51 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1611 : {
1612 39 : const double dfNum = GetSrcVal(papoSources[0], eSrcType, ii);
1613 39 : const double dfDenom = GetSrcVal(papoSources[1], eSrcType, ii);
1614 :
1615 39 : double dfPixVal = dfNoData;
1616 43 : if (!bHasNoData || (!IsNoData(dfNum, dfNoData) &&
1617 4 : !IsNoData(dfDenom, dfNoData)))
1618 : {
1619 : // coverity[divide_by_zero]
1620 35 : dfPixVal =
1621 : dfDenom == 0
1622 35 : ? 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 39 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1633 : static_cast<GByte *>(pData) +
1634 39 : static_cast<GSpacing>(nLineSpace) * iLine +
1635 39 : iCol * nPixelSpace,
1636 : eBufType, nPixelSpace, 1);
1637 : }
1638 : }
1639 : }
1640 :
1641 : /* ---- Return success ---- */
1642 8 : 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 13 : 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 13 : if (nSources != 1)
1729 1 : return CE_Failure;
1730 :
1731 12 : double dfK = 1.0;
1732 12 : if (FetchDoubleArg(papszArgs, "k", &dfK, &dfK) != CE_None)
1733 0 : return CE_Failure;
1734 :
1735 12 : double dfNoData{0};
1736 12 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1737 12 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1738 0 : return CE_Failure;
1739 :
1740 : /* ---- Set pixels ---- */
1741 12 : 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 10 : size_t ii = 0;
1775 58 : for (int iLine = 0; iLine < nYSize; ++iLine)
1776 : {
1777 860 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1778 : {
1779 : // Source raster pixels may be obtained with GetSrcVal macro.
1780 : // Not complex.
1781 812 : const double dfVal = GetSrcVal(papoSources[0], eSrcType, ii);
1782 812 : double dfPixVal = dfNoData;
1783 :
1784 812 : if (!bHasNoData || !IsNoData(dfVal, dfNoData))
1785 : {
1786 807 : dfPixVal =
1787 : dfVal == 0
1788 807 : ? std::numeric_limits<double>::infinity()
1789 806 : : dfK /
1790 : #ifdef __COVERITY__
1791 : (dfVal + std::numeric_limits<double>::min())
1792 : #else
1793 : dfVal
1794 : #endif
1795 : ;
1796 : }
1797 :
1798 812 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1799 : static_cast<GByte *>(pData) +
1800 812 : static_cast<GSpacing>(nLineSpace) * iLine +
1801 812 : iCol * nPixelSpace,
1802 : eBufType, nPixelSpace, 1);
1803 : }
1804 : }
1805 : }
1806 :
1807 : /* ---- Return success ---- */
1808 12 : 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 4 : 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 4 : if (nSources != 1)
1884 1 : return CE_Failure;
1885 3 : if (GDALDataTypeIsComplex(eSrcType))
1886 0 : return CE_Failure;
1887 :
1888 3 : double dfNoData{0};
1889 3 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1890 3 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1891 0 : return CE_Failure;
1892 :
1893 : /* ---- Set pixels ---- */
1894 3 : size_t ii = 0;
1895 25 : for (int iLine = 0; iLine < nYSize; ++iLine)
1896 : {
1897 425 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1898 : {
1899 : // Source raster pixels may be obtained with GetSrcVal macro.
1900 403 : double dfPixVal = GetSrcVal(papoSources[0], eSrcType, ii);
1901 :
1902 403 : if (bHasNoData && IsNoData(dfPixVal, dfNoData))
1903 : {
1904 2 : dfPixVal = dfNoData;
1905 : }
1906 : else
1907 : {
1908 401 : dfPixVal = std::sqrt(dfPixVal);
1909 : }
1910 :
1911 403 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1912 : static_cast<GByte *>(pData) +
1913 403 : static_cast<GSpacing>(nLineSpace) * iLine +
1914 403 : iCol * nPixelSpace,
1915 : eBufType, nPixelSpace, 1);
1916 : }
1917 : }
1918 :
1919 : /* ---- Return success ---- */
1920 3 : return CE_None;
1921 : } // SqrtPixelFunc
1922 :
1923 14 : 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 14 : if (nSources != 1)
1931 2 : return CE_Failure;
1932 :
1933 12 : double dfNoData{0};
1934 12 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1935 12 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1936 0 : return CE_Failure;
1937 :
1938 12 : 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 7 : size_t ii = 0;
1977 90 : for (int iLine = 0; iLine < nYSize; ++iLine)
1978 : {
1979 1688 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1980 : {
1981 : // Source raster pixels may be obtained with GetSrcVal macro.
1982 1605 : const double dfSrcVal = GetSrcVal(papoSources[0], eSrcType, ii);
1983 : const double dfPixVal =
1984 4 : bHasNoData && IsNoData(dfSrcVal, dfNoData)
1985 1609 : ? dfNoData
1986 1601 : : fact * std::log10(std::abs(dfSrcVal));
1987 :
1988 1605 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1989 : static_cast<GByte *>(pData) +
1990 1605 : static_cast<GSpacing>(nLineSpace) * iLine +
1991 1605 : iCol * nPixelSpace,
1992 : eBufType, nPixelSpace, 1);
1993 : }
1994 : }
1995 : }
1996 :
1997 : /* ---- Return success ---- */
1998 12 : 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 6 : 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 6 : return Log10PixelFuncHelper(papoSources, nSources, pData, nXSize, nYSize,
2012 : eSrcType, eBufType, nPixelSpace, nLineSpace,
2013 6 : 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 9 : 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 9 : if (nSources != 1)
2045 2 : return CE_Failure;
2046 7 : if (GDALDataTypeIsComplex(eSrcType))
2047 0 : return CE_Failure;
2048 :
2049 7 : double dfNoData{0};
2050 7 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
2051 7 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
2052 0 : return CE_Failure;
2053 :
2054 : /* ---- Set pixels ---- */
2055 7 : size_t ii = 0;
2056 109 : for (int iLine = 0; iLine < nYSize; ++iLine)
2057 : {
2058 2105 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
2059 : {
2060 : // Source raster pixels may be obtained with GetSrcVal macro.
2061 2003 : const double dfVal = GetSrcVal(papoSources[0], eSrcType, ii);
2062 2 : const double dfPixVal = bHasNoData && IsNoData(dfVal, dfNoData)
2063 2005 : ? dfNoData
2064 2001 : : pow(base, dfVal * fact);
2065 :
2066 2003 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
2067 : static_cast<GByte *>(pData) +
2068 2003 : static_cast<GSpacing>(nLineSpace) * iLine +
2069 2003 : iCol * nPixelSpace,
2070 : eBufType, nPixelSpace, 1);
2071 : }
2072 : }
2073 :
2074 : /* ---- Return success ---- */
2075 7 : 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 5 : 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 5 : double dfBase = 2.7182818284590452353602874713526624;
2092 5 : double dfFact = 1.;
2093 :
2094 5 : if (FetchDoubleArg(papszArgs, "base", &dfBase, &dfBase) != CE_None)
2095 0 : return CE_Failure;
2096 :
2097 5 : if (FetchDoubleArg(papszArgs, "fact", &dfFact, &dfFact) != CE_None)
2098 0 : return CE_Failure;
2099 :
2100 5 : return ExpPixelFuncHelper(papoSources, nSources, pData, nXSize, nYSize,
2101 : eSrcType, eBufType, nPixelSpace, nLineSpace,
2102 5 : 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 3 : 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 3 : if (nSources != 1)
2139 0 : return CE_Failure;
2140 3 : if (GDALDataTypeIsComplex(eSrcType))
2141 0 : return CE_Failure;
2142 :
2143 : double power;
2144 3 : if (FetchDoubleArg(papszArgs, "power", &power) != CE_None)
2145 0 : return CE_Failure;
2146 :
2147 3 : double dfNoData{0};
2148 3 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
2149 3 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
2150 0 : return CE_Failure;
2151 :
2152 : /* ---- Set pixels ---- */
2153 3 : size_t ii = 0;
2154 25 : for (int iLine = 0; iLine < nYSize; ++iLine)
2155 : {
2156 425 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
2157 : {
2158 403 : const double dfVal = GetSrcVal(papoSources[0], eSrcType, ii);
2159 :
2160 2 : const double dfPixVal = bHasNoData && IsNoData(dfVal, dfNoData)
2161 405 : ? dfNoData
2162 401 : : std::pow(dfVal, power);
2163 :
2164 403 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
2165 : static_cast<GByte *>(pData) +
2166 403 : static_cast<GSpacing>(nLineSpace) * iLine +
2167 403 : iCol * nPixelSpace,
2168 : eBufType, nPixelSpace, 1);
2169 : }
2170 : }
2171 :
2172 : /* ---- Return success ---- */
2173 3 : 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 32 : 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 32 : 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 32 : double dfNoData = std::numeric_limits<double>::quiet_NaN();
2517 32 : if (FetchDoubleArg(papszArgs, "NoData", &dfNoData, &dfNoData) != CE_None)
2518 0 : return CE_Failure;
2519 32 : const bool bPropagateNoData = CPLTestBool(
2520 : CSLFetchNameValueDef(papszArgs, "propagateNoData", "false"));
2521 :
2522 : /* ---- Set pixels ---- */
2523 32 : size_t ii = 0;
2524 260 : for (int iLine = 0; iLine < nYSize; ++iLine)
2525 : {
2526 10278 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
2527 : {
2528 10050 : double dfRes = std::numeric_limits<double>::quiet_NaN();
2529 10050 : double dfResSrc = std::numeric_limits<double>::quiet_NaN();
2530 :
2531 33568 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
2532 : {
2533 25580 : const double dfVal = GetSrcVal(papoSources[iSrc], eSrcType, ii);
2534 :
2535 25580 : 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 11155 : else if (Comparator::compare(dfVal, dfRes))
2548 : {
2549 6588 : 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 10038 : if (std::isnan(dfRes))
2565 : {
2566 3211 : dfRes = dfNoData;
2567 : }
2568 :
2569 10038 : if (IsNoData(dfRes, dfNoData))
2570 : {
2571 5269 : if (!bPropagateNoData && !std::isnan(dfK))
2572 : {
2573 6 : dfRes = dfK;
2574 : }
2575 : }
2576 4769 : else if (!std::isnan(dfK) && Comparator::compare(dfK, dfRes))
2577 : {
2578 6 : dfRes = dfK;
2579 : }
2580 : }
2581 :
2582 10050 : GDALCopyWords(&dfRes, GDT_Float64, 0,
2583 : static_cast<GByte *>(pData) +
2584 10050 : static_cast<GSpacing>(nLineSpace) * iLine +
2585 10050 : iCol * nPixelSpace,
2586 : eBufType, nPixelSpace, 1);
2587 : }
2588 : }
2589 :
2590 : /* ---- Return success ---- */
2591 32 : return CE_None;
2592 : } /* MinOrMaxPixelFunc */
2593 :
2594 : #ifdef USE_SSE2
2595 :
2596 : template <class T, class SSEWrapper>
2597 23 : static void OptimizedMinOrMaxSSE2(const void *const *papoSources, int nSources,
2598 : void *pData, int nXSize, int nYSize,
2599 : int nLineSpace)
2600 : {
2601 23 : assert(nSources >= 1);
2602 23 : constexpr int VALUES_PER_REG =
2603 : static_cast<int>(sizeof(typename SSEWrapper::Vec) / sizeof(T));
2604 585 : for (int iLine = 0; iLine < nYSize; ++iLine)
2605 : {
2606 562 : T *CPL_RESTRICT pDest =
2607 : reinterpret_cast<T *>(static_cast<GByte *>(pData) +
2608 562 : static_cast<GSpacing>(nLineSpace) * iLine);
2609 562 : const size_t iOffsetLine = static_cast<size_t>(iLine) * nXSize;
2610 562 : int iCol = 0;
2611 3124 : for (; iCol < nXSize - (2 * VALUES_PER_REG - 1);
2612 : iCol += 2 * VALUES_PER_REG)
2613 : {
2614 2562 : auto reg0 = SSEWrapper::LoadU(
2615 2562 : static_cast<const T * CPL_RESTRICT>(papoSources[0]) +
2616 2562 : iOffsetLine + iCol);
2617 2562 : auto reg1 = SSEWrapper::LoadU(
2618 2562 : static_cast<const T * CPL_RESTRICT>(papoSources[0]) +
2619 2562 : iOffsetLine + iCol + VALUES_PER_REG);
2620 5186 : for (int iSrc = 1; iSrc < nSources; ++iSrc)
2621 : {
2622 2624 : reg0 = SSEWrapper::MinOrMax(
2623 2624 : reg0, SSEWrapper::LoadU(static_cast<const T * CPL_RESTRICT>(
2624 2624 : papoSources[iSrc]) +
2625 2624 : iOffsetLine + iCol));
2626 2624 : reg1 = SSEWrapper::MinOrMax(
2627 : reg1,
2628 : SSEWrapper::LoadU(
2629 2624 : static_cast<const T * CPL_RESTRICT>(papoSources[iSrc]) +
2630 2624 : iOffsetLine + iCol + VALUES_PER_REG));
2631 : }
2632 2562 : SSEWrapper::StoreU(pDest + iCol, reg0);
2633 2562 : SSEWrapper::StoreU(pDest + iCol + VALUES_PER_REG, reg1);
2634 : }
2635 4072 : for (; iCol < nXSize; ++iCol)
2636 : {
2637 3510 : T v = static_cast<const T * CPL_RESTRICT>(
2638 3510 : papoSources[0])[iOffsetLine + iCol];
2639 7926 : for (int iSrc = 1; iSrc < nSources; ++iSrc)
2640 : {
2641 4416 : v = SSEWrapper::MinOrMax(
2642 4416 : v, static_cast<const T * CPL_RESTRICT>(
2643 4416 : papoSources[iSrc])[iOffsetLine + iCol]);
2644 : }
2645 3510 : pDest[iCol] = v;
2646 : }
2647 : }
2648 23 : }
2649 :
2650 : /************************************************************************/
2651 : /* NaNAwareMinOrMaxFloat/Double() */
2652 : /************************************************************************/
2653 :
2654 1216 : static inline __m128 NaNAwareMinOrMaxFloat(__m128 x, __m128 y, bool bMin)
2655 : {
2656 1216 : const __m128 xIsNaN = _mm_cmpunord_ps(x, x);
2657 1216 : const __m128 yIsNaN = _mm_cmpunord_ps(y, y);
2658 : const __m128 xx =
2659 3648 : _mm_or_ps(_mm_andnot_ps(xIsNaN, x), _mm_and_ps(xIsNaN, y));
2660 : const __m128 yy =
2661 2432 : _mm_or_ps(_mm_andnot_ps(yIsNaN, y), _mm_and_ps(yIsNaN, x));
2662 2432 : return bMin ? _mm_min_ps(xx, yy) : _mm_max_ps(xx, yy);
2663 : }
2664 :
2665 2432 : static inline __m128d NaNAwareMinOrMaxDouble(__m128d x, __m128d y, bool bMin)
2666 : {
2667 2432 : const __m128d xIsNaN = _mm_cmpunord_pd(x, x);
2668 2432 : const __m128d yIsNaN = _mm_cmpunord_pd(y, y);
2669 : const __m128d xx =
2670 7296 : _mm_or_pd(_mm_andnot_pd(xIsNaN, x), _mm_and_pd(xIsNaN, y));
2671 : const __m128d yy =
2672 4864 : _mm_or_pd(_mm_andnot_pd(yIsNaN, y), _mm_and_pd(yIsNaN, x));
2673 4864 : return bMin ? _mm_min_pd(xx, yy) : _mm_max_pd(xx, yy);
2674 : }
2675 :
2676 408 : template <class T> static inline T NaNAwareMinOrMax(T x, T y, bool bMin)
2677 : {
2678 408 : if (std::isnan(x))
2679 0 : return y;
2680 408 : if (std::isnan(y))
2681 0 : return x;
2682 408 : return bMin ? std::min(x, y) : std::max(x, y);
2683 : }
2684 :
2685 : // clang-format off
2686 : namespace
2687 : {
2688 : struct SSEWrapperMinByte
2689 : {
2690 : using T = uint8_t;
2691 : typedef __m128i Vec;
2692 :
2693 400 : static inline Vec LoadU(const T *x) { return _mm_loadu_si128(reinterpret_cast<const Vec*>(x)); }
2694 100 : static inline void StoreU(T *x, Vec y) { _mm_storeu_si128(reinterpret_cast<Vec*>(x), y); }
2695 200 : static inline Vec MinOrMax(Vec x, Vec y) { return _mm_min_epu8(x, y); }
2696 904 : static inline T MinOrMax(T x, T y) { return std::min(x, y); }
2697 : };
2698 :
2699 : struct SSEWrapperMaxByte
2700 : {
2701 : using T = uint8_t;
2702 : typedef __m128i Vec;
2703 :
2704 1000 : static inline Vec LoadU(const T *x) { return _mm_loadu_si128(reinterpret_cast<const Vec*>(x)); }
2705 200 : static inline void StoreU(T *x, Vec y) { _mm_storeu_si128(reinterpret_cast<Vec*>(x), y); }
2706 600 : static inline Vec MinOrMax(Vec x, Vec y) { return _mm_max_epu8(x, y); }
2707 2704 : static inline T MinOrMax(T x, T y) { return std::max(x, y); }
2708 : };
2709 :
2710 : struct SSEWrapperMinUInt16
2711 : {
2712 : using T = uint16_t;
2713 : typedef __m128i Vec;
2714 :
2715 1200 : static inline Vec LoadU(const T *x) { return _mm_loadu_si128(reinterpret_cast<const Vec*>(x)); }
2716 300 : static inline void StoreU(T *x, Vec y) { _mm_storeu_si128(reinterpret_cast<Vec*>(x), y); }
2717 : #if defined(__SSE4_1__) || defined(USE_NEON_OPTIMIZATIONS)
2718 : static inline Vec MinOrMax(Vec x, Vec y) { return _mm_min_epu16(x, y); }
2719 : #else
2720 300 : static inline Vec MinOrMax(Vec x, Vec y) { return
2721 1800 : _mm_add_epi16(
2722 : _mm_min_epi16(
2723 : _mm_add_epi16(x, _mm_set1_epi16(-32768)),
2724 : _mm_add_epi16(y, _mm_set1_epi16(-32768))),
2725 300 : _mm_set1_epi16(-32768)); }
2726 : #endif
2727 100 : static inline T MinOrMax(T x, T y) { return std::min(x, y); }
2728 : };
2729 :
2730 : struct SSEWrapperMaxUInt16
2731 : {
2732 : using T = uint16_t;
2733 : typedef __m128i Vec;
2734 :
2735 1200 : static inline Vec LoadU(const T *x) { return _mm_loadu_si128(reinterpret_cast<const Vec*>(x)); }
2736 300 : static inline void StoreU(T *x, Vec y) { _mm_storeu_si128(reinterpret_cast<Vec*>(x), y); }
2737 : #if defined(__SSE4_1__) || defined(USE_NEON_OPTIMIZATIONS)
2738 : static inline Vec MinOrMax(Vec x, Vec y) { return _mm_max_epu16(x, y); }
2739 : #else
2740 300 : static inline Vec MinOrMax(Vec x, Vec y) { return
2741 1800 : _mm_add_epi16(
2742 : _mm_max_epi16(
2743 : _mm_add_epi16(x, _mm_set1_epi16(-32768)),
2744 : _mm_add_epi16(y, _mm_set1_epi16(-32768))),
2745 300 : _mm_set1_epi16(-32768)); }
2746 : #endif
2747 100 : static inline T MinOrMax(T x, T y) { return std::max(x, y); }
2748 : };
2749 :
2750 : struct SSEWrapperMinInt16
2751 : {
2752 : using T = int16_t;
2753 : typedef __m128i Vec;
2754 :
2755 1200 : static inline Vec LoadU(const T *x) { return _mm_loadu_si128(reinterpret_cast<const Vec*>(x)); }
2756 300 : static inline void StoreU(T *x, Vec y) { _mm_storeu_si128(reinterpret_cast<Vec*>(x), y); }
2757 600 : static inline Vec MinOrMax(Vec x, Vec y) { return _mm_min_epi16(x, y); }
2758 100 : static inline T MinOrMax(T x, T y) { return std::min(x, y); }
2759 : };
2760 :
2761 : struct SSEWrapperMaxInt16
2762 : {
2763 : using T = int16_t;
2764 : typedef __m128i Vec;
2765 :
2766 1200 : static inline Vec LoadU(const T *x) { return _mm_loadu_si128(reinterpret_cast<const Vec*>(x)); }
2767 300 : static inline void StoreU(T *x, Vec y) { _mm_storeu_si128(reinterpret_cast<Vec*>(x), y); }
2768 600 : static inline Vec MinOrMax(Vec x, Vec y) { return _mm_max_epi16(x, y); }
2769 100 : static inline T MinOrMax(T x, T y) { return std::max(x, y); }
2770 : };
2771 :
2772 : struct SSEWrapperMinFloat
2773 : {
2774 : using T = float;
2775 : typedef __m128 Vec;
2776 :
2777 2424 : static inline Vec LoadU(const T *x) { return _mm_loadu_ps(x); }
2778 604 : static inline void StoreU(T *x, Vec y) { _mm_storeu_ps(x, y); }
2779 608 : static inline Vec MinOrMax(Vec x, Vec y) { return NaNAwareMinOrMaxFloat(x, y, true); }
2780 100 : static inline T MinOrMax(T x, T y) { return NaNAwareMinOrMax(x, y, true); }
2781 : };
2782 :
2783 : struct SSEWrapperMaxFloat
2784 : {
2785 : using T = float;
2786 : typedef __m128 Vec;
2787 :
2788 2424 : static inline Vec LoadU(const T *x) { return _mm_loadu_ps(x); }
2789 604 : static inline void StoreU(T *x, Vec y) { _mm_storeu_ps(x, y); }
2790 608 : static inline Vec MinOrMax(Vec x, Vec y) { return NaNAwareMinOrMaxFloat(x, y, false); }
2791 100 : static inline T MinOrMax(T x, T y) { return NaNAwareMinOrMax(x, y, false); }
2792 : };
2793 :
2794 : struct SSEWrapperMinDouble
2795 : {
2796 : using T = double;
2797 : typedef __m128d Vec;
2798 :
2799 4848 : static inline Vec LoadU(const T *x) { return _mm_loadu_pd(x); }
2800 1208 : static inline void StoreU(T *x, Vec y) { _mm_storeu_pd(x, y); }
2801 1216 : static inline Vec MinOrMax(Vec x, Vec y) { return NaNAwareMinOrMaxDouble(x, y, true); }
2802 104 : static inline T MinOrMax(T x, T y) { return NaNAwareMinOrMax(x, y, true); }
2803 : };
2804 :
2805 : struct SSEWrapperMaxDouble
2806 : {
2807 : using T = double;
2808 : typedef __m128d Vec;
2809 :
2810 4848 : static inline Vec LoadU(const T *x) { return _mm_loadu_pd(x); }
2811 1208 : static inline void StoreU(T *x, Vec y) { _mm_storeu_pd(x, y); }
2812 1216 : static inline Vec MinOrMax(Vec x, Vec y) { return NaNAwareMinOrMaxDouble(x, y, false); }
2813 104 : static inline T MinOrMax(T x, T y) { return NaNAwareMinOrMax(x, y, false); }
2814 : };
2815 :
2816 : } // namespace
2817 :
2818 : // clang-format on
2819 :
2820 : #endif // USE_SSE2
2821 :
2822 : template <typename ReturnType>
2823 30 : static CPLErr MinPixelFunc(void **papoSources, int nSources, void *pData,
2824 : int nXSize, int nYSize, GDALDataType eSrcType,
2825 : GDALDataType eBufType, int nPixelSpace,
2826 : int nLineSpace, CSLConstList papszArgs)
2827 : {
2828 : struct Comparator
2829 : {
2830 2738 : static bool compare(double x, double resVal)
2831 : {
2832 : // Written this way to deal with resVal being NaN
2833 2738 : return !(x >= resVal);
2834 : }
2835 : };
2836 :
2837 30 : double dfK = std::numeric_limits<double>::quiet_NaN();
2838 : if constexpr (std::is_same_v<ReturnType, ReturnValue>)
2839 : {
2840 27 : if (FetchDoubleArg(papszArgs, "k", &dfK, &dfK) != CE_None)
2841 0 : return CE_Failure;
2842 :
2843 : #ifdef USE_SSE2
2844 27 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
2845 46 : if (std::isnan(dfK) && nSources > 0 && !bHasNoData &&
2846 46 : eSrcType == eBufType &&
2847 11 : nPixelSpace == GDALGetDataTypeSizeBytes(eSrcType))
2848 : {
2849 11 : if (eSrcType == GDT_UInt8)
2850 : {
2851 3 : OptimizedMinOrMaxSSE2<uint8_t, SSEWrapperMinByte>(
2852 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
2853 3 : return CE_None;
2854 : }
2855 8 : else if (eSrcType == GDT_UInt16)
2856 : {
2857 1 : OptimizedMinOrMaxSSE2<uint16_t, SSEWrapperMinUInt16>(
2858 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
2859 1 : return CE_None;
2860 : }
2861 7 : else if (eSrcType == GDT_Int16)
2862 : {
2863 1 : OptimizedMinOrMaxSSE2<int16_t, SSEWrapperMinInt16>(
2864 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
2865 1 : return CE_None;
2866 : }
2867 6 : else if (eSrcType == GDT_Float32)
2868 : {
2869 2 : OptimizedMinOrMaxSSE2<float, SSEWrapperMinFloat>(
2870 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
2871 2 : return CE_None;
2872 : }
2873 4 : else if (eSrcType == GDT_Float64)
2874 : {
2875 4 : OptimizedMinOrMaxSSE2<double, SSEWrapperMinDouble>(
2876 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
2877 4 : return CE_None;
2878 : }
2879 : }
2880 : #endif
2881 : }
2882 :
2883 19 : return MinOrMaxPixelFunc<Comparator, ReturnType>(
2884 : dfK, papoSources, nSources, pData, nXSize, nYSize, eSrcType, eBufType,
2885 22 : nPixelSpace, nLineSpace, papszArgs);
2886 : }
2887 :
2888 : template <typename ReturnType>
2889 25 : static CPLErr MaxPixelFunc(void **papoSources, int nSources, void *pData,
2890 : int nXSize, int nYSize, GDALDataType eSrcType,
2891 : GDALDataType eBufType, int nPixelSpace,
2892 : int nLineSpace, CSLConstList papszArgs)
2893 : {
2894 : struct Comparator
2895 : {
2896 8429 : static bool compare(double x, double resVal)
2897 : {
2898 : // Written this way to deal with resVal being NaN
2899 8429 : return !(x <= resVal);
2900 : }
2901 : };
2902 :
2903 25 : double dfK = std::numeric_limits<double>::quiet_NaN();
2904 : if constexpr (std::is_same_v<ReturnType, ReturnValue>)
2905 : {
2906 22 : if (FetchDoubleArg(papszArgs, "k", &dfK, &dfK) != CE_None)
2907 0 : return CE_Failure;
2908 :
2909 : #ifdef USE_SSE2
2910 22 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
2911 38 : if (std::isnan(dfK) && nSources > 0 && !bHasNoData &&
2912 38 : eSrcType == eBufType &&
2913 12 : nPixelSpace == GDALGetDataTypeSizeBytes(eSrcType))
2914 : {
2915 12 : if (eSrcType == GDT_UInt8)
2916 : {
2917 4 : OptimizedMinOrMaxSSE2<uint8_t, SSEWrapperMaxByte>(
2918 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
2919 4 : return CE_None;
2920 : }
2921 8 : else if (eSrcType == GDT_UInt16)
2922 : {
2923 1 : OptimizedMinOrMaxSSE2<uint16_t, SSEWrapperMaxUInt16>(
2924 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
2925 1 : return CE_None;
2926 : }
2927 7 : else if (eSrcType == GDT_Int16)
2928 : {
2929 1 : OptimizedMinOrMaxSSE2<int16_t, SSEWrapperMaxInt16>(
2930 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
2931 1 : return CE_None;
2932 : }
2933 6 : else if (eSrcType == GDT_Float32)
2934 : {
2935 2 : OptimizedMinOrMaxSSE2<float, SSEWrapperMaxFloat>(
2936 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
2937 2 : return CE_None;
2938 : }
2939 4 : else if (eSrcType == GDT_Float64)
2940 : {
2941 4 : OptimizedMinOrMaxSSE2<double, SSEWrapperMaxDouble>(
2942 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
2943 4 : return CE_None;
2944 : }
2945 : }
2946 : #endif
2947 : }
2948 :
2949 13 : return MinOrMaxPixelFunc<Comparator, ReturnType>(
2950 : dfK, papoSources, nSources, pData, nXSize, nYSize, eSrcType, eBufType,
2951 16 : nPixelSpace, nLineSpace, papszArgs);
2952 : }
2953 :
2954 : static const char pszExprPixelFuncMetadata[] =
2955 : "<PixelFunctionArgumentsList>"
2956 : " <Argument type='builtin' value='NoData' optional='true' />"
2957 : " <Argument name='propagateNoData' description='Whether the output value "
2958 : "should be NoData as as soon as one source is NoData' type='boolean' "
2959 : "default='false' />"
2960 : " <Argument name='expression' "
2961 : " description='Expression to be evaluated' "
2962 : " type='string'></Argument>"
2963 : " <Argument name='dialect' "
2964 : " description='Expression dialect' "
2965 : " type='string-select'"
2966 : " default='muparser'>"
2967 : " <Value>exprtk</Value>"
2968 : " <Value>muparser</Value>"
2969 : " </Argument>"
2970 : " <Argument type='builtin' value='source_names' />"
2971 : " <Argument type='builtin' value='xoff' />"
2972 : " <Argument type='builtin' value='yoff' />"
2973 : " <Argument type='builtin' value='geotransform' />"
2974 : "</PixelFunctionArgumentsList>";
2975 :
2976 367 : static CPLErr ExprPixelFunc(void **papoSources, int nSources, void *pData,
2977 : int nXSize, int nYSize, GDALDataType eSrcType,
2978 : GDALDataType eBufType, int nPixelSpace,
2979 : int nLineSpace, CSLConstList papszArgs)
2980 : {
2981 : /* ---- Init ---- */
2982 367 : if (GDALDataTypeIsComplex(eSrcType))
2983 : {
2984 0 : CPLError(CE_Failure, CPLE_AppDefined,
2985 : "expression cannot by applied to complex data types");
2986 0 : return CE_Failure;
2987 : }
2988 :
2989 367 : double dfNoData{0};
2990 367 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
2991 367 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
2992 0 : return CE_Failure;
2993 :
2994 367 : const bool bPropagateNoData = CPLTestBool(
2995 : CSLFetchNameValueDef(papszArgs, "propagateNoData", "false"));
2996 :
2997 367 : const char *pszExpression = CSLFetchNameValue(papszArgs, "expression");
2998 367 : if (!pszExpression)
2999 : {
3000 1 : CPLError(CE_Failure, CPLE_AppDefined,
3001 : "Missing 'expression' pixel function argument");
3002 1 : return CE_Failure;
3003 : }
3004 :
3005 366 : const char *pszSourceNames = CSLFetchNameValue(papszArgs, "source_names");
3006 : const CPLStringList aosSourceNames(
3007 732 : CSLTokenizeString2(pszSourceNames, "|", 0));
3008 366 : if (aosSourceNames.size() != nSources)
3009 : {
3010 7 : CPLError(CE_Failure, CPLE_AppDefined,
3011 : "The source_names variable passed to ExprPixelFunc() has %d "
3012 : "values, whereas %d were expected. An invalid variable name "
3013 : "has likely been used",
3014 : aosSourceNames.size(), nSources);
3015 7 : return CE_Failure;
3016 : }
3017 :
3018 718 : std::vector<double> adfValuesForPixel(nSources);
3019 :
3020 359 : const char *pszDialect = CSLFetchNameValue(papszArgs, "dialect");
3021 359 : if (!pszDialect)
3022 : {
3023 220 : pszDialect = "muparser";
3024 : }
3025 :
3026 718 : auto poExpression = gdal::MathExpression::Create(pszExpression, pszDialect);
3027 :
3028 : // cppcheck-suppress knownConditionTrueFalse
3029 359 : if (!poExpression)
3030 : {
3031 0 : return CE_Failure;
3032 : }
3033 :
3034 359 : int nXOff = 0;
3035 359 : int nYOff = 0;
3036 359 : GDALGeoTransform gt;
3037 359 : double dfCenterX = 0;
3038 359 : double dfCenterY = 0;
3039 :
3040 359 : bool includeCenterCoords = false;
3041 359 : if (strstr(pszExpression, "_CENTER_X_") ||
3042 357 : strstr(pszExpression, "_CENTER_Y_"))
3043 : {
3044 2 : includeCenterCoords = true;
3045 :
3046 2 : const char *pszXOff = CSLFetchNameValue(papszArgs, "xoff");
3047 2 : nXOff = std::atoi(pszXOff);
3048 :
3049 2 : const char *pszYOff = CSLFetchNameValue(papszArgs, "yoff");
3050 2 : nYOff = std::atoi(pszYOff);
3051 :
3052 2 : const char *pszGT = CSLFetchNameValue(papszArgs, "geotransform");
3053 2 : if (pszGT == nullptr)
3054 : {
3055 1 : CPLError(CE_Failure, CPLE_AppDefined,
3056 : "To use _CENTER_X_ or _CENTER_Y_ in an expression, "
3057 : "VRTDataset must have a <GeoTransform> element.");
3058 1 : return CE_Failure;
3059 : }
3060 :
3061 1 : if (!gt.Init(pszGT))
3062 : {
3063 0 : CPLError(CE_Failure, CPLE_AppDefined,
3064 : "Invalid GeoTransform argument");
3065 0 : return CE_Failure;
3066 : }
3067 : }
3068 :
3069 : {
3070 358 : int iSource = 0;
3071 939 : for (const auto &osName : aosSourceNames)
3072 : {
3073 1162 : poExpression->RegisterVariable(osName,
3074 581 : &adfValuesForPixel[iSource++]);
3075 : }
3076 : }
3077 :
3078 358 : if (includeCenterCoords)
3079 : {
3080 1 : poExpression->RegisterVariable("_CENTER_X_", &dfCenterX);
3081 1 : poExpression->RegisterVariable("_CENTER_Y_", &dfCenterY);
3082 : }
3083 :
3084 358 : if (bHasNoData)
3085 : {
3086 10 : poExpression->RegisterVariable("NODATA", &dfNoData);
3087 : }
3088 :
3089 358 : if (strstr(pszExpression, "BANDS"))
3090 : {
3091 2 : poExpression->RegisterVector("BANDS", &adfValuesForPixel);
3092 : }
3093 :
3094 : std::unique_ptr<double, VSIFreeReleaser> padfResults(
3095 716 : static_cast<double *>(VSI_MALLOC2_VERBOSE(nXSize, sizeof(double))));
3096 358 : if (!padfResults)
3097 0 : return CE_Failure;
3098 :
3099 : /* ---- Set pixels ---- */
3100 358 : size_t ii = 0;
3101 5343 : for (int iLine = 0; iLine < nYSize; ++iLine)
3102 : {
3103 13007500 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
3104 : {
3105 13002500 : double &dfResult = padfResults.get()[iCol];
3106 13002500 : bool resultIsNoData = false;
3107 :
3108 39020600 : for (int iSrc = 0; iSrc < nSources; iSrc++)
3109 : {
3110 : // cppcheck-suppress unreadVariable
3111 26018100 : double dfVal = GetSrcVal(papoSources[iSrc], eSrcType, ii);
3112 :
3113 26018100 : if (bHasNoData && bPropagateNoData && IsNoData(dfVal, dfNoData))
3114 : {
3115 1 : resultIsNoData = true;
3116 : }
3117 :
3118 26018100 : adfValuesForPixel[iSrc] = dfVal;
3119 : }
3120 :
3121 13002500 : if (includeCenterCoords)
3122 : {
3123 : // Add 0.5 to pixel / line to move from pixel corner to cell center
3124 400 : gt.Apply(static_cast<double>(iCol + nXOff) + 0.5,
3125 400 : static_cast<double>(iLine + nYOff) + 0.5, &dfCenterX,
3126 : &dfCenterY);
3127 : }
3128 :
3129 13002500 : if (resultIsNoData)
3130 : {
3131 1 : dfResult = dfNoData;
3132 : }
3133 : else
3134 : {
3135 13002500 : if (auto eErr = poExpression->Evaluate(); eErr != CE_None)
3136 : {
3137 5 : return CE_Failure;
3138 : }
3139 :
3140 13002500 : dfResult = poExpression->Results()[0];
3141 : }
3142 : }
3143 :
3144 4985 : GDALCopyWords(padfResults.get(), GDT_Float64, sizeof(double),
3145 : static_cast<GByte *>(pData) +
3146 4985 : static_cast<GSpacing>(nLineSpace) * iLine,
3147 : eBufType, nPixelSpace, nXSize);
3148 : }
3149 :
3150 : /* ---- Return success ---- */
3151 353 : return CE_None;
3152 : } // ExprPixelFunc
3153 :
3154 : static constexpr char pszAreaPixelFuncMetadata[] =
3155 : "<PixelFunctionArgumentsList>"
3156 : " <Argument type='builtin' value='crs' />"
3157 : " <Argument type='builtin' value='xoff' />"
3158 : " <Argument type='builtin' value='yoff' />"
3159 : " <Argument type='builtin' value='geotransform' />"
3160 : "</PixelFunctionArgumentsList>";
3161 :
3162 4 : static CPLErr AreaPixelFunc(void ** /*papoSources*/, int nSources, void *pData,
3163 : int nXSize, int nYSize, GDALDataType /* eSrcType */,
3164 : GDALDataType eBufType, int nPixelSpace,
3165 : int nLineSpace, CSLConstList papszArgs)
3166 : {
3167 4 : if (nSources)
3168 : {
3169 1 : CPLError(CE_Failure, CPLE_AppDefined,
3170 : "area: unexpected source band(s)");
3171 1 : return CE_Failure;
3172 : }
3173 :
3174 3 : const char *pszGT = CSLFetchNameValue(papszArgs, "geotransform");
3175 3 : if (pszGT == nullptr)
3176 : {
3177 1 : CPLError(CE_Failure, CPLE_AppDefined,
3178 : "area: VRTDataset has no <GeoTransform>");
3179 1 : return CE_Failure;
3180 : }
3181 :
3182 2 : GDALGeoTransform gt;
3183 2 : if (!gt.Init(pszGT))
3184 : {
3185 0 : CPLError(CE_Failure, CPLE_AppDefined,
3186 : "area: Invalid GeoTransform argument");
3187 0 : return CE_Failure;
3188 : }
3189 :
3190 2 : const char *pszXOff = CSLFetchNameValue(papszArgs, "xoff");
3191 2 : const int nXOff = std::atoi(pszXOff);
3192 :
3193 2 : const char *pszYOff = CSLFetchNameValue(papszArgs, "yoff");
3194 2 : const int nYOff = std::atoi(pszYOff);
3195 :
3196 2 : const char *pszCrsPtr = CSLFetchNameValue(papszArgs, "crs");
3197 :
3198 2 : if (!pszCrsPtr)
3199 : {
3200 0 : CPLError(CE_Failure, CPLE_AppDefined, "area: VRTDataset has no <SRS>");
3201 0 : return CE_Failure;
3202 : }
3203 :
3204 2 : std::uintptr_t nCrsPtr = 0;
3205 2 : if (auto [end, ec] =
3206 2 : std::from_chars(pszCrsPtr, pszCrsPtr + strlen(pszCrsPtr), nCrsPtr);
3207 : ec != std::errc())
3208 : {
3209 : // Since "crs" is populated by GDAL, this should never happen.
3210 0 : CPLError(CE_Failure, CPLE_AppDefined, "Failed to read CRS");
3211 0 : return CE_Failure;
3212 : }
3213 :
3214 2 : const OGRSpatialReference *poCRS =
3215 : reinterpret_cast<const OGRSpatialReference *>(nCrsPtr);
3216 :
3217 2 : if (!poCRS)
3218 : {
3219 : // can't get here, but cppcheck doesn't know that
3220 0 : return CE_Failure;
3221 : }
3222 :
3223 2 : const OGRSpatialReference *poGeographicCRS = nullptr;
3224 2 : std::unique_ptr<OGRSpatialReference> poGeographicCRSHolder;
3225 2 : std::unique_ptr<OGRCoordinateTransformation> poTransform;
3226 :
3227 2 : if (!poCRS->IsGeographic())
3228 : {
3229 1 : poGeographicCRSHolder = std::make_unique<OGRSpatialReference>();
3230 1 : if (poGeographicCRSHolder->CopyGeogCSFrom(poCRS) != OGRERR_NONE)
3231 : {
3232 0 : CPLError(CE_Failure, CPLE_AppDefined,
3233 : "Cannot reproject geometry to geographic CRS");
3234 0 : return CE_Failure;
3235 : }
3236 1 : poGeographicCRSHolder->SetAxisMappingStrategy(
3237 : OAMS_TRADITIONAL_GIS_ORDER);
3238 :
3239 1 : poTransform.reset(OGRCreateCoordinateTransformation(
3240 1 : poCRS, poGeographicCRSHolder.get()));
3241 :
3242 1 : if (!poTransform)
3243 : {
3244 0 : CPLError(CE_Failure, CPLE_AppDefined,
3245 : "Cannot reproject geometry to geographic CRS");
3246 0 : return CE_Failure;
3247 : }
3248 :
3249 1 : poGeographicCRS = poGeographicCRSHolder.get();
3250 : }
3251 : else
3252 : {
3253 1 : poGeographicCRS = poCRS;
3254 : }
3255 :
3256 2 : geod_geodesic g{};
3257 2 : OGRErr eErr = OGRERR_NONE;
3258 2 : double dfSemiMajor = poGeographicCRS->GetSemiMajor(&eErr);
3259 2 : if (eErr != OGRERR_NONE)
3260 0 : return CE_Failure;
3261 2 : const double dfInvFlattening = poGeographicCRS->GetInvFlattening(&eErr);
3262 2 : if (eErr != OGRERR_NONE)
3263 0 : return CE_Failure;
3264 2 : geod_init(&g, dfSemiMajor,
3265 : dfInvFlattening != 0 ? 1.0 / dfInvFlattening : 0.0);
3266 :
3267 2 : std::array<double, 5> adfLon{};
3268 2 : std::array<double, 5> adfLat{};
3269 :
3270 22 : for (int iLine = 0; iLine < nYSize; ++iLine)
3271 : {
3272 420 : for (int iCol = 0; iCol < nXSize; ++iCol)
3273 : {
3274 400 : gt.Apply(static_cast<double>(iCol + nXOff),
3275 400 : static_cast<double>(iLine + nYOff), &adfLon[0],
3276 400 : &adfLat[0]);
3277 400 : gt.Apply(static_cast<double>(iCol + nXOff + 1),
3278 400 : static_cast<double>(iLine + nYOff), &adfLon[1],
3279 400 : &adfLat[1]);
3280 400 : gt.Apply(static_cast<double>(iCol + nXOff + 1),
3281 400 : static_cast<double>(iLine + nYOff + 1), &adfLon[2],
3282 400 : &adfLat[2]);
3283 400 : gt.Apply(static_cast<double>(iCol + nXOff),
3284 400 : static_cast<double>(iLine + nYOff + 1), &adfLon[3],
3285 400 : &adfLat[3]);
3286 400 : adfLon[4] = adfLon[0];
3287 400 : adfLat[4] = adfLat[0];
3288 :
3289 600 : if (poTransform &&
3290 200 : !poTransform->Transform(adfLon.size(), adfLon.data(),
3291 400 : adfLat.data(), nullptr))
3292 : {
3293 0 : CPLError(CE_Failure, CPLE_AppDefined,
3294 : "Failed to reproject cell corners to geographic CRS");
3295 0 : return CE_Failure;
3296 : }
3297 :
3298 400 : double dfArea = -1.0;
3299 400 : geod_polygonarea(&g, adfLat.data(), adfLon.data(),
3300 400 : static_cast<int>(adfLat.size()), &dfArea, nullptr);
3301 400 : dfArea = std::fabs(dfArea);
3302 :
3303 400 : GDALCopyWords(&dfArea, GDT_Float64, 0,
3304 : static_cast<GByte *>(pData) +
3305 400 : static_cast<GSpacing>(nLineSpace) * iLine +
3306 400 : iCol * nPixelSpace,
3307 : eBufType, nPixelSpace, 1);
3308 : }
3309 : }
3310 :
3311 2 : return CE_None;
3312 : } // AreaPixelFunc
3313 :
3314 : static const char pszReclassifyPixelFuncMetadata[] =
3315 : "<PixelFunctionArgumentsList>"
3316 : " <Argument name='mapping' "
3317 : " description='Lookup table for mapping, in format "
3318 : "from=to,from=to' "
3319 : " type='string'></Argument>"
3320 : " <Argument type='builtin' value='NoData' optional='true' />"
3321 : "</PixelFunctionArgumentsList>";
3322 :
3323 34 : static CPLErr ReclassifyPixelFunc(void **papoSources, int nSources, void *pData,
3324 : int nXSize, int nYSize, GDALDataType eSrcType,
3325 : GDALDataType eBufType, int nPixelSpace,
3326 : int nLineSpace, CSLConstList papszArgs)
3327 : {
3328 34 : if (GDALDataTypeIsComplex(eSrcType))
3329 : {
3330 0 : CPLError(CE_Failure, CPLE_AppDefined,
3331 : "reclassify cannot by applied to complex data types");
3332 0 : return CE_Failure;
3333 : }
3334 :
3335 34 : if (nSources != 1)
3336 : {
3337 0 : CPLError(CE_Failure, CPLE_AppDefined,
3338 : "reclassify only be applied to a single source at a time");
3339 0 : return CE_Failure;
3340 : }
3341 34 : std::optional<double> noDataValue{};
3342 :
3343 34 : const char *pszNoData = CSLFetchNameValue(papszArgs, "NoData");
3344 34 : if (pszNoData != nullptr)
3345 : {
3346 10 : noDataValue = CPLAtof(pszNoData);
3347 : }
3348 :
3349 34 : const char *pszMappings = CSLFetchNameValue(papszArgs, "mapping");
3350 34 : if (pszMappings == nullptr)
3351 : {
3352 0 : CPLError(CE_Failure, CPLE_AppDefined,
3353 : "reclassify must be called with 'mapping' argument");
3354 0 : return CE_Failure;
3355 : }
3356 :
3357 68 : gdal::Reclassifier oReclassifier;
3358 34 : if (auto eErr = oReclassifier.Init(pszMappings, noDataValue, eBufType);
3359 : eErr != CE_None)
3360 : {
3361 14 : return eErr;
3362 : }
3363 :
3364 : std::unique_ptr<double, VSIFreeReleaser> padfResults(
3365 40 : static_cast<double *>(VSI_MALLOC2_VERBOSE(nXSize, sizeof(double))));
3366 20 : if (!padfResults)
3367 0 : return CE_Failure;
3368 :
3369 20 : size_t ii = 0;
3370 20 : bool bSuccess = false;
3371 436 : for (int iLine = 0; iLine < nYSize; ++iLine)
3372 : {
3373 20808 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
3374 : {
3375 20392 : double srcVal = GetSrcVal(papoSources[0], eSrcType, ii);
3376 40784 : padfResults.get()[iCol] =
3377 20392 : oReclassifier.Reclassify(srcVal, bSuccess);
3378 20392 : if (!bSuccess)
3379 : {
3380 3 : CPLError(CE_Failure, CPLE_AppDefined,
3381 : "Encountered value %g with no specified mapping",
3382 : srcVal);
3383 3 : return CE_Failure;
3384 : }
3385 : }
3386 :
3387 416 : GDALCopyWords(padfResults.get(), GDT_Float64, sizeof(double),
3388 : static_cast<GByte *>(pData) +
3389 416 : static_cast<GSpacing>(nLineSpace) * iLine,
3390 : eBufType, nPixelSpace, nXSize);
3391 : }
3392 :
3393 17 : return CE_None;
3394 : } // ReclassifyPixelFunc
3395 :
3396 : struct MeanKernel
3397 : {
3398 : static constexpr const char *pszName = "mean";
3399 :
3400 : double dfMean = 0;
3401 : int nValidSources = 0;
3402 :
3403 568 : void Reset()
3404 : {
3405 568 : dfMean = 0;
3406 568 : nValidSources = 0;
3407 568 : }
3408 :
3409 94 : static CPLErr ProcessArguments(CSLConstList)
3410 : {
3411 94 : return CE_None;
3412 : }
3413 :
3414 1375 : void ProcessPixel(double dfVal)
3415 : {
3416 1375 : ++nValidSources;
3417 :
3418 1375 : if (CPL_UNLIKELY(std::isinf(dfVal)))
3419 : {
3420 620 : if (nValidSources == 1)
3421 : {
3422 310 : dfMean = dfVal;
3423 : }
3424 310 : else if (dfVal == -dfMean)
3425 : {
3426 62 : dfMean = std::numeric_limits<double>::quiet_NaN();
3427 : }
3428 : }
3429 755 : else if (CPL_UNLIKELY(std::isinf(dfMean)))
3430 : {
3431 186 : if (!std::isfinite(dfVal))
3432 : {
3433 62 : dfMean = std::numeric_limits<double>::quiet_NaN();
3434 : }
3435 : }
3436 : else
3437 : {
3438 569 : const double delta = dfVal - dfMean;
3439 569 : if (CPL_UNLIKELY(std::isinf(delta)))
3440 0 : dfMean += dfVal / nValidSources - dfMean / nValidSources;
3441 : else
3442 569 : dfMean += delta / nValidSources;
3443 : }
3444 1375 : }
3445 :
3446 566 : bool HasValue() const
3447 : {
3448 566 : return nValidSources > 0;
3449 : }
3450 :
3451 563 : double GetValue() const
3452 : {
3453 563 : return dfMean;
3454 : }
3455 : };
3456 :
3457 : struct GeoMeanKernel
3458 : {
3459 : static constexpr const char *pszName = "geometric_mean";
3460 :
3461 : double dfProduct = 1;
3462 : int nValidSources = 0;
3463 :
3464 6 : void Reset()
3465 : {
3466 6 : dfProduct = 1;
3467 6 : nValidSources = 0;
3468 6 : }
3469 :
3470 3 : static CPLErr ProcessArguments(CSLConstList)
3471 : {
3472 3 : return CE_None;
3473 : }
3474 :
3475 3 : void ProcessPixel(double dfVal)
3476 : {
3477 3 : dfProduct *= dfVal;
3478 3 : nValidSources++;
3479 3 : }
3480 :
3481 4 : bool HasValue() const
3482 : {
3483 4 : return nValidSources > 0;
3484 : }
3485 :
3486 1 : double GetValue() const
3487 : {
3488 1 : return std::pow(dfProduct, 1.0 / nValidSources);
3489 : }
3490 : };
3491 :
3492 : struct HarmonicMeanKernel
3493 : {
3494 : static constexpr const char *pszName = "harmonic_mean";
3495 :
3496 : double dfDenom = 0;
3497 : int nValidSources = 0;
3498 : bool bValueIsZero = false;
3499 : bool bPropagateZero = false;
3500 :
3501 10 : void Reset()
3502 : {
3503 10 : dfDenom = 0;
3504 10 : nValidSources = 0;
3505 10 : bValueIsZero = false;
3506 10 : }
3507 :
3508 7 : void ProcessPixel(double dfVal)
3509 : {
3510 7 : if (dfVal == 0)
3511 : {
3512 2 : bValueIsZero = true;
3513 : }
3514 : else
3515 : {
3516 5 : dfDenom += 1 / dfVal;
3517 : }
3518 7 : nValidSources++;
3519 7 : }
3520 :
3521 5 : CPLErr ProcessArguments(CSLConstList papszArgs)
3522 : {
3523 5 : bPropagateZero =
3524 5 : CPLTestBool(CSLFetchNameValueDef(papszArgs, "propagateZero", "0"));
3525 5 : return CE_None;
3526 : }
3527 :
3528 8 : bool HasValue() const
3529 : {
3530 8 : return dfDenom > 0 && (bPropagateZero || !bValueIsZero);
3531 : }
3532 :
3533 2 : double GetValue() const
3534 : {
3535 2 : if (bPropagateZero && bValueIsZero)
3536 : {
3537 1 : return 0;
3538 : }
3539 1 : return static_cast<double>(nValidSources) / dfDenom;
3540 : }
3541 : };
3542 :
3543 : struct MedianKernel
3544 : {
3545 : static constexpr const char *pszName = "median";
3546 :
3547 : mutable std::vector<double> values{};
3548 :
3549 9 : void Reset()
3550 : {
3551 9 : values.clear();
3552 9 : }
3553 :
3554 5 : static CPLErr ProcessArguments(CSLConstList)
3555 : {
3556 5 : return CE_None;
3557 : }
3558 :
3559 9 : void ProcessPixel(double dfVal)
3560 : {
3561 9 : if (!std::isnan(dfVal))
3562 : {
3563 9 : values.push_back(dfVal);
3564 : }
3565 9 : }
3566 :
3567 7 : bool HasValue() const
3568 : {
3569 7 : return !values.empty();
3570 : }
3571 :
3572 3 : double GetValue() const
3573 : {
3574 3 : std::sort(values.begin(), values.end());
3575 3 : if (values.size() % 2 == 0)
3576 : {
3577 : return 0.5 *
3578 1 : (values[values.size() / 2 - 1] + values[values.size() / 2]);
3579 : }
3580 :
3581 2 : return values[values.size() / 2];
3582 : }
3583 : };
3584 :
3585 : struct QuantileKernel
3586 : {
3587 : static constexpr const char *pszName = "quantile";
3588 :
3589 : mutable std::vector<double> values{};
3590 : double q = 0.5;
3591 :
3592 14 : void Reset()
3593 : {
3594 14 : values.clear();
3595 : // q intentionally preserved (set via ProcessArguments)
3596 14 : }
3597 :
3598 15 : CPLErr ProcessArguments(CSLConstList papszArgs)
3599 : {
3600 15 : const char *pszQ = CSLFetchNameValue(papszArgs, "q");
3601 15 : if (pszQ)
3602 : {
3603 : char *end;
3604 14 : const double dq = CPLStrtod(pszQ, &end);
3605 14 : while (isspace(*end))
3606 : {
3607 0 : end++;
3608 : }
3609 14 : if (*end != '\0' || dq < 0.0 || dq > 1.0 || std::isnan(dq))
3610 : {
3611 4 : CPLError(CE_Failure, CPLE_AppDefined,
3612 : "quantile: q must be between 0 and 1");
3613 4 : return CE_Failure;
3614 : }
3615 10 : q = dq;
3616 : }
3617 : else
3618 : {
3619 1 : CPLError(CE_Failure, CPLE_AppDefined,
3620 : "quantile: q must be specified");
3621 1 : return CE_Failure;
3622 : }
3623 10 : return CE_None;
3624 : }
3625 :
3626 15 : void ProcessPixel(double dfVal)
3627 : {
3628 15 : if (!std::isnan(dfVal))
3629 : {
3630 15 : values.push_back(dfVal);
3631 : }
3632 15 : }
3633 :
3634 12 : bool HasValue() const
3635 : {
3636 12 : return !values.empty();
3637 : }
3638 :
3639 7 : double GetValue() const
3640 : {
3641 7 : if (values.empty())
3642 : {
3643 0 : return std::numeric_limits<double>::quiet_NaN();
3644 : }
3645 :
3646 7 : std::sort(values.begin(), values.end());
3647 7 : const double loc = q * static_cast<double>(values.size() - 1);
3648 :
3649 : // Use formula from NumPy docs with default linear interpolation
3650 : // g: fractional component of loc
3651 : // j: integral component of loc
3652 : double j;
3653 7 : const double g = std::modf(loc, &j);
3654 :
3655 7 : if (static_cast<size_t>(j) + 1 == values.size())
3656 : {
3657 3 : return values[static_cast<size_t>(j)];
3658 : }
3659 :
3660 4 : return (1 - g) * values[static_cast<size_t>(j)] +
3661 4 : g * values[static_cast<size_t>(j) + 1];
3662 : }
3663 : };
3664 :
3665 : struct ModeKernel
3666 : {
3667 : static constexpr const char *pszName = "mode";
3668 :
3669 : std::map<double, size_t> counts{};
3670 : std::size_t nanCount{0};
3671 : double dfMax = std::numeric_limits<double>::quiet_NaN();
3672 : decltype(counts.begin()) oMax = counts.end();
3673 :
3674 7 : void Reset()
3675 : {
3676 7 : nanCount = 0;
3677 7 : counts.clear();
3678 7 : oMax = counts.end();
3679 7 : }
3680 :
3681 4 : static CPLErr ProcessArguments(CSLConstList)
3682 : {
3683 4 : return CE_None;
3684 : }
3685 :
3686 11 : void ProcessPixel(double dfVal)
3687 : {
3688 11 : if (std::isnan(dfVal))
3689 : {
3690 2 : nanCount += 1;
3691 2 : return;
3692 : }
3693 :
3694 : // if dfVal is NaN, try_emplace will return an entry for a different key!
3695 9 : auto [it, inserted] = counts.try_emplace(dfVal, 0);
3696 :
3697 9 : it->second += 1;
3698 :
3699 9 : if (oMax == counts.end() || it->second > oMax->second)
3700 : {
3701 5 : oMax = it;
3702 : }
3703 : }
3704 :
3705 5 : bool HasValue() const
3706 : {
3707 5 : return nanCount > 0 || oMax != counts.end();
3708 : }
3709 :
3710 3 : double GetValue() const
3711 : {
3712 3 : double ret = std::numeric_limits<double>::quiet_NaN();
3713 3 : if (oMax != counts.end())
3714 : {
3715 3 : const size_t nCount = oMax->second;
3716 3 : if (nCount > nanCount)
3717 2 : ret = oMax->first;
3718 : }
3719 3 : return ret;
3720 : }
3721 : };
3722 :
3723 : static const char pszBasicPixelFuncMetadata[] =
3724 : "<PixelFunctionArgumentsList>"
3725 : " <Argument type='builtin' value='NoData' optional='true' />"
3726 : " <Argument name='propagateNoData' description='Whether the output value "
3727 : "should be NoData as as soon as one source is NoData' type='boolean' "
3728 : "default='false' />"
3729 : "</PixelFunctionArgumentsList>";
3730 :
3731 : static const char pszQuantilePixelFuncMetadata[] =
3732 : "<PixelFunctionArgumentsList>"
3733 : " <Argument name='q' type='float' description='Quantile in [0,1]' />"
3734 : " <Argument type='builtin' value='NoData' optional='true' />"
3735 : " <Argument name='propagateNoData' type='boolean' default='false' />"
3736 : "</PixelFunctionArgumentsList>";
3737 :
3738 : #if defined(USE_SSE2) && !defined(USE_NEON_OPTIMIZATIONS)
3739 636 : inline __m128i packus_epi32(__m128i low, __m128i high)
3740 : {
3741 : #if __SSE4_1__
3742 : return _mm_packus_epi32(low, high); // Pack uint32 to uint16
3743 : #else
3744 1272 : low = _mm_add_epi32(low, _mm_set1_epi32(-32768));
3745 1272 : high = _mm_add_epi32(high, _mm_set1_epi32(-32768));
3746 1908 : return _mm_sub_epi16(_mm_packs_epi32(low, high), _mm_set1_epi16(-32768));
3747 : #endif
3748 : }
3749 : #endif
3750 :
3751 : #ifdef USE_SSE2
3752 :
3753 : template <class T, class SSEWrapper>
3754 46 : static void OptimizedMeanFloatSSE2(const void *const *papoSources, int nSources,
3755 : void *pData, int nXSize, int nYSize,
3756 : int nLineSpace)
3757 : {
3758 46 : assert(nSources >= 1);
3759 46 : constexpr int VALUES_PER_REG =
3760 : static_cast<int>(sizeof(typename SSEWrapper::Vec) / sizeof(T));
3761 46 : const T invSources = static_cast<T>(1.0) / static_cast<T>(nSources);
3762 46 : const auto invSourcesSSE = SSEWrapper::Set1(invSources);
3763 46 : const auto signMaskSSE = SSEWrapper::Set1(static_cast<T>(-0.0));
3764 46 : const auto infSSE = SSEWrapper::Set1(std::numeric_limits<T>::infinity());
3765 307 : for (int iLine = 0; iLine < nYSize; ++iLine)
3766 : {
3767 261 : T *CPL_RESTRICT pDest =
3768 : reinterpret_cast<T *>(static_cast<GByte *>(pData) +
3769 261 : static_cast<GSpacing>(nLineSpace) * iLine);
3770 261 : const size_t iOffsetLine = static_cast<size_t>(iLine) * nXSize;
3771 261 : int iCol = 0;
3772 2571 : for (; iCol < nXSize - (2 * VALUES_PER_REG - 1);
3773 : iCol += 2 * VALUES_PER_REG)
3774 : {
3775 2325 : auto reg0 = SSEWrapper::LoadU(
3776 2325 : static_cast<const T * CPL_RESTRICT>(papoSources[0]) +
3777 2325 : iOffsetLine + iCol);
3778 2325 : auto reg1 = SSEWrapper::LoadU(
3779 2325 : static_cast<const T * CPL_RESTRICT>(papoSources[0]) +
3780 2325 : iOffsetLine + iCol + VALUES_PER_REG);
3781 6700 : for (int iSrc = 1; iSrc < nSources; ++iSrc)
3782 : {
3783 4375 : const auto inputVal0 = SSEWrapper::LoadU(
3784 4375 : static_cast<const T * CPL_RESTRICT>(papoSources[iSrc]) +
3785 4375 : iOffsetLine + iCol);
3786 4375 : const auto inputVal1 = SSEWrapper::LoadU(
3787 4375 : static_cast<const T * CPL_RESTRICT>(papoSources[iSrc]) +
3788 4375 : iOffsetLine + iCol + VALUES_PER_REG);
3789 4375 : reg0 = SSEWrapper::Add(reg0, inputVal0);
3790 4375 : reg1 = SSEWrapper::Add(reg1, inputVal1);
3791 : }
3792 2325 : reg0 = SSEWrapper::Mul(reg0, invSourcesSSE);
3793 2325 : reg1 = SSEWrapper::Mul(reg1, invSourcesSSE);
3794 :
3795 : // Detect infinity that could happen when summing huge
3796 : // values
3797 2325 : if (SSEWrapper::MoveMask(SSEWrapper::Or(
3798 : SSEWrapper::CmpEq(SSEWrapper::AndNot(signMaskSSE, reg0),
3799 : infSSE),
3800 : SSEWrapper::CmpEq(SSEWrapper::AndNot(signMaskSSE, reg1),
3801 : infSSE))))
3802 : {
3803 15 : break;
3804 : }
3805 :
3806 2310 : SSEWrapper::StoreU(pDest + iCol, reg0);
3807 2310 : SSEWrapper::StoreU(pDest + iCol + VALUES_PER_REG, reg1);
3808 : }
3809 :
3810 : // Use numerically stable mean computation
3811 1240 : for (; iCol < nXSize; ++iCol)
3812 : {
3813 979 : T mean = static_cast<const T * CPL_RESTRICT>(
3814 979 : papoSources[0])[iOffsetLine + iCol];
3815 979 : if (nSources >= 2)
3816 : {
3817 969 : T new_val = static_cast<const T * CPL_RESTRICT>(
3818 969 : papoSources[1])[iOffsetLine + iCol];
3819 969 : if (CPL_UNLIKELY(std::isinf(new_val)))
3820 : {
3821 268 : if (new_val == -mean)
3822 : {
3823 10 : pDest[iCol] = std::numeric_limits<T>::quiet_NaN();
3824 10 : continue;
3825 : }
3826 : }
3827 701 : else if (CPL_UNLIKELY(std::isinf(mean)))
3828 : {
3829 144 : if (!std::isfinite(new_val))
3830 : {
3831 10 : pDest[iCol] = std::numeric_limits<T>::quiet_NaN();
3832 10 : continue;
3833 : }
3834 : }
3835 : else
3836 : {
3837 557 : const T delta = new_val - mean;
3838 557 : if (CPL_UNLIKELY(std::isinf(delta)))
3839 10 : mean += new_val * static_cast<T>(0.5) -
3840 10 : mean * static_cast<T>(0.5);
3841 : else
3842 547 : mean += delta * static_cast<T>(0.5);
3843 : }
3844 :
3845 1681 : for (int iSrc = 2; iSrc < nSources; ++iSrc)
3846 : {
3847 752 : new_val = static_cast<const T * CPL_RESTRICT>(
3848 752 : papoSources[iSrc])[iOffsetLine + iCol];
3849 752 : if (CPL_UNLIKELY(std::isinf(new_val)))
3850 : {
3851 196 : if (new_val == -mean)
3852 : {
3853 10 : mean = std::numeric_limits<T>::quiet_NaN();
3854 10 : break;
3855 : }
3856 : }
3857 556 : else if (CPL_UNLIKELY(std::isinf(mean)))
3858 : {
3859 72 : if (!std::isfinite(new_val))
3860 : {
3861 10 : mean = std::numeric_limits<T>::quiet_NaN();
3862 10 : break;
3863 : }
3864 : }
3865 : else
3866 : {
3867 484 : const T delta = new_val - mean;
3868 484 : if (CPL_UNLIKELY(std::isinf(delta)))
3869 62 : mean += new_val / static_cast<T>(iSrc + 1) -
3870 62 : mean / static_cast<T>(iSrc + 1);
3871 : else
3872 422 : mean += delta / static_cast<T>(iSrc + 1);
3873 : }
3874 : }
3875 : }
3876 959 : pDest[iCol] = mean;
3877 : }
3878 : }
3879 46 : }
3880 :
3881 : // clang-format off
3882 : namespace
3883 : {
3884 : #ifdef __AVX2__
3885 : struct SSEWrapperFloat
3886 : {
3887 : typedef __m256 Vec;
3888 :
3889 : static inline Vec Set1(float x) { return _mm256_set1_ps(x); }
3890 : static inline Vec LoadU(const float *x) { return _mm256_loadu_ps(x); }
3891 : static inline void StoreU(float *x, Vec y) { _mm256_storeu_ps(x, y); }
3892 : static inline Vec Add(Vec x, Vec y) { return _mm256_add_ps(x, y); }
3893 : static inline Vec Mul(Vec x, Vec y) { return _mm256_mul_ps(x, y); }
3894 : static inline Vec Or(Vec x, Vec y) { return _mm256_or_ps(x, y); }
3895 : static inline Vec AndNot(Vec x, Vec y) { return _mm256_andnot_ps(x, y); }
3896 : static inline Vec CmpEq(Vec x, Vec y) { return _mm256_cmp_ps(x, y, _CMP_EQ_OQ); }
3897 : static inline int MoveMask(Vec x) { return _mm256_movemask_ps(x); }
3898 : };
3899 :
3900 : struct SSEWrapperDouble
3901 : {
3902 : typedef __m256d Vec;
3903 :
3904 : static inline Vec Set1(double x) { return _mm256_set1_pd(x); }
3905 : static inline Vec LoadU(const double *x) { return _mm256_loadu_pd(x); }
3906 : static inline void StoreU(double *x, Vec y) { _mm256_storeu_pd(x, y); }
3907 : static inline Vec Add(Vec x, Vec y) { return _mm256_add_pd(x, y); }
3908 : static inline Vec Mul(Vec x, Vec y) { return _mm256_mul_pd(x, y); }
3909 : static inline Vec Or(Vec x, Vec y) { return _mm256_or_pd(x, y); }
3910 : static inline Vec AndNot(Vec x, Vec y) { return _mm256_andnot_pd(x, y); }
3911 : static inline Vec CmpEq(Vec x, Vec y) { return _mm256_cmp_pd(x, y, _CMP_EQ_OQ); }
3912 : static inline int MoveMask(Vec x) { return _mm256_movemask_pd(x); }
3913 : };
3914 :
3915 : #else
3916 :
3917 : struct SSEWrapperFloat
3918 : {
3919 : typedef __m128 Vec;
3920 :
3921 114 : static inline Vec Set1(float x) { return _mm_set1_ps(x); }
3922 3984 : static inline Vec LoadU(const float *x) { return _mm_loadu_ps(x); }
3923 666 : static inline void StoreU(float *x, Vec y) { _mm_storeu_ps(x, y); }
3924 2624 : static inline Vec Add(Vec x, Vec y) { return _mm_add_ps(x, y); }
3925 1360 : static inline Vec Mul(Vec x, Vec y) { return _mm_mul_ps(x, y); }
3926 680 : static inline Vec Or(Vec x, Vec y) { return _mm_or_ps(x, y); }
3927 1360 : static inline Vec AndNot(Vec x, Vec y) { return _mm_andnot_ps(x, y); }
3928 1360 : static inline Vec CmpEq(Vec x, Vec y) { return _mm_cmpeq_ps(x, y); }
3929 680 : static inline int MoveMask(Vec x) { return _mm_movemask_ps(x); }
3930 : };
3931 :
3932 : struct SSEWrapperDouble
3933 : {
3934 : typedef __m128d Vec;
3935 :
3936 162 : static inline Vec Set1(double x) { return _mm_set1_pd(x); }
3937 22816 : static inline Vec LoadU(const double *x) { return _mm_loadu_pd(x); }
3938 3954 : static inline void StoreU(double *x, Vec y) { _mm_storeu_pd(x, y); }
3939 14876 : static inline Vec Add(Vec x, Vec y) { return _mm_add_pd(x, y); }
3940 7940 : static inline Vec Mul(Vec x, Vec y) { return _mm_mul_pd(x, y); }
3941 3970 : static inline Vec Or(Vec x, Vec y) { return _mm_or_pd(x, y); }
3942 7940 : static inline Vec AndNot(Vec x, Vec y) { return _mm_andnot_pd(x, y); }
3943 7940 : static inline Vec CmpEq(Vec x, Vec y) { return _mm_cmpeq_pd(x, y); }
3944 3970 : static inline int MoveMask(Vec x) { return _mm_movemask_pd(x); }
3945 : };
3946 : #endif
3947 : } // namespace
3948 :
3949 : // clang-format on
3950 :
3951 : #endif // USE_SSE2
3952 :
3953 : template <typename Kernel>
3954 127 : static CPLErr BasicPixelFunc(void **papoSources, int nSources, void *pData,
3955 : int nXSize, int nYSize, GDALDataType eSrcType,
3956 : GDALDataType eBufType, int nPixelSpace,
3957 : int nLineSpace, CSLConstList papszArgs)
3958 : {
3959 : /* ---- Init ---- */
3960 152 : Kernel oKernel;
3961 :
3962 127 : if (GDALDataTypeIsComplex(eSrcType))
3963 : {
3964 1 : CPLError(CE_Failure, CPLE_AppDefined,
3965 : "Complex data types not supported by %s", oKernel.pszName);
3966 1 : return CE_Failure;
3967 : }
3968 :
3969 126 : double dfNoData{0};
3970 126 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
3971 126 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
3972 0 : return CE_Failure;
3973 :
3974 126 : const bool bPropagateNoData = CPLTestBool(
3975 : CSLFetchNameValueDef(papszArgs, "propagateNoData", "false"));
3976 :
3977 126 : if (oKernel.ProcessArguments(papszArgs) == CE_Failure)
3978 : {
3979 5 : return CE_Failure;
3980 : }
3981 :
3982 : #if defined(USE_SSE2) && !defined(USE_NEON_OPTIMIZATIONS)
3983 : if constexpr (std::is_same_v<Kernel, MeanKernel>)
3984 : {
3985 91 : if (!bHasNoData && eSrcType == GDT_UInt8 && eBufType == GDT_UInt8 &&
3986 185 : nPixelSpace == 1 &&
3987 : // We use signed int16 to accumulate
3988 11 : nSources <= std::numeric_limits<int16_t>::max() /
3989 11 : std::numeric_limits<uint8_t>::max())
3990 : {
3991 : using T = uint8_t;
3992 10 : constexpr int VALUES_PER_REG = 16;
3993 10 : if (nSources == 2)
3994 : {
3995 207 : for (int iLine = 0; iLine < nYSize; ++iLine)
3996 : {
3997 203 : T *CPL_RESTRICT pDest = reinterpret_cast<T *>(
3998 : static_cast<GByte *>(pData) +
3999 203 : static_cast<GSpacing>(nLineSpace) * iLine);
4000 203 : const size_t iOffsetLine =
4001 203 : static_cast<size_t>(iLine) * nXSize;
4002 203 : int iCol = 0;
4003 5209 : for (; iCol < nXSize - (VALUES_PER_REG - 1);
4004 : iCol += VALUES_PER_REG)
4005 : {
4006 : const __m128i inputVal0 =
4007 10012 : _mm_loadu_si128(reinterpret_cast<const __m128i *>(
4008 : static_cast<const T * CPL_RESTRICT>(
4009 : papoSources[0]) +
4010 5006 : iOffsetLine + iCol));
4011 : const __m128i inputVal1 =
4012 5006 : _mm_loadu_si128(reinterpret_cast<const __m128i *>(
4013 5006 : static_cast<const T * CPL_RESTRICT>(
4014 : papoSources[1]) +
4015 5006 : iOffsetLine + iCol));
4016 5006 : _mm_storeu_si128(
4017 5006 : reinterpret_cast<__m128i *>(pDest + iCol),
4018 : _mm_avg_epu8(inputVal0, inputVal1));
4019 : }
4020 235 : for (; iCol < nXSize; ++iCol)
4021 : {
4022 32 : uint32_t acc = 1 +
4023 32 : static_cast<const T * CPL_RESTRICT>(
4024 32 : papoSources[0])[iOffsetLine + iCol] +
4025 32 : static_cast<const T * CPL_RESTRICT>(
4026 32 : papoSources[1])[iOffsetLine + iCol];
4027 32 : pDest[iCol] = static_cast<T>(acc / 2);
4028 : }
4029 : }
4030 : }
4031 : else
4032 : {
4033 6 : libdivide::divider<uint16_t> fast_d(
4034 : static_cast<uint16_t>(nSources));
4035 : const auto halfConstant =
4036 6 : _mm_set1_epi16(static_cast<int16_t>(nSources / 2));
4037 211 : for (int iLine = 0; iLine < nYSize; ++iLine)
4038 : {
4039 205 : T *CPL_RESTRICT pDest =
4040 : static_cast<GByte *>(pData) +
4041 205 : static_cast<GSpacing>(nLineSpace) * iLine;
4042 205 : const size_t iOffsetLine =
4043 205 : static_cast<size_t>(iLine) * nXSize;
4044 205 : int iCol = 0;
4045 5220 : for (; iCol < nXSize - (VALUES_PER_REG - 1);
4046 : iCol += VALUES_PER_REG)
4047 : {
4048 5015 : __m128i reg0 = halfConstant;
4049 5015 : __m128i reg1 = halfConstant;
4050 20435 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
4051 : {
4052 15420 : const __m128i inputVal = _mm_loadu_si128(
4053 : reinterpret_cast<const __m128i *>(
4054 15420 : static_cast<const T * CPL_RESTRICT>(
4055 15420 : papoSources[iSrc]) +
4056 15420 : iOffsetLine + iCol));
4057 46260 : reg0 = _mm_add_epi16(
4058 : reg0, _mm_unpacklo_epi8(inputVal,
4059 : _mm_setzero_si128()));
4060 46260 : reg1 = _mm_add_epi16(
4061 : reg1, _mm_unpackhi_epi8(inputVal,
4062 : _mm_setzero_si128()));
4063 : }
4064 : reg0 /= fast_d;
4065 : reg1 /= fast_d;
4066 5015 : _mm_storeu_si128(
4067 5015 : reinterpret_cast<__m128i *>(pDest + iCol),
4068 : _mm_packus_epi16(reg0, reg1));
4069 : }
4070 280 : for (; iCol < nXSize; ++iCol)
4071 : {
4072 75 : uint32_t acc = nSources / 2;
4073 2175 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
4074 : {
4075 2100 : acc += static_cast<const T * CPL_RESTRICT>(
4076 2100 : papoSources[iSrc])[iOffsetLine + iCol];
4077 : }
4078 75 : pDest[iCol] = static_cast<T>(acc / nSources);
4079 : }
4080 : }
4081 : }
4082 10 : return CE_None;
4083 : }
4084 :
4085 81 : if (!bHasNoData && eSrcType == GDT_UInt8 && eBufType == GDT_UInt8 &&
4086 165 : nPixelSpace == 1 &&
4087 : // We use signed int32 to accumulate
4088 1 : nSources <= std::numeric_limits<int32_t>::max() /
4089 1 : std::numeric_limits<uint8_t>::max())
4090 : {
4091 : using T = uint8_t;
4092 1 : constexpr int VALUES_PER_REG = 16;
4093 1 : libdivide::divider<uint32_t> fast_d(nSources);
4094 1 : const auto halfConstant = _mm_set1_epi32(nSources / 2);
4095 2 : for (int iLine = 0; iLine < nYSize; ++iLine)
4096 : {
4097 1 : T *CPL_RESTRICT pDest =
4098 : static_cast<GByte *>(pData) +
4099 1 : static_cast<GSpacing>(nLineSpace) * iLine;
4100 1 : const size_t iOffsetLine = static_cast<size_t>(iLine) * nXSize;
4101 1 : int iCol = 0;
4102 4 : for (; iCol < nXSize - (VALUES_PER_REG - 1);
4103 : iCol += VALUES_PER_REG)
4104 : {
4105 3 : __m128i reg0 = halfConstant;
4106 3 : __m128i reg1 = halfConstant;
4107 3 : __m128i reg2 = halfConstant;
4108 3 : __m128i reg3 = halfConstant;
4109 98307 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
4110 : {
4111 : const __m128i inputVal =
4112 98304 : _mm_loadu_si128(reinterpret_cast<const __m128i *>(
4113 98304 : static_cast<const T * CPL_RESTRICT>(
4114 98304 : papoSources[iSrc]) +
4115 98304 : iOffsetLine + iCol));
4116 : const __m128i low =
4117 196608 : _mm_unpacklo_epi8(inputVal, _mm_setzero_si128());
4118 : const __m128i high =
4119 196608 : _mm_unpackhi_epi8(inputVal, _mm_setzero_si128());
4120 294912 : reg0 = _mm_add_epi32(
4121 : reg0, _mm_unpacklo_epi16(low, _mm_setzero_si128()));
4122 294912 : reg1 = _mm_add_epi32(
4123 : reg1, _mm_unpackhi_epi16(low, _mm_setzero_si128()));
4124 294912 : reg2 = _mm_add_epi32(
4125 : reg2,
4126 : _mm_unpacklo_epi16(high, _mm_setzero_si128()));
4127 294912 : reg3 = _mm_add_epi32(
4128 : reg3,
4129 : _mm_unpackhi_epi16(high, _mm_setzero_si128()));
4130 : }
4131 : reg0 /= fast_d;
4132 : reg1 /= fast_d;
4133 : reg2 /= fast_d;
4134 : reg3 /= fast_d;
4135 3 : _mm_storeu_si128(
4136 3 : reinterpret_cast<__m128i *>(pDest + iCol),
4137 : _mm_packus_epi16(packus_epi32(reg0, reg1),
4138 : packus_epi32(reg2, reg3)));
4139 : }
4140 16 : for (; iCol < nXSize; ++iCol)
4141 : {
4142 15 : uint32_t acc = nSources / 2;
4143 491535 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
4144 : {
4145 491520 : acc += static_cast<const T * CPL_RESTRICT>(
4146 491520 : papoSources[iSrc])[iOffsetLine + iCol];
4147 : }
4148 15 : pDest[iCol] = static_cast<T>(acc / nSources);
4149 : }
4150 : }
4151 1 : return CE_None;
4152 : }
4153 :
4154 80 : if (!bHasNoData && eSrcType == GDT_UInt16 && eBufType == GDT_UInt16 &&
4155 163 : nPixelSpace == 2 &&
4156 5 : nSources <= std::numeric_limits<int32_t>::max() /
4157 5 : std::numeric_limits<uint16_t>::max())
4158 : {
4159 5 : libdivide::divider<uint32_t> fast_d(nSources);
4160 : using T = uint16_t;
4161 5 : const auto halfConstant = _mm_set1_epi32(nSources / 2);
4162 5 : constexpr int VALUES_PER_REG = 8;
4163 59 : for (int iLine = 0; iLine < nYSize; ++iLine)
4164 : {
4165 54 : T *CPL_RESTRICT pDest = reinterpret_cast<T *>(
4166 : static_cast<GByte *>(pData) +
4167 54 : static_cast<GSpacing>(nLineSpace) * iLine);
4168 54 : const size_t iOffsetLine = static_cast<size_t>(iLine) * nXSize;
4169 54 : int iCol = 0;
4170 366 : for (; iCol < nXSize - (VALUES_PER_REG - 1);
4171 : iCol += VALUES_PER_REG)
4172 : {
4173 312 : __m128i reg0 = halfConstant;
4174 312 : __m128i reg1 = halfConstant;
4175 99534 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
4176 : {
4177 : const __m128i inputVal =
4178 99222 : _mm_loadu_si128(reinterpret_cast<const __m128i *>(
4179 99222 : static_cast<const T * CPL_RESTRICT>(
4180 99222 : papoSources[iSrc]) +
4181 99222 : iOffsetLine + iCol));
4182 297666 : reg0 = _mm_add_epi32(
4183 : reg0,
4184 : _mm_unpacklo_epi16(inputVal, _mm_setzero_si128()));
4185 297666 : reg1 = _mm_add_epi32(
4186 : reg1,
4187 : _mm_unpackhi_epi16(inputVal, _mm_setzero_si128()));
4188 : }
4189 : reg0 /= fast_d;
4190 : reg1 /= fast_d;
4191 312 : _mm_storeu_si128(reinterpret_cast<__m128i *>(pDest + iCol),
4192 : packus_epi32(reg0, reg1));
4193 : }
4194 182 : for (; iCol < nXSize; ++iCol)
4195 : {
4196 128 : uint32_t acc = nSources / 2;
4197 229846 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
4198 : {
4199 229718 : acc += static_cast<const T * CPL_RESTRICT>(
4200 229718 : papoSources[iSrc])[iOffsetLine + iCol];
4201 : }
4202 128 : pDest[iCol] = static_cast<T>(acc / nSources);
4203 : }
4204 : }
4205 5 : return CE_None;
4206 : }
4207 :
4208 75 : if (!bHasNoData && eSrcType == GDT_Int16 && eBufType == GDT_Int16 &&
4209 153 : nPixelSpace == 2 &&
4210 7 : nSources <= std::numeric_limits<int32_t>::max() /
4211 7 : std::numeric_limits<uint16_t>::max())
4212 : {
4213 7 : libdivide::divider<uint32_t> fast_d(nSources);
4214 : using T = int16_t;
4215 7 : const auto halfConstant = _mm_set1_epi32(nSources / 2);
4216 7 : const auto shift = _mm_set1_epi16(std::numeric_limits<T>::min());
4217 7 : constexpr int VALUES_PER_REG = 8;
4218 63 : for (int iLine = 0; iLine < nYSize; ++iLine)
4219 : {
4220 56 : T *CPL_RESTRICT pDest = reinterpret_cast<T *>(
4221 : static_cast<GByte *>(pData) +
4222 56 : static_cast<GSpacing>(nLineSpace) * iLine);
4223 56 : const size_t iOffsetLine = static_cast<size_t>(iLine) * nXSize;
4224 56 : int iCol = 0;
4225 374 : for (; iCol < nXSize - (VALUES_PER_REG - 1);
4226 : iCol += VALUES_PER_REG)
4227 : {
4228 318 : __m128i reg0 = halfConstant;
4229 318 : __m128i reg1 = halfConstant;
4230 99555 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
4231 : {
4232 : // Shift input values by 32768 to get unsigned values
4233 198474 : const __m128i inputVal = _mm_add_epi16(
4234 : _mm_loadu_si128(reinterpret_cast<const __m128i *>(
4235 99237 : static_cast<const T * CPL_RESTRICT>(
4236 99237 : papoSources[iSrc]) +
4237 99237 : iOffsetLine + iCol)),
4238 : shift);
4239 297711 : reg0 = _mm_add_epi32(
4240 : reg0,
4241 : _mm_unpacklo_epi16(inputVal, _mm_setzero_si128()));
4242 297711 : reg1 = _mm_add_epi32(
4243 : reg1,
4244 : _mm_unpackhi_epi16(inputVal, _mm_setzero_si128()));
4245 : }
4246 : reg0 /= fast_d;
4247 : reg1 /= fast_d;
4248 318 : _mm_storeu_si128(
4249 318 : reinterpret_cast<__m128i *>(pDest + iCol),
4250 : _mm_add_epi16(packus_epi32(reg0, reg1), shift));
4251 : }
4252 198 : for (; iCol < nXSize; ++iCol)
4253 : {
4254 142 : int32_t acc = (-std::numeric_limits<T>::min()) * nSources +
4255 142 : nSources / 2;
4256 229895 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
4257 : {
4258 229753 : acc += static_cast<const T * CPL_RESTRICT>(
4259 229753 : papoSources[iSrc])[iOffsetLine + iCol];
4260 : }
4261 284 : pDest[iCol] = static_cast<T>(acc / nSources +
4262 142 : std::numeric_limits<T>::min());
4263 : }
4264 : }
4265 7 : return CE_None;
4266 : }
4267 : }
4268 : #endif // defined(USE_SSE2) && !defined(USE_NEON_OPTIMIZATIONS)
4269 :
4270 : #if defined(USE_SSE2)
4271 : if constexpr (std::is_same_v<Kernel, MeanKernel>)
4272 : {
4273 71 : if (!bHasNoData && eSrcType == GDT_Float32 && eBufType == GDT_Float32 &&
4274 19 : nPixelSpace == 4 && nSources > 0)
4275 : {
4276 19 : OptimizedMeanFloatSSE2<float, SSEWrapperFloat>(
4277 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
4278 19 : return CE_None;
4279 : }
4280 :
4281 52 : if (!bHasNoData && eSrcType == GDT_Float64 && eBufType == GDT_Float64 &&
4282 27 : nPixelSpace == 8 && nSources > 0)
4283 : {
4284 27 : OptimizedMeanFloatSSE2<double, SSEWrapperDouble>(
4285 : papoSources, nSources, pData, nXSize, nYSize, nLineSpace);
4286 27 : return CE_None;
4287 : }
4288 : }
4289 : #endif // USE_SSE2
4290 :
4291 : /* ---- Set pixels ---- */
4292 52 : size_t ii = 0;
4293 104 : for (int iLine = 0; iLine < nYSize; ++iLine)
4294 : {
4295 666 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
4296 : {
4297 614 : oKernel.Reset();
4298 614 : bool bWriteNoData = false;
4299 :
4300 2107 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
4301 : {
4302 1505 : const double dfVal = GetSrcVal(papoSources[iSrc], eSrcType, ii);
4303 :
4304 1505 : if (bHasNoData && IsNoData(dfVal, dfNoData))
4305 : {
4306 85 : if (bPropagateNoData)
4307 : {
4308 12 : bWriteNoData = true;
4309 12 : break;
4310 : }
4311 : }
4312 : else
4313 : {
4314 1420 : oKernel.ProcessPixel(dfVal);
4315 : }
4316 : }
4317 :
4318 614 : double dfPixVal{dfNoData};
4319 614 : if (!bWriteNoData && oKernel.HasValue())
4320 : {
4321 579 : dfPixVal = oKernel.GetValue();
4322 : }
4323 :
4324 614 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
4325 : static_cast<GByte *>(pData) +
4326 614 : static_cast<GSpacing>(nLineSpace) * iLine +
4327 614 : iCol * nPixelSpace,
4328 : eBufType, nPixelSpace, 1);
4329 : }
4330 : }
4331 :
4332 : /* ---- Return success ---- */
4333 52 : return CE_None;
4334 : } // BasicPixelFunc
4335 :
4336 : /************************************************************************/
4337 : /* GDALRegisterDefaultPixelFunc() */
4338 : /************************************************************************/
4339 :
4340 : /**
4341 : * This adds a default set of pixel functions to the global list of
4342 : * available pixel functions for derived bands:
4343 : *
4344 : * - "real": extract real part from a single raster band (just a copy if the
4345 : * input is non-complex)
4346 : * - "imag": extract imaginary part from a single raster band (0 for
4347 : * non-complex)
4348 : * - "complex": make a complex band merging two bands used as real and
4349 : * imag values
4350 : * - "polar": make a complex band using input bands for amplitude and
4351 : * phase values (b1 * exp( j * b2 ))
4352 : * - "mod": extract module from a single raster band (real or complex)
4353 : * - "phase": extract phase from a single raster band [-PI,PI] (0 or PI for
4354 : non-complex)
4355 : * - "conj": computes the complex conjugate of a single raster band (just a
4356 : * copy if the input is non-complex)
4357 : * - "sum": sum 2 or more raster bands
4358 : * - "diff": computes the difference between 2 raster bands (b1 - b2)
4359 : * - "mul": multiply 2 or more raster bands
4360 : * - "div": divide one raster band by another (b1 / b2).
4361 : * - "min": minimum value of 2 or more raster bands
4362 : * - "max": maximum value of 2 or more raster bands
4363 : * - "norm_diff": computes the normalized difference between two raster bands:
4364 : * ``(b1 - b2)/(b1 + b2)``.
4365 : * - "cmul": multiply the first band for the complex conjugate of the second
4366 : * - "inv": inverse (1./x).
4367 : * - "intensity": computes the intensity Re(x*conj(x)) of a single raster band
4368 : * (real or complex)
4369 : * - "sqrt": perform the square root of a single raster band (real only)
4370 : * - "log10": compute the logarithm (base 10) of the abs of a single raster
4371 : * band (real or complex): log10( abs( x ) )
4372 : * - "dB": perform conversion to dB of the abs of a single raster
4373 : * band (real or complex): 20. * log10( abs( x ) ).
4374 : * Note: the optional fact parameter can be set to 10. to get the
4375 : * alternative formula: 10. * log10( abs( x ) )
4376 : * - "exp": computes the exponential of each element in the input band ``x``
4377 : * (of real values): ``e ^ x``.
4378 : * The function also accepts two optional parameters: ``base`` and
4379 : ``fact``
4380 : * that allow to compute the generalized formula: ``base ^ ( fact *
4381 : x)``.
4382 : * Note: this function is the recommended one to perform conversion
4383 : * form logarithmic scale (dB): `` 10. ^ (x / 20.)``, in this case
4384 : * ``base = 10.`` and ``fact = 1./20``
4385 : * - "dB2amp": perform scale conversion from logarithmic to linear
4386 : * (amplitude) (i.e. 10 ^ ( x / 20 ) ) of a single raster
4387 : * band (real only).
4388 : * Deprecated in GDAL v3.5. Please use the ``exp`` pixel function
4389 : with
4390 : * ``base = 10.`` and ``fact = 0.05`` i.e. ``1./20``
4391 : * - "dB2pow": perform scale conversion from logarithmic to linear
4392 : * (power) (i.e. 10 ^ ( x / 10 ) ) of a single raster
4393 : * band (real only)
4394 : * Deprecated in GDAL v3.5. Please use the ``exp`` pixel function
4395 : with
4396 : * ``base = 10.`` and ``fact = 0.1`` i.e. ``1./10``
4397 : * - "pow": raise a single raster band to a constant power
4398 : * - "interpolate_linear": interpolate values between two raster bands
4399 : * using linear interpolation
4400 : * - "interpolate_exp": interpolate values between two raster bands using
4401 : * exponential interpolation
4402 : * - "scale": Apply the RasterBand metadata values of "offset" and "scale"
4403 : * - "reclassify": Reclassify values matching ranges in a table
4404 : * - "nan": Convert incoming NoData values to IEEE 754 nan
4405 : *
4406 : * @see GDALAddDerivedBandPixelFunc
4407 : *
4408 : * @return CE_None
4409 : */
4410 1599 : CPLErr GDALRegisterDefaultPixelFunc()
4411 : {
4412 1599 : GDALAddDerivedBandPixelFunc("real", RealPixelFunc);
4413 1599 : GDALAddDerivedBandPixelFunc("imag", ImagPixelFunc);
4414 1599 : GDALAddDerivedBandPixelFunc("complex", ComplexPixelFunc);
4415 1599 : GDALAddDerivedBandPixelFuncWithArgs("polar", PolarPixelFunc,
4416 : pszPolarPixelFuncMetadata);
4417 1599 : GDALAddDerivedBandPixelFuncWithArgs("mod", ModulePixelFunc,
4418 : pszModulePixelFuncMetadata);
4419 1599 : GDALAddDerivedBandPixelFuncWithArgs("abs", ModulePixelFunc,
4420 : pszModulePixelFuncMetadata);
4421 1599 : GDALAddDerivedBandPixelFunc("phase", PhasePixelFunc);
4422 1599 : GDALAddDerivedBandPixelFunc("conj", ConjPixelFunc);
4423 1599 : GDALAddDerivedBandPixelFuncWithArgs("sum", SumPixelFunc,
4424 : pszSumPixelFuncMetadata);
4425 1599 : GDALAddDerivedBandPixelFuncWithArgs("diff", DiffPixelFunc,
4426 : pszDiffPixelFuncMetadata);
4427 1599 : GDALAddDerivedBandPixelFuncWithArgs("mul", MulPixelFunc,
4428 : pszMulPixelFuncMetadata);
4429 1599 : GDALAddDerivedBandPixelFuncWithArgs("div", DivPixelFunc,
4430 : pszDivPixelFuncMetadata);
4431 1599 : GDALAddDerivedBandPixelFunc("cmul", CMulPixelFunc);
4432 1599 : GDALAddDerivedBandPixelFuncWithArgs("inv", InvPixelFunc,
4433 : pszInvPixelFuncMetadata);
4434 1599 : GDALAddDerivedBandPixelFunc("intensity", IntensityPixelFunc);
4435 1599 : GDALAddDerivedBandPixelFuncWithArgs("sqrt", SqrtPixelFunc,
4436 : pszSqrtPixelFuncMetadata);
4437 1599 : GDALAddDerivedBandPixelFuncWithArgs("log10", Log10PixelFunc,
4438 : pszLog10PixelFuncMetadata);
4439 1599 : GDALAddDerivedBandPixelFuncWithArgs("dB", DBPixelFunc,
4440 : pszDBPixelFuncMetadata);
4441 1599 : GDALAddDerivedBandPixelFuncWithArgs("exp", ExpPixelFunc,
4442 : pszExpPixelFuncMetadata);
4443 1599 : GDALAddDerivedBandPixelFunc("dB2amp",
4444 : dB2AmpPixelFunc); // deprecated in v3.5
4445 1599 : GDALAddDerivedBandPixelFunc("dB2pow",
4446 : dB2PowPixelFunc); // deprecated in v3.5
4447 1599 : GDALAddDerivedBandPixelFuncWithArgs("pow", PowPixelFunc,
4448 : pszPowPixelFuncMetadata);
4449 1599 : GDALAddDerivedBandPixelFuncWithArgs("interpolate_linear",
4450 : InterpolatePixelFunc<InterpolateLinear>,
4451 : pszInterpolatePixelFuncMetadata);
4452 1599 : GDALAddDerivedBandPixelFuncWithArgs(
4453 : "interpolate_exp", InterpolatePixelFunc<InterpolateExponential>,
4454 : pszInterpolatePixelFuncMetadata);
4455 1599 : GDALAddDerivedBandPixelFuncWithArgs("replace_nodata",
4456 : ReplaceNoDataPixelFunc,
4457 : pszReplaceNoDataPixelFuncMetadata);
4458 1599 : GDALAddDerivedBandPixelFuncWithArgs("scale", ScalePixelFunc,
4459 : pszScalePixelFuncMetadata);
4460 1599 : GDALAddDerivedBandPixelFuncWithArgs("norm_diff", NormDiffPixelFunc,
4461 : pszNormDiffPixelFuncMetadata);
4462 1599 : GDALAddDerivedBandPixelFuncWithArgs("min", MinPixelFunc<ReturnValue>,
4463 : pszMinMaxFuncMetadataNodata);
4464 1599 : GDALAddDerivedBandPixelFuncWithArgs("argmin", MinPixelFunc<ReturnIndex>,
4465 : pszArgMinMaxFuncMetadataNodata);
4466 1599 : GDALAddDerivedBandPixelFuncWithArgs("max", MaxPixelFunc<ReturnValue>,
4467 : pszMinMaxFuncMetadataNodata);
4468 1599 : GDALAddDerivedBandPixelFuncWithArgs("argmax", MaxPixelFunc<ReturnIndex>,
4469 : pszArgMinMaxFuncMetadataNodata);
4470 1599 : GDALAddDerivedBandPixelFuncWithArgs("expression", ExprPixelFunc,
4471 : pszExprPixelFuncMetadata);
4472 1599 : GDALAddDerivedBandPixelFuncWithArgs("reclassify", ReclassifyPixelFunc,
4473 : pszReclassifyPixelFuncMetadata);
4474 1599 : GDALAddDerivedBandPixelFuncWithArgs("round", RoundPixelFunc,
4475 : pszRoundPixelFuncMetadata);
4476 1599 : GDALAddDerivedBandPixelFuncWithArgs("mean", BasicPixelFunc<MeanKernel>,
4477 : pszBasicPixelFuncMetadata);
4478 1599 : GDALAddDerivedBandPixelFuncWithArgs("geometric_mean",
4479 : BasicPixelFunc<GeoMeanKernel>,
4480 : pszBasicPixelFuncMetadata);
4481 1599 : GDALAddDerivedBandPixelFuncWithArgs("harmonic_mean",
4482 : BasicPixelFunc<HarmonicMeanKernel>,
4483 : pszBasicPixelFuncMetadata);
4484 1599 : GDALAddDerivedBandPixelFuncWithArgs("median", BasicPixelFunc<MedianKernel>,
4485 : pszBasicPixelFuncMetadata);
4486 1599 : GDALAddDerivedBandPixelFuncWithArgs("quantile",
4487 : BasicPixelFunc<QuantileKernel>,
4488 : pszQuantilePixelFuncMetadata);
4489 1599 : GDALAddDerivedBandPixelFuncWithArgs("mode", BasicPixelFunc<ModeKernel>,
4490 : pszBasicPixelFuncMetadata);
4491 1599 : GDALAddDerivedBandPixelFuncWithArgs("area", AreaPixelFunc,
4492 : pszAreaPixelFuncMetadata);
4493 1599 : return CE_None;
4494 : }
|
