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 <cmath>
15 : #include "gdal.h"
16 : #include "vrtdataset.h"
17 : #include "vrtexpression.h"
18 : #include "vrtreclassifier.h"
19 : #include "cpl_float.h"
20 :
21 : #if defined(__x86_64) || defined(_M_X64) || defined(USE_NEON_OPTIMIZATIONS)
22 : #define USE_SSE2
23 : #include "gdalsse_priv.h"
24 : #endif
25 :
26 : #include "gdal_priv_templates.hpp"
27 :
28 : #include <limits>
29 :
30 : template <typename T>
31 177300 : inline double GetSrcVal(const void *pSource, GDALDataType eSrcType, T ii)
32 : {
33 177300 : switch (eSrcType)
34 : {
35 0 : case GDT_Unknown:
36 0 : return 0;
37 125890 : case GDT_Byte:
38 125890 : return static_cast<const GByte *>(pSource)[ii];
39 476 : case GDT_Int8:
40 476 : return static_cast<const GInt8 *>(pSource)[ii];
41 1156 : case GDT_UInt16:
42 1156 : return static_cast<const GUInt16 *>(pSource)[ii];
43 1158 : case GDT_Int16:
44 1158 : return static_cast<const GInt16 *>(pSource)[ii];
45 8092 : case GDT_UInt32:
46 8092 : return static_cast<const GUInt32 *>(pSource)[ii];
47 2406 : case GDT_Int32:
48 2406 : return static_cast<const GInt32 *>(pSource)[ii];
49 : // Precision loss currently for int64/uint64
50 476 : case GDT_UInt64:
51 : return static_cast<double>(
52 476 : static_cast<const uint64_t *>(pSource)[ii]);
53 476 : case GDT_Int64:
54 : return static_cast<double>(
55 476 : static_cast<const int64_t *>(pSource)[ii]);
56 0 : case GDT_Float16:
57 0 : return static_cast<const GFloat16 *>(pSource)[ii];
58 4850 : case GDT_Float32:
59 4850 : return static_cast<const float *>(pSource)[ii];
60 22058 : case GDT_Float64:
61 22058 : return static_cast<const double *>(pSource)[ii];
62 1552 : case GDT_CInt16:
63 1552 : return static_cast<const GInt16 *>(pSource)[2 * ii];
64 952 : case GDT_CInt32:
65 952 : return static_cast<const GInt32 *>(pSource)[2 * ii];
66 0 : case GDT_CFloat16:
67 0 : return static_cast<const GFloat16 *>(pSource)[2 * ii];
68 5904 : case GDT_CFloat32:
69 5904 : return static_cast<const float *>(pSource)[2 * ii];
70 1854 : case GDT_CFloat64:
71 1854 : return static_cast<const double *>(pSource)[2 * ii];
72 0 : case GDT_TypeCount:
73 0 : break;
74 : }
75 0 : return 0;
76 : }
77 :
78 208 : static bool IsNoData(double dfVal, double dfNoData)
79 : {
80 208 : return dfVal == dfNoData || (std::isnan(dfVal) && std::isnan(dfNoData));
81 : }
82 :
83 1107 : static CPLErr FetchDoubleArg(CSLConstList papszArgs, const char *pszName,
84 : double *pdfX, double *pdfDefault = nullptr)
85 : {
86 1107 : const char *pszVal = CSLFetchNameValue(papszArgs, pszName);
87 :
88 1107 : if (pszVal == nullptr)
89 : {
90 70 : if (pdfDefault == nullptr)
91 : {
92 0 : CPLError(CE_Failure, CPLE_AppDefined,
93 : "Missing pixel function argument: %s", pszName);
94 0 : return CE_Failure;
95 : }
96 : else
97 : {
98 70 : *pdfX = *pdfDefault;
99 70 : return CE_None;
100 : }
101 : }
102 :
103 1037 : char *pszEnd = nullptr;
104 1037 : *pdfX = std::strtod(pszVal, &pszEnd);
105 1037 : if (pszEnd == pszVal)
106 : {
107 0 : CPLError(CE_Failure, CPLE_AppDefined,
108 : "Failed to parse pixel function argument: %s", pszName);
109 0 : return CE_Failure;
110 : }
111 :
112 1037 : return CE_None;
113 : }
114 :
115 7 : static CPLErr RealPixelFunc(void **papoSources, int nSources, void *pData,
116 : int nXSize, int nYSize, GDALDataType eSrcType,
117 : GDALDataType eBufType, int nPixelSpace,
118 : int nLineSpace)
119 : {
120 : /* ---- Init ---- */
121 7 : if (nSources != 1)
122 1 : return CE_Failure;
123 :
124 6 : const int nPixelSpaceSrc = GDALGetDataTypeSizeBytes(eSrcType);
125 6 : const size_t nLineSpaceSrc = static_cast<size_t>(nPixelSpaceSrc) * nXSize;
126 :
127 : /* ---- Set pixels ---- */
128 98 : for (int iLine = 0; iLine < nYSize; ++iLine)
129 : {
130 92 : GDALCopyWords(static_cast<GByte *>(papoSources[0]) +
131 92 : nLineSpaceSrc * iLine,
132 : eSrcType, nPixelSpaceSrc,
133 : static_cast<GByte *>(pData) +
134 92 : static_cast<GSpacing>(nLineSpace) * iLine,
135 : eBufType, nPixelSpace, nXSize);
136 : }
137 :
138 : /* ---- Return success ---- */
139 6 : return CE_None;
140 : } // RealPixelFunc
141 :
142 8 : static CPLErr ImagPixelFunc(void **papoSources, int nSources, void *pData,
143 : int nXSize, int nYSize, GDALDataType eSrcType,
144 : GDALDataType eBufType, int nPixelSpace,
145 : int nLineSpace)
146 : {
147 : /* ---- Init ---- */
148 8 : if (nSources != 1)
149 1 : return CE_Failure;
150 :
151 7 : if (GDALDataTypeIsComplex(eSrcType))
152 : {
153 6 : const GDALDataType eSrcBaseType = GDALGetNonComplexDataType(eSrcType);
154 6 : const int nPixelSpaceSrc = GDALGetDataTypeSizeBytes(eSrcType);
155 6 : const size_t nLineSpaceSrc =
156 6 : static_cast<size_t>(nPixelSpaceSrc) * nXSize;
157 :
158 6 : const void *const pImag = static_cast<GByte *>(papoSources[0]) +
159 6 : GDALGetDataTypeSizeBytes(eSrcType) / 2;
160 :
161 : /* ---- Set pixels ---- */
162 56 : for (int iLine = 0; iLine < nYSize; ++iLine)
163 : {
164 50 : GDALCopyWords(static_cast<const GByte *>(pImag) +
165 50 : nLineSpaceSrc * iLine,
166 : eSrcBaseType, nPixelSpaceSrc,
167 : static_cast<GByte *>(pData) +
168 50 : static_cast<GSpacing>(nLineSpace) * iLine,
169 : eBufType, nPixelSpace, nXSize);
170 : }
171 : }
172 : else
173 : {
174 1 : const double dfImag = 0;
175 :
176 : /* ---- Set pixels ---- */
177 21 : for (int iLine = 0; iLine < nYSize; ++iLine)
178 : {
179 : // Always copy from the same location.
180 20 : GDALCopyWords(&dfImag, eSrcType, 0,
181 : static_cast<GByte *>(pData) +
182 20 : static_cast<GSpacing>(nLineSpace) * iLine,
183 : eBufType, nPixelSpace, nXSize);
184 : }
185 : }
186 :
187 : /* ---- Return success ---- */
188 7 : return CE_None;
189 : } // ImagPixelFunc
190 :
191 6 : static CPLErr ComplexPixelFunc(void **papoSources, int nSources, void *pData,
192 : int nXSize, int nYSize, GDALDataType eSrcType,
193 : GDALDataType eBufType, int nPixelSpace,
194 : int nLineSpace)
195 : {
196 : /* ---- Init ---- */
197 6 : if (nSources != 2)
198 1 : return CE_Failure;
199 :
200 5 : const void *const pReal = papoSources[0];
201 5 : const void *const pImag = papoSources[1];
202 :
203 : /* ---- Set pixels ---- */
204 5 : size_t ii = 0;
205 281 : for (int iLine = 0; iLine < nYSize; ++iLine)
206 : {
207 17060 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
208 : {
209 : // Source raster pixels may be obtained with GetSrcVal macro.
210 : const double adfPixVal[2] = {
211 16784 : GetSrcVal(pReal, eSrcType, ii), // re
212 33568 : GetSrcVal(pImag, eSrcType, ii) // im
213 16784 : };
214 :
215 16784 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
216 : static_cast<GByte *>(pData) +
217 16784 : static_cast<GSpacing>(nLineSpace) * iLine +
218 16784 : iCol * nPixelSpace,
219 : eBufType, nPixelSpace, 1);
220 : }
221 : }
222 :
223 : /* ---- Return success ---- */
224 5 : return CE_None;
225 : } // ComplexPixelFunc
226 :
227 : typedef enum
228 : {
229 : GAT_amplitude,
230 : GAT_intensity,
231 : GAT_dB
232 : } PolarAmplitudeType;
233 :
234 : static const char pszPolarPixelFuncMetadata[] =
235 : "<PixelFunctionArgumentsList>"
236 : " <Argument name='amplitude_type' description='Amplitude Type' "
237 : "type='string-select' default='AMPLITUDE'>"
238 : " <Value>INTENSITY</Value>"
239 : " <Value>dB</Value>"
240 : " <Value>AMPLITUDE</Value>"
241 : " </Argument>"
242 : "</PixelFunctionArgumentsList>";
243 :
244 4 : static CPLErr PolarPixelFunc(void **papoSources, int nSources, void *pData,
245 : int nXSize, int nYSize, GDALDataType eSrcType,
246 : GDALDataType eBufType, int nPixelSpace,
247 : int nLineSpace, CSLConstList papszArgs)
248 : {
249 : /* ---- Init ---- */
250 4 : if (nSources != 2)
251 0 : return CE_Failure;
252 :
253 4 : const char pszName[] = "amplitude_type";
254 4 : const char *pszVal = CSLFetchNameValue(papszArgs, pszName);
255 4 : PolarAmplitudeType amplitudeType = GAT_amplitude;
256 4 : if (pszVal != nullptr)
257 : {
258 3 : if (strcmp(pszVal, "INTENSITY") == 0)
259 1 : amplitudeType = GAT_intensity;
260 2 : else if (strcmp(pszVal, "dB") == 0)
261 1 : amplitudeType = GAT_dB;
262 1 : else if (strcmp(pszVal, "AMPLITUDE") != 0)
263 : {
264 0 : CPLError(CE_Failure, CPLE_AppDefined,
265 : "Invalid value for pixel function argument '%s': %s",
266 : pszName, pszVal);
267 0 : return CE_Failure;
268 : }
269 : }
270 :
271 4 : const void *const pAmp = papoSources[0];
272 4 : const void *const pPhase = papoSources[1];
273 :
274 : /* ---- Set pixels ---- */
275 4 : size_t ii = 0;
276 84 : for (int iLine = 0; iLine < nYSize; ++iLine)
277 : {
278 1680 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
279 : {
280 : // Source raster pixels may be obtained with GetSrcVal macro.
281 1600 : double dfAmp = GetSrcVal(pAmp, eSrcType, ii);
282 1600 : switch (amplitudeType)
283 : {
284 400 : case GAT_intensity:
285 : // clip to zero
286 400 : dfAmp = dfAmp <= 0 ? 0 : std::sqrt(dfAmp);
287 400 : break;
288 400 : case GAT_dB:
289 400 : dfAmp = dfAmp <= 0
290 400 : ? -std::numeric_limits<double>::infinity()
291 400 : : pow(10, dfAmp / 20.);
292 400 : break;
293 800 : case GAT_amplitude:
294 800 : break;
295 : }
296 1600 : const double dfPhase = GetSrcVal(pPhase, eSrcType, ii);
297 : const double adfPixVal[2] = {
298 1600 : dfAmp * std::cos(dfPhase), // re
299 1600 : dfAmp * std::sin(dfPhase) // im
300 1600 : };
301 :
302 1600 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
303 : static_cast<GByte *>(pData) +
304 1600 : static_cast<GSpacing>(nLineSpace) * iLine +
305 1600 : iCol * nPixelSpace,
306 : eBufType, nPixelSpace, 1);
307 : }
308 : }
309 :
310 : /* ---- Return success ---- */
311 4 : return CE_None;
312 : } // PolarPixelFunc
313 :
314 4 : static CPLErr ModulePixelFunc(void **papoSources, int nSources, void *pData,
315 : int nXSize, int nYSize, GDALDataType eSrcType,
316 : GDALDataType eBufType, int nPixelSpace,
317 : int nLineSpace)
318 : {
319 : /* ---- Init ---- */
320 4 : if (nSources != 1)
321 1 : return CE_Failure;
322 :
323 3 : if (GDALDataTypeIsComplex(eSrcType))
324 : {
325 2 : const void *pReal = papoSources[0];
326 2 : const void *pImag = static_cast<GByte *>(papoSources[0]) +
327 2 : GDALGetDataTypeSizeBytes(eSrcType) / 2;
328 :
329 : /* ---- Set pixels ---- */
330 2 : size_t ii = 0;
331 14 : for (int iLine = 0; iLine < nYSize; ++iLine)
332 : {
333 72 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
334 : {
335 : // Source raster pixels may be obtained with GetSrcVal macro.
336 60 : const double dfReal = GetSrcVal(pReal, eSrcType, ii);
337 60 : const double dfImag = GetSrcVal(pImag, eSrcType, ii);
338 :
339 60 : const double dfPixVal = sqrt(dfReal * dfReal + dfImag * dfImag);
340 :
341 60 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
342 : static_cast<GByte *>(pData) +
343 60 : static_cast<GSpacing>(nLineSpace) * iLine +
344 60 : iCol * nPixelSpace,
345 : eBufType, nPixelSpace, 1);
346 : }
347 : }
348 : }
349 : else
350 : {
351 : /* ---- Set pixels ---- */
352 1 : size_t ii = 0;
353 21 : for (int iLine = 0; iLine < nYSize; ++iLine)
354 : {
355 420 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
356 : {
357 : // Source raster pixels may be obtained with GetSrcVal macro.
358 : const double dfPixVal =
359 400 : fabs(GetSrcVal(papoSources[0], eSrcType, ii));
360 :
361 400 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
362 : static_cast<GByte *>(pData) +
363 400 : static_cast<GSpacing>(nLineSpace) * iLine +
364 400 : iCol * nPixelSpace,
365 : eBufType, nPixelSpace, 1);
366 : }
367 : }
368 : }
369 :
370 : /* ---- Return success ---- */
371 3 : return CE_None;
372 : } // ModulePixelFunc
373 :
374 5 : static CPLErr PhasePixelFunc(void **papoSources, int nSources, void *pData,
375 : int nXSize, int nYSize, GDALDataType eSrcType,
376 : GDALDataType eBufType, int nPixelSpace,
377 : int nLineSpace)
378 : {
379 : /* ---- Init ---- */
380 5 : if (nSources != 1)
381 1 : return CE_Failure;
382 :
383 4 : if (GDALDataTypeIsComplex(eSrcType))
384 : {
385 2 : const void *const pReal = papoSources[0];
386 2 : const void *const pImag = static_cast<GByte *>(papoSources[0]) +
387 2 : GDALGetDataTypeSizeBytes(eSrcType) / 2;
388 :
389 : /* ---- Set pixels ---- */
390 2 : size_t ii = 0;
391 14 : for (int iLine = 0; iLine < nYSize; ++iLine)
392 : {
393 72 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
394 : {
395 : // Source raster pixels may be obtained with GetSrcVal macro.
396 60 : const double dfReal = GetSrcVal(pReal, eSrcType, ii);
397 60 : const double dfImag = GetSrcVal(pImag, eSrcType, ii);
398 :
399 60 : const double dfPixVal = atan2(dfImag, dfReal);
400 :
401 60 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
402 : static_cast<GByte *>(pData) +
403 60 : static_cast<GSpacing>(nLineSpace) * iLine +
404 60 : iCol * nPixelSpace,
405 : eBufType, nPixelSpace, 1);
406 : }
407 : }
408 : }
409 2 : else if (GDALDataTypeIsInteger(eSrcType) && !GDALDataTypeIsSigned(eSrcType))
410 : {
411 1 : constexpr double dfZero = 0;
412 7 : for (int iLine = 0; iLine < nYSize; ++iLine)
413 : {
414 6 : GDALCopyWords(&dfZero, GDT_Float64, 0,
415 : static_cast<GByte *>(pData) +
416 6 : static_cast<GSpacing>(nLineSpace) * iLine,
417 : eBufType, nPixelSpace, nXSize);
418 : }
419 : }
420 : else
421 : {
422 : /* ---- Set pixels ---- */
423 1 : size_t ii = 0;
424 7 : for (int iLine = 0; iLine < nYSize; ++iLine)
425 : {
426 36 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
427 : {
428 30 : const void *const pReal = papoSources[0];
429 :
430 : // Source raster pixels may be obtained with GetSrcVal macro.
431 30 : const double dfReal = GetSrcVal(pReal, eSrcType, ii);
432 30 : const double dfPixVal = (dfReal < 0) ? M_PI : 0.0;
433 :
434 30 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
435 : static_cast<GByte *>(pData) +
436 30 : static_cast<GSpacing>(nLineSpace) * iLine +
437 30 : iCol * nPixelSpace,
438 : eBufType, nPixelSpace, 1);
439 : }
440 : }
441 : }
442 :
443 : /* ---- Return success ---- */
444 4 : return CE_None;
445 : } // PhasePixelFunc
446 :
447 4 : static CPLErr ConjPixelFunc(void **papoSources, int nSources, void *pData,
448 : int nXSize, int nYSize, GDALDataType eSrcType,
449 : GDALDataType eBufType, int nPixelSpace,
450 : int nLineSpace)
451 : {
452 : /* ---- Init ---- */
453 4 : if (nSources != 1)
454 1 : return CE_Failure;
455 :
456 3 : if (GDALDataTypeIsComplex(eSrcType) && GDALDataTypeIsComplex(eBufType))
457 : {
458 2 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
459 2 : const void *const pReal = papoSources[0];
460 2 : const void *const pImag =
461 2 : static_cast<GByte *>(papoSources[0]) + nOffset;
462 :
463 : /* ---- Set pixels ---- */
464 2 : size_t ii = 0;
465 14 : for (int iLine = 0; iLine < nYSize; ++iLine)
466 : {
467 72 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
468 : {
469 : // Source raster pixels may be obtained with GetSrcVal macro.
470 : const double adfPixVal[2] = {
471 60 : +GetSrcVal(pReal, eSrcType, ii), // re
472 120 : -GetSrcVal(pImag, eSrcType, ii) // im
473 60 : };
474 :
475 60 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
476 : static_cast<GByte *>(pData) +
477 60 : static_cast<GSpacing>(nLineSpace) * iLine +
478 60 : iCol * nPixelSpace,
479 : eBufType, nPixelSpace, 1);
480 : }
481 : }
482 : }
483 : else
484 : {
485 : // No complex data type.
486 1 : return RealPixelFunc(papoSources, nSources, pData, nXSize, nYSize,
487 1 : eSrcType, eBufType, nPixelSpace, nLineSpace);
488 : }
489 :
490 : /* ---- Return success ---- */
491 2 : return CE_None;
492 : } // ConjPixelFunc
493 :
494 : #ifdef USE_SSE2
495 :
496 : /************************************************************************/
497 : /* OptimizedSumToFloat_SSE2() */
498 : /************************************************************************/
499 :
500 : template <typename Tsrc>
501 87 : static void OptimizedSumToFloat_SSE2(double dfK, void *pOutBuffer,
502 : int nLineSpace, int nXSize, int nYSize,
503 : int nSources,
504 : const void *const *papoSources)
505 : {
506 87 : const XMMReg4Float cst = XMMReg4Float::Set1(static_cast<float>(dfK));
507 :
508 279 : for (int iLine = 0; iLine < nYSize; ++iLine)
509 : {
510 192 : float *CPL_RESTRICT const pDest = reinterpret_cast<float *>(
511 : static_cast<GByte *>(pOutBuffer) +
512 192 : static_cast<GSpacing>(nLineSpace) * iLine);
513 192 : const size_t iOffsetLine = static_cast<size_t>(iLine) * nXSize;
514 :
515 192 : constexpr int VALUES_PER_REG = 4;
516 192 : constexpr int UNROLLING = 4 * VALUES_PER_REG;
517 192 : int iCol = 0;
518 384 : for (; iCol < nXSize - (UNROLLING - 1); iCol += UNROLLING)
519 : {
520 192 : XMMReg4Float d0(cst);
521 192 : XMMReg4Float d1(cst);
522 192 : XMMReg4Float d2(cst);
523 192 : XMMReg4Float d3(cst);
524 556 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
525 : {
526 : XMMReg4Float t0, t1, t2, t3;
527 364 : XMMReg4Float::Load16Val(
528 364 : static_cast<const Tsrc * CPL_RESTRICT>(papoSources[iSrc]) +
529 364 : iOffsetLine + iCol,
530 : t0, t1, t2, t3);
531 364 : d0 += t0;
532 364 : d1 += t1;
533 364 : d2 += t2;
534 364 : d3 += t3;
535 : }
536 192 : d0.Store4Val(pDest + iCol + VALUES_PER_REG * 0);
537 192 : d1.Store4Val(pDest + iCol + VALUES_PER_REG * 1);
538 192 : d2.Store4Val(pDest + iCol + VALUES_PER_REG * 2);
539 192 : d3.Store4Val(pDest + iCol + VALUES_PER_REG * 3);
540 : }
541 :
542 444 : for (; iCol < nXSize; iCol++)
543 : {
544 252 : float d = static_cast<float>(dfK);
545 676 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
546 : {
547 424 : d += static_cast<const Tsrc * CPL_RESTRICT>(
548 424 : papoSources[iSrc])[iOffsetLine + iCol];
549 : }
550 252 : pDest[iCol] = d;
551 : }
552 : }
553 87 : }
554 :
555 : /************************************************************************/
556 : /* OptimizedSumToDouble_SSE2() */
557 : /************************************************************************/
558 :
559 : template <typename Tsrc>
560 141 : static void OptimizedSumToDouble_SSE2(double dfK, void *pOutBuffer,
561 : int nLineSpace, int nXSize, int nYSize,
562 : int nSources,
563 : const void *const *papoSources)
564 : {
565 141 : const XMMReg4Double cst = XMMReg4Double::Set1(dfK);
566 :
567 4041 : for (int iLine = 0; iLine < nYSize; ++iLine)
568 : {
569 3900 : double *CPL_RESTRICT const pDest = reinterpret_cast<double *>(
570 : static_cast<GByte *>(pOutBuffer) +
571 3900 : static_cast<GSpacing>(nLineSpace) * iLine);
572 3900 : const size_t iOffsetLine = static_cast<size_t>(iLine) * nXSize;
573 :
574 3900 : constexpr int VALUES_PER_REG = 4;
575 3900 : constexpr int UNROLLING = 2 * VALUES_PER_REG;
576 3900 : int iCol = 0;
577 1624476 : for (; iCol < nXSize - (UNROLLING - 1); iCol += UNROLLING)
578 : {
579 1620574 : XMMReg4Double d0(cst);
580 1620574 : XMMReg4Double d1(cst);
581 4861552 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
582 : {
583 : XMMReg4Double t0, t1;
584 3240988 : XMMReg4Double::Load8Val(
585 3240988 : static_cast<const Tsrc * CPL_RESTRICT>(papoSources[iSrc]) +
586 3240988 : iOffsetLine + iCol,
587 : t0, t1);
588 3240988 : d0 += t0;
589 3240988 : d1 += t1;
590 : }
591 1620574 : d0.Store4Val(pDest + iCol + VALUES_PER_REG * 0);
592 1620574 : d1.Store4Val(pDest + iCol + VALUES_PER_REG * 1);
593 : }
594 :
595 4478 : for (; iCol < nXSize; iCol++)
596 : {
597 578 : double d = dfK;
598 1450 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
599 : {
600 872 : d += static_cast<const Tsrc * CPL_RESTRICT>(
601 872 : papoSources[iSrc])[iOffsetLine + iCol];
602 : }
603 578 : pDest[iCol] = d;
604 : }
605 : }
606 141 : }
607 :
608 : #endif // USE_SSE2
609 :
610 : /************************************************************************/
611 : /* OptimizedSumPackedOutput() */
612 : /************************************************************************/
613 :
614 : template <typename Tsrc, typename Tdest>
615 264 : static void OptimizedSumPackedOutput(double dfK, void *pOutBuffer,
616 : int nLineSpace, int nXSize, int nYSize,
617 : int nSources,
618 : const void *const *papoSources)
619 : {
620 : #ifdef USE_SSE2
621 : if constexpr (std::is_same_v<Tdest, float> && !std::is_same_v<Tsrc, double>)
622 : {
623 87 : OptimizedSumToFloat_SSE2<Tsrc>(dfK, pOutBuffer, nLineSpace, nXSize,
624 : nYSize, nSources, papoSources);
625 : }
626 : else if constexpr (std::is_same_v<Tdest, double>)
627 : {
628 141 : OptimizedSumToDouble_SSE2<Tsrc>(dfK, pOutBuffer, nLineSpace, nXSize,
629 : nYSize, nSources, papoSources);
630 : }
631 : else
632 : #endif // USE_SSE2
633 : {
634 36 : const Tdest nCst = static_cast<Tdest>(dfK);
635 144 : for (int iLine = 0; iLine < nYSize; ++iLine)
636 : {
637 108 : Tdest *CPL_RESTRICT const pDest = reinterpret_cast<Tdest *>(
638 : static_cast<GByte *>(pOutBuffer) +
639 108 : static_cast<GSpacing>(nLineSpace) * iLine);
640 108 : const size_t iOffsetLine = static_cast<size_t>(iLine) * nXSize;
641 :
642 : #define LOAD_SRCVAL(iSrc_, j_) \
643 : static_cast<Tdest>(static_cast<const Tsrc * CPL_RESTRICT>( \
644 : papoSources[(iSrc_)])[iOffsetLine + iCol + (j_)])
645 :
646 108 : constexpr int UNROLLING = 4;
647 108 : int iCol = 0;
648 580 : for (; iCol < nXSize - (UNROLLING - 1); iCol += UNROLLING)
649 : {
650 472 : Tdest d[4] = {nCst, nCst, nCst, nCst};
651 1616 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
652 : {
653 1144 : d[0] += LOAD_SRCVAL(iSrc, 0);
654 1144 : d[1] += LOAD_SRCVAL(iSrc, 1);
655 1144 : d[2] += LOAD_SRCVAL(iSrc, 2);
656 1144 : d[3] += LOAD_SRCVAL(iSrc, 3);
657 : }
658 472 : pDest[iCol + 0] = d[0];
659 472 : pDest[iCol + 1] = d[1];
660 472 : pDest[iCol + 2] = d[2];
661 472 : pDest[iCol + 3] = d[3];
662 : }
663 176 : for (; iCol < nXSize; iCol++)
664 : {
665 68 : Tdest d0 = nCst;
666 204 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
667 : {
668 136 : d0 += LOAD_SRCVAL(iSrc, 0);
669 : }
670 68 : pDest[iCol] = d0;
671 : }
672 : #undef LOAD_SRCVAL
673 : }
674 : }
675 264 : }
676 :
677 : /************************************************************************/
678 : /* OptimizedSumPackedOutput() */
679 : /************************************************************************/
680 :
681 : template <typename Tdest>
682 287 : static bool OptimizedSumPackedOutput(GDALDataType eSrcType, double dfK,
683 : void *pOutBuffer, int nLineSpace,
684 : int nXSize, int nYSize, int nSources,
685 : const void *const *papoSources)
686 : {
687 287 : switch (eSrcType)
688 : {
689 73 : case GDT_Byte:
690 73 : OptimizedSumPackedOutput<uint8_t, Tdest>(dfK, pOutBuffer,
691 : nLineSpace, nXSize, nYSize,
692 : nSources, papoSources);
693 73 : return true;
694 :
695 34 : case GDT_UInt16:
696 34 : OptimizedSumPackedOutput<uint16_t, Tdest>(
697 : dfK, pOutBuffer, nLineSpace, nXSize, nYSize, nSources,
698 : papoSources);
699 34 : return true;
700 :
701 34 : case GDT_Int16:
702 34 : OptimizedSumPackedOutput<int16_t, Tdest>(dfK, pOutBuffer,
703 : nLineSpace, nXSize, nYSize,
704 : nSources, papoSources);
705 34 : return true;
706 :
707 34 : case GDT_Int32:
708 34 : OptimizedSumPackedOutput<int32_t, Tdest>(dfK, pOutBuffer,
709 : nLineSpace, nXSize, nYSize,
710 : nSources, papoSources);
711 34 : return true;
712 :
713 36 : case GDT_Float32:
714 36 : OptimizedSumPackedOutput<float, Tdest>(dfK, pOutBuffer, nLineSpace,
715 : nXSize, nYSize, nSources,
716 : papoSources);
717 36 : return true;
718 :
719 36 : case GDT_Float64:
720 36 : OptimizedSumPackedOutput<double, Tdest>(dfK, pOutBuffer, nLineSpace,
721 : nXSize, nYSize, nSources,
722 : papoSources);
723 36 : return true;
724 :
725 40 : default:
726 40 : break;
727 : }
728 40 : return false;
729 : }
730 :
731 : /************************************************************************/
732 : /* OptimizedSumThroughLargerType() */
733 : /************************************************************************/
734 :
735 : namespace
736 : {
737 : template <typename Tsrc, typename Tdest, typename Enable = void>
738 : struct TintermediateS
739 : {
740 : using type = double;
741 : };
742 :
743 : template <typename Tsrc, typename Tdest>
744 : struct TintermediateS<
745 : Tsrc, Tdest,
746 : std::enable_if_t<
747 : (std::is_same_v<Tsrc, uint8_t> || std::is_same_v<Tsrc, int16_t> ||
748 : std::is_same_v<Tsrc, uint16_t>)&&(std::is_same_v<Tdest, uint8_t> ||
749 : std::is_same_v<Tdest, int16_t> ||
750 : std::is_same_v<Tdest, uint16_t>),
751 : bool>>
752 : {
753 : using type = int32_t;
754 : };
755 :
756 : } // namespace
757 :
758 : template <typename Tsrc, typename Tdest>
759 397 : static bool OptimizedSumThroughLargerType(double dfK, void *pOutBuffer,
760 : int nPixelSpace, int nLineSpace,
761 : int nXSize, int nYSize, int nSources,
762 : const void *const *papoSources)
763 : {
764 : using Tintermediate = typename TintermediateS<Tsrc, Tdest>::type;
765 397 : const Tintermediate k = static_cast<Tintermediate>(dfK);
766 :
767 397 : size_t ii = 0;
768 1189 : for (int iLine = 0; iLine < nYSize; ++iLine)
769 : {
770 792 : GByte *CPL_RESTRICT pDest = static_cast<GByte *>(pOutBuffer) +
771 792 : static_cast<GSpacing>(nLineSpace) * iLine;
772 :
773 792 : constexpr int UNROLLING = 4;
774 792 : int iCol = 0;
775 3952 : for (; iCol < nXSize - (UNROLLING - 1);
776 : iCol += UNROLLING, ii += UNROLLING)
777 : {
778 3160 : Tintermediate aSum[4] = {k, k, k, k};
779 :
780 9480 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
781 : {
782 6320 : aSum[0] += static_cast<const Tsrc *>(papoSources[iSrc])[ii + 0];
783 6320 : aSum[1] += static_cast<const Tsrc *>(papoSources[iSrc])[ii + 1];
784 6320 : aSum[2] += static_cast<const Tsrc *>(papoSources[iSrc])[ii + 2];
785 6320 : aSum[3] += static_cast<const Tsrc *>(papoSources[iSrc])[ii + 3];
786 : }
787 :
788 3160 : GDALCopyWord(aSum[0], *reinterpret_cast<Tdest *>(pDest));
789 3160 : pDest += nPixelSpace;
790 3160 : GDALCopyWord(aSum[1], *reinterpret_cast<Tdest *>(pDest));
791 3160 : pDest += nPixelSpace;
792 3160 : GDALCopyWord(aSum[2], *reinterpret_cast<Tdest *>(pDest));
793 3160 : pDest += nPixelSpace;
794 3160 : GDALCopyWord(aSum[3], *reinterpret_cast<Tdest *>(pDest));
795 3160 : pDest += nPixelSpace;
796 : }
797 1584 : for (; iCol < nXSize; ++iCol, ++ii, pDest += nPixelSpace)
798 : {
799 792 : Tintermediate sum = k;
800 2374 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
801 : {
802 1582 : sum += static_cast<const Tsrc *>(papoSources[iSrc])[ii];
803 : }
804 :
805 792 : auto pDst = reinterpret_cast<Tdest *>(pDest);
806 792 : GDALCopyWord(sum, *pDst);
807 : }
808 : }
809 397 : return true;
810 : }
811 :
812 : /************************************************************************/
813 : /* OptimizedSumThroughLargerType() */
814 : /************************************************************************/
815 :
816 : template <typename Tsrc>
817 500 : static bool OptimizedSumThroughLargerType(GDALDataType eBufType, double dfK,
818 : void *pOutBuffer, int nPixelSpace,
819 : int nLineSpace, int nXSize,
820 : int nYSize, int nSources,
821 : const void *const *papoSources)
822 : {
823 500 : switch (eBufType)
824 : {
825 105 : case GDT_Byte:
826 105 : return OptimizedSumThroughLargerType<Tsrc, uint8_t>(
827 : dfK, pOutBuffer, nPixelSpace, nLineSpace, nXSize, nYSize,
828 105 : nSources, papoSources);
829 :
830 103 : case GDT_UInt16:
831 103 : return OptimizedSumThroughLargerType<Tsrc, uint16_t>(
832 : dfK, pOutBuffer, nPixelSpace, nLineSpace, nXSize, nYSize,
833 103 : nSources, papoSources);
834 :
835 103 : case GDT_Int16:
836 103 : return OptimizedSumThroughLargerType<Tsrc, int16_t>(
837 : dfK, pOutBuffer, nPixelSpace, nLineSpace, nXSize, nYSize,
838 103 : nSources, papoSources);
839 :
840 86 : case GDT_Int32:
841 86 : return OptimizedSumThroughLargerType<Tsrc, int32_t>(
842 : dfK, pOutBuffer, nPixelSpace, nLineSpace, nXSize, nYSize,
843 86 : nSources, papoSources);
844 :
845 : // Float32 and Float64 already covered by OptimizedSum() for packed case
846 103 : default:
847 103 : break;
848 : }
849 103 : return false;
850 : }
851 :
852 : /************************************************************************/
853 : /* OptimizedSumThroughLargerType() */
854 : /************************************************************************/
855 :
856 640 : static bool OptimizedSumThroughLargerType(GDALDataType eSrcType,
857 : GDALDataType eBufType, double dfK,
858 : void *pOutBuffer, int nPixelSpace,
859 : int nLineSpace, int nXSize,
860 : int nYSize, int nSources,
861 : const void *const *papoSources)
862 : {
863 640 : switch (eSrcType)
864 : {
865 70 : case GDT_Byte:
866 70 : return OptimizedSumThroughLargerType<uint8_t>(
867 : eBufType, dfK, pOutBuffer, nPixelSpace, nLineSpace, nXSize,
868 70 : nYSize, nSources, papoSources);
869 :
870 85 : case GDT_UInt16:
871 85 : return OptimizedSumThroughLargerType<uint16_t>(
872 : eBufType, dfK, pOutBuffer, nPixelSpace, nLineSpace, nXSize,
873 85 : nYSize, nSources, papoSources);
874 :
875 85 : case GDT_Int16:
876 85 : return OptimizedSumThroughLargerType<int16_t>(
877 : eBufType, dfK, pOutBuffer, nPixelSpace, nLineSpace, nXSize,
878 85 : nYSize, nSources, papoSources);
879 :
880 85 : case GDT_Int32:
881 85 : return OptimizedSumThroughLargerType<int32_t>(
882 : eBufType, dfK, pOutBuffer, nPixelSpace, nLineSpace, nXSize,
883 85 : nYSize, nSources, papoSources);
884 :
885 90 : case GDT_Float32:
886 90 : return OptimizedSumThroughLargerType<float>(
887 : eBufType, dfK, pOutBuffer, nPixelSpace, nLineSpace, nXSize,
888 90 : nYSize, nSources, papoSources);
889 :
890 85 : case GDT_Float64:
891 85 : return OptimizedSumThroughLargerType<double>(
892 : eBufType, dfK, pOutBuffer, nPixelSpace, nLineSpace, nXSize,
893 85 : nYSize, nSources, papoSources);
894 :
895 140 : default:
896 140 : break;
897 : }
898 :
899 140 : return false;
900 : }
901 :
902 : /************************************************************************/
903 : /* SumPixelFunc() */
904 : /************************************************************************/
905 :
906 : static const char pszSumPixelFuncMetadata[] =
907 : "<PixelFunctionArgumentsList>"
908 : " <Argument name='k' description='Optional constant term' type='double' "
909 : "default='0.0' />"
910 : " <Argument name='propagateNoData' description='Whether the output value "
911 : "should be NoData as as soon as one source is NoData' type='boolean' "
912 : "default='false' />"
913 : " <Argument type='builtin' value='NoData' optional='true' />"
914 : "</PixelFunctionArgumentsList>";
915 :
916 945 : static CPLErr SumPixelFunc(void **papoSources, int nSources, void *pData,
917 : int nXSize, int nYSize, GDALDataType eSrcType,
918 : GDALDataType eBufType, int nPixelSpace,
919 : int nLineSpace, CSLConstList papszArgs)
920 : {
921 : /* ---- Init ---- */
922 945 : if (nSources < 1)
923 : {
924 1 : CPLError(CE_Failure, CPLE_AppDefined,
925 : "sum requires at least one source");
926 1 : return CE_Failure;
927 : }
928 :
929 944 : double dfK = 0.0;
930 944 : if (FetchDoubleArg(papszArgs, "k", &dfK, &dfK) != CE_None)
931 0 : return CE_Failure;
932 :
933 944 : double dfNoData{0};
934 944 : bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
935 944 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
936 0 : return CE_Failure;
937 :
938 944 : const bool bPropagateNoData = CPLTestBool(
939 : CSLFetchNameValueDef(papszArgs, "propagateNoData", "false"));
940 :
941 944 : if (dfNoData == 0 && !bPropagateNoData)
942 940 : bHasNoData = false;
943 :
944 : /* ---- Set pixels ---- */
945 944 : if (GDALDataTypeIsComplex(eSrcType))
946 : {
947 36 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
948 :
949 : /* ---- Set pixels ---- */
950 36 : size_t ii = 0;
951 112 : for (int iLine = 0; iLine < nYSize; ++iLine)
952 : {
953 1296 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
954 : {
955 1220 : double adfSum[2] = {dfK, 0.0};
956 :
957 3690 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
958 : {
959 2470 : const void *const pReal = papoSources[iSrc];
960 2470 : const void *const pImag =
961 2470 : static_cast<const GByte *>(pReal) + nOffset;
962 :
963 : // Source raster pixels may be obtained with GetSrcVal
964 : // macro.
965 2470 : adfSum[0] += GetSrcVal(pReal, eSrcType, ii);
966 2470 : adfSum[1] += GetSrcVal(pImag, eSrcType, ii);
967 : }
968 :
969 1220 : GDALCopyWords(adfSum, GDT_CFloat64, 0,
970 : static_cast<GByte *>(pData) +
971 1220 : static_cast<GSpacing>(nLineSpace) * iLine +
972 1220 : iCol * nPixelSpace,
973 : eBufType, nPixelSpace, 1);
974 : }
975 : }
976 : }
977 : else
978 : {
979 : /* ---- Set pixels ---- */
980 908 : bool bGeneralCase = true;
981 908 : if (dfNoData == 0 && !bPropagateNoData)
982 : {
983 904 : if (eBufType == GDT_Float32 && nPixelSpace == sizeof(float))
984 : {
985 126 : bGeneralCase = !OptimizedSumPackedOutput<float>(
986 : eSrcType, dfK, pData, nLineSpace, nXSize, nYSize, nSources,
987 : papoSources);
988 : }
989 778 : else if (eBufType == GDT_Float64 && nPixelSpace == sizeof(double))
990 : {
991 161 : bGeneralCase = !OptimizedSumPackedOutput<double>(
992 : eSrcType, dfK, pData, nLineSpace, nXSize, nYSize, nSources,
993 : papoSources);
994 : }
995 617 : else if (
996 617 : dfK >= 0 && dfK <= INT_MAX && eBufType == GDT_Int32 &&
997 1251 : nPixelSpace == sizeof(int32_t) && eSrcType == GDT_Byte &&
998 : // Limitation to avoid overflow of int32 if all source values are at the max of their data type
999 17 : nSources <=
1000 17 : (INT_MAX - dfK) / std::numeric_limits<uint8_t>::max())
1001 : {
1002 17 : bGeneralCase = false;
1003 17 : OptimizedSumPackedOutput<uint8_t, int32_t>(
1004 : dfK, pData, nLineSpace, nXSize, nYSize, nSources,
1005 : papoSources);
1006 : }
1007 :
1008 1544 : if (bGeneralCase && dfK >= 0 && dfK <= INT_MAX &&
1009 640 : nSources <=
1010 640 : (INT_MAX - dfK) / std::numeric_limits<uint16_t>::max())
1011 : {
1012 640 : bGeneralCase = !OptimizedSumThroughLargerType(
1013 : eSrcType, eBufType, dfK, pData, nPixelSpace, nLineSpace,
1014 : nXSize, nYSize, nSources, papoSources);
1015 : }
1016 : }
1017 :
1018 908 : if (bGeneralCase)
1019 : {
1020 247 : size_t ii = 0;
1021 737 : for (int iLine = 0; iLine < nYSize; ++iLine)
1022 : {
1023 8765 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1024 : {
1025 8275 : double dfSum = dfK;
1026 24824 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
1027 : {
1028 : const double dfVal =
1029 16551 : GetSrcVal(papoSources[iSrc], eSrcType, ii);
1030 :
1031 16551 : if (bHasNoData && IsNoData(dfVal, dfNoData))
1032 : {
1033 12 : if (bPropagateNoData)
1034 : {
1035 2 : dfSum = dfNoData;
1036 2 : break;
1037 : }
1038 : }
1039 : else
1040 : {
1041 16539 : dfSum += dfVal;
1042 : }
1043 : }
1044 :
1045 8275 : GDALCopyWords(&dfSum, GDT_Float64, 0,
1046 : static_cast<GByte *>(pData) +
1047 8275 : static_cast<GSpacing>(nLineSpace) *
1048 8275 : iLine +
1049 8275 : iCol * nPixelSpace,
1050 : eBufType, nPixelSpace, 1);
1051 : }
1052 : }
1053 : }
1054 : }
1055 :
1056 : /* ---- Return success ---- */
1057 944 : return CE_None;
1058 : } /* SumPixelFunc */
1059 :
1060 : static const char pszDiffPixelFuncMetadata[] =
1061 : "<PixelFunctionArgumentsList>"
1062 : " <Argument type='builtin' value='NoData' optional='true' />"
1063 : "</PixelFunctionArgumentsList>";
1064 :
1065 5 : static CPLErr DiffPixelFunc(void **papoSources, int nSources, void *pData,
1066 : int nXSize, int nYSize, GDALDataType eSrcType,
1067 : GDALDataType eBufType, int nPixelSpace,
1068 : int nLineSpace, CSLConstList papszArgs)
1069 : {
1070 : /* ---- Init ---- */
1071 5 : if (nSources != 2)
1072 1 : return CE_Failure;
1073 :
1074 4 : double dfNoData{0};
1075 4 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1076 4 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1077 0 : return CE_Failure;
1078 :
1079 4 : if (GDALDataTypeIsComplex(eSrcType))
1080 : {
1081 1 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
1082 1 : const void *const pReal0 = papoSources[0];
1083 1 : const void *const pImag0 =
1084 1 : static_cast<GByte *>(papoSources[0]) + nOffset;
1085 1 : const void *const pReal1 = papoSources[1];
1086 1 : const void *const pImag1 =
1087 1 : static_cast<GByte *>(papoSources[1]) + nOffset;
1088 :
1089 : /* ---- Set pixels ---- */
1090 1 : size_t ii = 0;
1091 7 : for (int iLine = 0; iLine < nYSize; ++iLine)
1092 : {
1093 36 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1094 : {
1095 : // Source raster pixels may be obtained with GetSrcVal macro.
1096 30 : double adfPixVal[2] = {GetSrcVal(pReal0, eSrcType, ii) -
1097 30 : GetSrcVal(pReal1, eSrcType, ii),
1098 90 : GetSrcVal(pImag0, eSrcType, ii) -
1099 30 : GetSrcVal(pImag1, eSrcType, ii)};
1100 :
1101 30 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
1102 : static_cast<GByte *>(pData) +
1103 30 : static_cast<GSpacing>(nLineSpace) * iLine +
1104 30 : iCol * nPixelSpace,
1105 : eBufType, nPixelSpace, 1);
1106 : }
1107 : }
1108 : }
1109 : else
1110 : {
1111 : /* ---- Set pixels ---- */
1112 3 : size_t ii = 0;
1113 11 : for (int iLine = 0; iLine < nYSize; ++iLine)
1114 : {
1115 42 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1116 : {
1117 34 : const double dfA = GetSrcVal(papoSources[0], eSrcType, ii);
1118 34 : const double dfB = GetSrcVal(papoSources[1], eSrcType, ii);
1119 :
1120 : const double dfPixVal =
1121 38 : bHasNoData &&
1122 4 : (IsNoData(dfA, dfNoData) || IsNoData(dfB, dfNoData))
1123 38 : ? dfNoData
1124 34 : : dfA - dfB;
1125 :
1126 34 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1127 : static_cast<GByte *>(pData) +
1128 34 : static_cast<GSpacing>(nLineSpace) * iLine +
1129 34 : iCol * nPixelSpace,
1130 : eBufType, nPixelSpace, 1);
1131 : }
1132 : }
1133 : }
1134 :
1135 : /* ---- Return success ---- */
1136 4 : return CE_None;
1137 : } // DiffPixelFunc
1138 :
1139 : static const char pszMulPixelFuncMetadata[] =
1140 : "<PixelFunctionArgumentsList>"
1141 : " <Argument name='k' description='Optional constant factor' "
1142 : "type='double' default='1.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 11 : static CPLErr MulPixelFunc(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 11 : if (nSources < 2 && CSLFetchNameValue(papszArgs, "k") == nullptr)
1156 : {
1157 1 : CPLError(CE_Failure, CPLE_AppDefined,
1158 : "mul requires at least two sources or a specified constant k");
1159 1 : return CE_Failure;
1160 : }
1161 :
1162 10 : double dfK = 1.0;
1163 10 : if (FetchDoubleArg(papszArgs, "k", &dfK, &dfK) != CE_None)
1164 0 : return CE_Failure;
1165 :
1166 10 : double dfNoData{0};
1167 10 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1168 10 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1169 0 : return CE_Failure;
1170 :
1171 10 : const bool bPropagateNoData = CPLTestBool(
1172 : CSLFetchNameValueDef(papszArgs, "propagateNoData", "false"));
1173 :
1174 : /* ---- Set pixels ---- */
1175 10 : if (GDALDataTypeIsComplex(eSrcType))
1176 : {
1177 1 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
1178 :
1179 : /* ---- Set pixels ---- */
1180 1 : size_t ii = 0;
1181 7 : for (int iLine = 0; iLine < nYSize; ++iLine)
1182 : {
1183 36 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1184 : {
1185 30 : double adfPixVal[2] = {dfK, 0.0};
1186 :
1187 90 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
1188 : {
1189 60 : const void *const pReal = papoSources[iSrc];
1190 60 : const void *const pImag =
1191 60 : static_cast<const GByte *>(pReal) + nOffset;
1192 :
1193 60 : const double dfOldR = adfPixVal[0];
1194 60 : const double dfOldI = adfPixVal[1];
1195 :
1196 : // Source raster pixels may be obtained with GetSrcVal
1197 : // macro.
1198 60 : const double dfNewR = GetSrcVal(pReal, eSrcType, ii);
1199 60 : const double dfNewI = GetSrcVal(pImag, eSrcType, ii);
1200 :
1201 60 : adfPixVal[0] = dfOldR * dfNewR - dfOldI * dfNewI;
1202 60 : adfPixVal[1] = dfOldR * dfNewI + dfOldI * dfNewR;
1203 : }
1204 :
1205 30 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
1206 : static_cast<GByte *>(pData) +
1207 30 : static_cast<GSpacing>(nLineSpace) * iLine +
1208 30 : iCol * nPixelSpace,
1209 : eBufType, nPixelSpace, 1);
1210 : }
1211 : }
1212 : }
1213 : else
1214 : {
1215 : /* ---- Set pixels ---- */
1216 9 : size_t ii = 0;
1217 75 : for (int iLine = 0; iLine < nYSize; ++iLine)
1218 : {
1219 1278 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1220 : {
1221 1212 : double dfPixVal = dfK; // Not complex.
1222 :
1223 4033 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
1224 : {
1225 : const double dfVal =
1226 2827 : GetSrcVal(papoSources[iSrc], eSrcType, ii);
1227 :
1228 2827 : if (bHasNoData && IsNoData(dfVal, dfNoData))
1229 : {
1230 18 : if (bPropagateNoData)
1231 : {
1232 6 : dfPixVal = dfNoData;
1233 6 : break;
1234 : }
1235 : }
1236 : else
1237 : {
1238 2809 : dfPixVal *= dfVal;
1239 : }
1240 : }
1241 :
1242 1212 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1243 : static_cast<GByte *>(pData) +
1244 1212 : static_cast<GSpacing>(nLineSpace) * iLine +
1245 1212 : iCol * nPixelSpace,
1246 : eBufType, nPixelSpace, 1);
1247 : }
1248 : }
1249 : }
1250 :
1251 : /* ---- Return success ---- */
1252 10 : return CE_None;
1253 : } // MulPixelFunc
1254 :
1255 : static const char pszDivPixelFuncMetadata[] =
1256 : "<PixelFunctionArgumentsList>"
1257 : " "
1258 : "<Argument type='builtin' value='NoData' optional='true' />"
1259 : "</PixelFunctionArgumentsList>";
1260 :
1261 7 : static CPLErr DivPixelFunc(void **papoSources, int nSources, void *pData,
1262 : int nXSize, int nYSize, GDALDataType eSrcType,
1263 : GDALDataType eBufType, int nPixelSpace,
1264 : int nLineSpace, CSLConstList papszArgs)
1265 : {
1266 : /* ---- Init ---- */
1267 7 : if (nSources != 2)
1268 0 : return CE_Failure;
1269 :
1270 7 : double dfNoData{0};
1271 7 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1272 7 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1273 0 : return CE_Failure;
1274 :
1275 : /* ---- Set pixels ---- */
1276 7 : if (GDALDataTypeIsComplex(eSrcType))
1277 : {
1278 1 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
1279 1 : const void *const pReal0 = papoSources[0];
1280 1 : const void *const pImag0 =
1281 1 : static_cast<GByte *>(papoSources[0]) + nOffset;
1282 1 : const void *const pReal1 = papoSources[1];
1283 1 : const void *const pImag1 =
1284 1 : static_cast<GByte *>(papoSources[1]) + nOffset;
1285 :
1286 : /* ---- Set pixels ---- */
1287 1 : size_t ii = 0;
1288 7 : for (int iLine = 0; iLine < nYSize; ++iLine)
1289 : {
1290 36 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1291 : {
1292 : // Source raster pixels may be obtained with GetSrcVal macro.
1293 30 : const double dfReal0 = GetSrcVal(pReal0, eSrcType, ii);
1294 30 : const double dfReal1 = GetSrcVal(pReal1, eSrcType, ii);
1295 30 : const double dfImag0 = GetSrcVal(pImag0, eSrcType, ii);
1296 30 : const double dfImag1 = GetSrcVal(pImag1, eSrcType, ii);
1297 30 : const double dfAux = dfReal1 * dfReal1 + dfImag1 * dfImag1;
1298 :
1299 : const double adfPixVal[2] = {
1300 : dfAux == 0
1301 30 : ? std::numeric_limits<double>::infinity()
1302 30 : : dfReal0 * dfReal1 / dfAux + dfImag0 * dfImag1 / dfAux,
1303 0 : dfAux == 0 ? std::numeric_limits<double>::infinity()
1304 30 : : dfReal1 / dfAux * dfImag0 -
1305 30 : dfReal0 * dfImag1 / dfAux};
1306 :
1307 30 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
1308 : static_cast<GByte *>(pData) +
1309 30 : static_cast<GSpacing>(nLineSpace) * iLine +
1310 30 : iCol * nPixelSpace,
1311 : eBufType, nPixelSpace, 1);
1312 : }
1313 : }
1314 : }
1315 : else
1316 : {
1317 : /* ---- Set pixels ---- */
1318 6 : size_t ii = 0;
1319 17 : for (int iLine = 0; iLine < nYSize; ++iLine)
1320 : {
1321 48 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1322 : {
1323 37 : const double dfNum = GetSrcVal(papoSources[0], eSrcType, ii);
1324 37 : const double dfDenom = GetSrcVal(papoSources[1], eSrcType, ii);
1325 :
1326 37 : double dfPixVal = dfNoData;
1327 41 : if (!bHasNoData || (!IsNoData(dfNum, dfNoData) &&
1328 4 : !IsNoData(dfDenom, dfNoData)))
1329 : {
1330 : // coverity[divide_by_zero]
1331 33 : dfPixVal = dfDenom == 0
1332 33 : ? std::numeric_limits<double>::infinity()
1333 : : dfNum / dfDenom;
1334 : }
1335 :
1336 37 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1337 : static_cast<GByte *>(pData) +
1338 37 : static_cast<GSpacing>(nLineSpace) * iLine +
1339 37 : iCol * nPixelSpace,
1340 : eBufType, nPixelSpace, 1);
1341 : }
1342 : }
1343 : }
1344 :
1345 : /* ---- Return success ---- */
1346 7 : return CE_None;
1347 : } // DivPixelFunc
1348 :
1349 3 : static CPLErr CMulPixelFunc(void **papoSources, int nSources, void *pData,
1350 : int nXSize, int nYSize, GDALDataType eSrcType,
1351 : GDALDataType eBufType, int nPixelSpace,
1352 : int nLineSpace)
1353 : {
1354 : /* ---- Init ---- */
1355 3 : if (nSources != 2)
1356 1 : return CE_Failure;
1357 :
1358 : /* ---- Set pixels ---- */
1359 2 : if (GDALDataTypeIsComplex(eSrcType))
1360 : {
1361 1 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
1362 1 : const void *const pReal0 = papoSources[0];
1363 1 : const void *const pImag0 =
1364 1 : static_cast<GByte *>(papoSources[0]) + nOffset;
1365 1 : const void *const pReal1 = papoSources[1];
1366 1 : const void *const pImag1 =
1367 1 : static_cast<GByte *>(papoSources[1]) + nOffset;
1368 :
1369 1 : size_t ii = 0;
1370 7 : for (int iLine = 0; iLine < nYSize; ++iLine)
1371 : {
1372 36 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1373 : {
1374 : // Source raster pixels may be obtained with GetSrcVal macro.
1375 30 : const double dfReal0 = GetSrcVal(pReal0, eSrcType, ii);
1376 30 : const double dfReal1 = GetSrcVal(pReal1, eSrcType, ii);
1377 30 : const double dfImag0 = GetSrcVal(pImag0, eSrcType, ii);
1378 30 : const double dfImag1 = GetSrcVal(pImag1, eSrcType, ii);
1379 : const double adfPixVal[2] = {
1380 30 : dfReal0 * dfReal1 + dfImag0 * dfImag1,
1381 30 : dfReal1 * dfImag0 - dfReal0 * dfImag1};
1382 :
1383 30 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
1384 : static_cast<GByte *>(pData) +
1385 30 : static_cast<GSpacing>(nLineSpace) * iLine +
1386 30 : iCol * nPixelSpace,
1387 : eBufType, nPixelSpace, 1);
1388 : }
1389 : }
1390 : }
1391 : else
1392 : {
1393 1 : size_t ii = 0;
1394 21 : for (int iLine = 0; iLine < nYSize; ++iLine)
1395 : {
1396 420 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1397 : {
1398 : // Source raster pixels may be obtained with GetSrcVal macro.
1399 : // Not complex.
1400 400 : const double adfPixVal[2] = {
1401 400 : GetSrcVal(papoSources[0], eSrcType, ii) *
1402 400 : GetSrcVal(papoSources[1], eSrcType, ii),
1403 400 : 0.0};
1404 :
1405 400 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
1406 : static_cast<GByte *>(pData) +
1407 400 : static_cast<GSpacing>(nLineSpace) * iLine +
1408 400 : iCol * nPixelSpace,
1409 : eBufType, nPixelSpace, 1);
1410 : }
1411 : }
1412 : }
1413 :
1414 : /* ---- Return success ---- */
1415 2 : return CE_None;
1416 : } // CMulPixelFunc
1417 :
1418 : static const char pszInvPixelFuncMetadata[] =
1419 : "<PixelFunctionArgumentsList>"
1420 : " <Argument name='k' description='Optional constant factor' "
1421 : "type='double' default='1.0' />"
1422 : " "
1423 : "<Argument type='builtin' value='NoData' optional='true' />"
1424 : "</PixelFunctionArgumentsList>";
1425 :
1426 9 : static CPLErr InvPixelFunc(void **papoSources, int nSources, void *pData,
1427 : int nXSize, int nYSize, GDALDataType eSrcType,
1428 : GDALDataType eBufType, int nPixelSpace,
1429 : int nLineSpace, CSLConstList papszArgs)
1430 : {
1431 : /* ---- Init ---- */
1432 9 : if (nSources != 1)
1433 1 : return CE_Failure;
1434 :
1435 8 : double dfK = 1.0;
1436 8 : if (FetchDoubleArg(papszArgs, "k", &dfK, &dfK) != CE_None)
1437 0 : return CE_Failure;
1438 :
1439 8 : double dfNoData{0};
1440 8 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1441 8 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1442 0 : return CE_Failure;
1443 :
1444 : /* ---- Set pixels ---- */
1445 8 : if (GDALDataTypeIsComplex(eSrcType))
1446 : {
1447 2 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
1448 2 : const void *const pReal = papoSources[0];
1449 2 : const void *const pImag =
1450 2 : static_cast<GByte *>(papoSources[0]) + nOffset;
1451 :
1452 2 : size_t ii = 0;
1453 9 : for (int iLine = 0; iLine < nYSize; ++iLine)
1454 : {
1455 38 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1456 : {
1457 : // Source raster pixels may be obtained with GetSrcVal macro.
1458 31 : const double dfReal = GetSrcVal(pReal, eSrcType, ii);
1459 31 : const double dfImag = GetSrcVal(pImag, eSrcType, ii);
1460 31 : const double dfAux = dfReal * dfReal + dfImag * dfImag;
1461 : const double adfPixVal[2] = {
1462 31 : dfAux == 0 ? std::numeric_limits<double>::infinity()
1463 30 : : dfK * dfReal / dfAux,
1464 1 : dfAux == 0 ? std::numeric_limits<double>::infinity()
1465 31 : : -dfK * dfImag / dfAux};
1466 :
1467 31 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
1468 : static_cast<GByte *>(pData) +
1469 31 : static_cast<GSpacing>(nLineSpace) * iLine +
1470 31 : iCol * nPixelSpace,
1471 : eBufType, nPixelSpace, 1);
1472 : }
1473 : }
1474 : }
1475 : else
1476 : {
1477 : /* ---- Set pixels ---- */
1478 6 : size_t ii = 0;
1479 50 : for (int iLine = 0; iLine < nYSize; ++iLine)
1480 : {
1481 851 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1482 : {
1483 : // Source raster pixels may be obtained with GetSrcVal macro.
1484 : // Not complex.
1485 807 : const double dfVal = GetSrcVal(papoSources[0], eSrcType, ii);
1486 807 : double dfPixVal = dfNoData;
1487 :
1488 807 : if (!bHasNoData || !IsNoData(dfVal, dfNoData))
1489 : {
1490 802 : dfPixVal = dfVal == 0
1491 802 : ? std::numeric_limits<double>::infinity()
1492 801 : : dfK / dfVal;
1493 : }
1494 :
1495 807 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1496 : static_cast<GByte *>(pData) +
1497 807 : static_cast<GSpacing>(nLineSpace) * iLine +
1498 807 : iCol * nPixelSpace,
1499 : eBufType, nPixelSpace, 1);
1500 : }
1501 : }
1502 : }
1503 :
1504 : /* ---- Return success ---- */
1505 8 : return CE_None;
1506 : } // InvPixelFunc
1507 :
1508 4 : static CPLErr IntensityPixelFunc(void **papoSources, int nSources, void *pData,
1509 : int nXSize, int nYSize, GDALDataType eSrcType,
1510 : GDALDataType eBufType, int nPixelSpace,
1511 : int nLineSpace)
1512 : {
1513 : /* ---- Init ---- */
1514 4 : if (nSources != 1)
1515 1 : return CE_Failure;
1516 :
1517 3 : if (GDALDataTypeIsComplex(eSrcType))
1518 : {
1519 2 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
1520 2 : const void *const pReal = papoSources[0];
1521 2 : const void *const pImag =
1522 2 : static_cast<GByte *>(papoSources[0]) + nOffset;
1523 :
1524 : /* ---- Set pixels ---- */
1525 2 : size_t ii = 0;
1526 14 : for (int iLine = 0; iLine < nYSize; ++iLine)
1527 : {
1528 72 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1529 : {
1530 : // Source raster pixels may be obtained with GetSrcVal macro.
1531 60 : const double dfReal = GetSrcVal(pReal, eSrcType, ii);
1532 60 : const double dfImag = GetSrcVal(pImag, eSrcType, ii);
1533 :
1534 60 : const double dfPixVal = dfReal * dfReal + dfImag * dfImag;
1535 :
1536 60 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1537 : static_cast<GByte *>(pData) +
1538 60 : static_cast<GSpacing>(nLineSpace) * iLine +
1539 60 : iCol * nPixelSpace,
1540 : eBufType, nPixelSpace, 1);
1541 : }
1542 : }
1543 : }
1544 : else
1545 : {
1546 : /* ---- Set pixels ---- */
1547 1 : size_t ii = 0;
1548 21 : for (int iLine = 0; iLine < nYSize; ++iLine)
1549 : {
1550 420 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1551 : {
1552 : // Source raster pixels may be obtained with GetSrcVal macro.
1553 400 : double dfPixVal = GetSrcVal(papoSources[0], eSrcType, ii);
1554 400 : dfPixVal *= dfPixVal;
1555 :
1556 400 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1557 : static_cast<GByte *>(pData) +
1558 400 : static_cast<GSpacing>(nLineSpace) * iLine +
1559 400 : iCol * nPixelSpace,
1560 : eBufType, nPixelSpace, 1);
1561 : }
1562 : }
1563 : }
1564 :
1565 : /* ---- Return success ---- */
1566 3 : return CE_None;
1567 : } // IntensityPixelFunc
1568 :
1569 : static const char pszSqrtPixelFuncMetadata[] =
1570 : "<PixelFunctionArgumentsList>"
1571 : " <Argument type='builtin' value='NoData' optional='true'/>"
1572 : "</PixelFunctionArgumentsList>";
1573 :
1574 3 : static CPLErr SqrtPixelFunc(void **papoSources, int nSources, void *pData,
1575 : int nXSize, int nYSize, GDALDataType eSrcType,
1576 : GDALDataType eBufType, int nPixelSpace,
1577 : int nLineSpace, CSLConstList papszArgs)
1578 : {
1579 : /* ---- Init ---- */
1580 3 : if (nSources != 1)
1581 1 : return CE_Failure;
1582 2 : if (GDALDataTypeIsComplex(eSrcType))
1583 0 : return CE_Failure;
1584 :
1585 2 : double dfNoData{0};
1586 2 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1587 2 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1588 0 : return CE_Failure;
1589 :
1590 : /* ---- Set pixels ---- */
1591 2 : size_t ii = 0;
1592 23 : for (int iLine = 0; iLine < nYSize; ++iLine)
1593 : {
1594 423 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1595 : {
1596 : // Source raster pixels may be obtained with GetSrcVal macro.
1597 402 : double dfPixVal = GetSrcVal(papoSources[0], eSrcType, ii);
1598 :
1599 402 : if (bHasNoData && IsNoData(dfPixVal, dfNoData))
1600 : {
1601 2 : dfPixVal = dfNoData;
1602 : }
1603 : else
1604 : {
1605 400 : dfPixVal = std::sqrt(dfPixVal);
1606 : }
1607 :
1608 402 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1609 : static_cast<GByte *>(pData) +
1610 402 : static_cast<GSpacing>(nLineSpace) * iLine +
1611 402 : iCol * nPixelSpace,
1612 : eBufType, nPixelSpace, 1);
1613 : }
1614 : }
1615 :
1616 : /* ---- Return success ---- */
1617 2 : return CE_None;
1618 : } // SqrtPixelFunc
1619 :
1620 13 : static CPLErr Log10PixelFuncHelper(void **papoSources, int nSources,
1621 : void *pData, int nXSize, int nYSize,
1622 : GDALDataType eSrcType, GDALDataType eBufType,
1623 : int nPixelSpace, int nLineSpace,
1624 : CSLConstList papszArgs, double fact)
1625 : {
1626 : /* ---- Init ---- */
1627 13 : if (nSources != 1)
1628 2 : return CE_Failure;
1629 :
1630 11 : double dfNoData{0};
1631 11 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1632 11 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1633 0 : return CE_Failure;
1634 :
1635 11 : if (GDALDataTypeIsComplex(eSrcType))
1636 : {
1637 : // Complex input datatype.
1638 5 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
1639 5 : const void *const pReal = papoSources[0];
1640 5 : const void *const pImag =
1641 5 : static_cast<GByte *>(papoSources[0]) + nOffset;
1642 :
1643 : /* We should compute fact * log10( sqrt( dfReal * dfReal + dfImag *
1644 : * dfImag ) ) */
1645 : /* Given that log10(sqrt(x)) = 0.5 * log10(x) */
1646 : /* we can remove the sqrt() by multiplying fact by 0.5 */
1647 5 : fact *= 0.5;
1648 :
1649 : /* ---- Set pixels ---- */
1650 5 : size_t ii = 0;
1651 35 : for (int iLine = 0; iLine < nYSize; ++iLine)
1652 : {
1653 180 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1654 : {
1655 : // Source raster pixels may be obtained with GetSrcVal macro.
1656 150 : const double dfReal = GetSrcVal(pReal, eSrcType, ii);
1657 150 : const double dfImag = GetSrcVal(pImag, eSrcType, ii);
1658 :
1659 : const double dfPixVal =
1660 150 : fact * log10(dfReal * dfReal + dfImag * dfImag);
1661 :
1662 150 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1663 : static_cast<GByte *>(pData) +
1664 150 : static_cast<GSpacing>(nLineSpace) * iLine +
1665 150 : iCol * nPixelSpace,
1666 : eBufType, nPixelSpace, 1);
1667 : }
1668 : }
1669 : }
1670 : else
1671 : {
1672 : /* ---- Set pixels ---- */
1673 6 : size_t ii = 0;
1674 88 : for (int iLine = 0; iLine < nYSize; ++iLine)
1675 : {
1676 1686 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1677 : {
1678 : // Source raster pixels may be obtained with GetSrcVal macro.
1679 1604 : const double dfSrcVal = GetSrcVal(papoSources[0], eSrcType, ii);
1680 : const double dfPixVal =
1681 4 : bHasNoData && IsNoData(dfSrcVal, dfNoData)
1682 1608 : ? dfNoData
1683 1600 : : fact * std::log10(std::abs(dfSrcVal));
1684 :
1685 1604 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1686 : static_cast<GByte *>(pData) +
1687 1604 : static_cast<GSpacing>(nLineSpace) * iLine +
1688 1604 : iCol * nPixelSpace,
1689 : eBufType, nPixelSpace, 1);
1690 : }
1691 : }
1692 : }
1693 :
1694 : /* ---- Return success ---- */
1695 11 : return CE_None;
1696 : } // Log10PixelFuncHelper
1697 :
1698 : static const char pszLog10PixelFuncMetadata[] =
1699 : "<PixelFunctionArgumentsList>"
1700 : " <Argument type='builtin' value='NoData' optional='true'/>"
1701 : "</PixelFunctionArgumentsList>";
1702 :
1703 5 : static CPLErr Log10PixelFunc(void **papoSources, int nSources, void *pData,
1704 : int nXSize, int nYSize, GDALDataType eSrcType,
1705 : GDALDataType eBufType, int nPixelSpace,
1706 : int nLineSpace, CSLConstList papszArgs)
1707 : {
1708 5 : return Log10PixelFuncHelper(papoSources, nSources, pData, nXSize, nYSize,
1709 : eSrcType, eBufType, nPixelSpace, nLineSpace,
1710 5 : papszArgs, 1.0);
1711 : } // Log10PixelFunc
1712 :
1713 : static const char pszDBPixelFuncMetadata[] =
1714 : "<PixelFunctionArgumentsList>"
1715 : " <Argument name='fact' description='Factor' type='double' "
1716 : "default='20.0' />"
1717 : " <Argument type='builtin' value='NoData' optional='true' />"
1718 : "</PixelFunctionArgumentsList>";
1719 :
1720 8 : static CPLErr DBPixelFunc(void **papoSources, int nSources, void *pData,
1721 : int nXSize, int nYSize, GDALDataType eSrcType,
1722 : GDALDataType eBufType, int nPixelSpace,
1723 : int nLineSpace, CSLConstList papszArgs)
1724 : {
1725 8 : double dfFact = 20.;
1726 8 : if (FetchDoubleArg(papszArgs, "fact", &dfFact, &dfFact) != CE_None)
1727 0 : return CE_Failure;
1728 :
1729 8 : return Log10PixelFuncHelper(papoSources, nSources, pData, nXSize, nYSize,
1730 : eSrcType, eBufType, nPixelSpace, nLineSpace,
1731 8 : papszArgs, dfFact);
1732 : } // DBPixelFunc
1733 :
1734 8 : static CPLErr ExpPixelFuncHelper(void **papoSources, int nSources, void *pData,
1735 : int nXSize, int nYSize, GDALDataType eSrcType,
1736 : GDALDataType eBufType, int nPixelSpace,
1737 : int nLineSpace, CSLConstList papszArgs,
1738 : double base, double fact)
1739 : {
1740 : /* ---- Init ---- */
1741 8 : if (nSources != 1)
1742 2 : return CE_Failure;
1743 6 : if (GDALDataTypeIsComplex(eSrcType))
1744 0 : return CE_Failure;
1745 :
1746 6 : double dfNoData{0};
1747 6 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1748 6 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1749 0 : return CE_Failure;
1750 :
1751 : /* ---- Set pixels ---- */
1752 6 : size_t ii = 0;
1753 107 : for (int iLine = 0; iLine < nYSize; ++iLine)
1754 : {
1755 2103 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1756 : {
1757 : // Source raster pixels may be obtained with GetSrcVal macro.
1758 2002 : const double dfVal = GetSrcVal(papoSources[0], eSrcType, ii);
1759 2 : const double dfPixVal = bHasNoData && IsNoData(dfVal, dfNoData)
1760 2004 : ? dfNoData
1761 2000 : : pow(base, dfVal * fact);
1762 :
1763 2002 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1764 : static_cast<GByte *>(pData) +
1765 2002 : static_cast<GSpacing>(nLineSpace) * iLine +
1766 2002 : iCol * nPixelSpace,
1767 : eBufType, nPixelSpace, 1);
1768 : }
1769 : }
1770 :
1771 : /* ---- Return success ---- */
1772 6 : return CE_None;
1773 : } // ExpPixelFuncHelper
1774 :
1775 : static const char pszExpPixelFuncMetadata[] =
1776 : "<PixelFunctionArgumentsList>"
1777 : " <Argument name='base' description='Base' type='double' "
1778 : "default='2.7182818284590452353602874713526624' />"
1779 : " <Argument name='fact' description='Factor' type='double' default='1' />"
1780 : " <Argument type='builtin' value='NoData' optional='true' />"
1781 : "</PixelFunctionArgumentsList>";
1782 :
1783 4 : static CPLErr ExpPixelFunc(void **papoSources, int nSources, void *pData,
1784 : int nXSize, int nYSize, GDALDataType eSrcType,
1785 : GDALDataType eBufType, int nPixelSpace,
1786 : int nLineSpace, CSLConstList papszArgs)
1787 : {
1788 4 : double dfBase = 2.7182818284590452353602874713526624;
1789 4 : double dfFact = 1.;
1790 :
1791 4 : if (FetchDoubleArg(papszArgs, "base", &dfBase, &dfBase) != CE_None)
1792 0 : return CE_Failure;
1793 :
1794 4 : if (FetchDoubleArg(papszArgs, "fact", &dfFact, &dfFact) != CE_None)
1795 0 : return CE_Failure;
1796 :
1797 4 : return ExpPixelFuncHelper(papoSources, nSources, pData, nXSize, nYSize,
1798 : eSrcType, eBufType, nPixelSpace, nLineSpace,
1799 4 : papszArgs, dfBase, dfFact);
1800 : } // ExpPixelFunc
1801 :
1802 2 : static CPLErr dB2AmpPixelFunc(void **papoSources, int nSources, void *pData,
1803 : int nXSize, int nYSize, GDALDataType eSrcType,
1804 : GDALDataType eBufType, int nPixelSpace,
1805 : int nLineSpace)
1806 : {
1807 2 : return ExpPixelFuncHelper(papoSources, nSources, pData, nXSize, nYSize,
1808 : eSrcType, eBufType, nPixelSpace, nLineSpace,
1809 2 : nullptr, 10.0, 1. / 20);
1810 : } // dB2AmpPixelFunc
1811 :
1812 2 : static CPLErr dB2PowPixelFunc(void **papoSources, int nSources, void *pData,
1813 : int nXSize, int nYSize, GDALDataType eSrcType,
1814 : GDALDataType eBufType, int nPixelSpace,
1815 : int nLineSpace)
1816 : {
1817 2 : return ExpPixelFuncHelper(papoSources, nSources, pData, nXSize, nYSize,
1818 : eSrcType, eBufType, nPixelSpace, nLineSpace,
1819 2 : nullptr, 10.0, 1. / 10);
1820 : } // dB2PowPixelFunc
1821 :
1822 : static const char pszPowPixelFuncMetadata[] =
1823 : "<PixelFunctionArgumentsList>"
1824 : " <Argument name='power' description='Exponent' type='double' "
1825 : "mandatory='1' />"
1826 : " <Argument type='builtin' value='NoData' optional='true' />"
1827 : "</PixelFunctionArgumentsList>";
1828 :
1829 2 : static CPLErr PowPixelFunc(void **papoSources, int nSources, void *pData,
1830 : int nXSize, int nYSize, GDALDataType eSrcType,
1831 : GDALDataType eBufType, int nPixelSpace,
1832 : int nLineSpace, CSLConstList papszArgs)
1833 : {
1834 : /* ---- Init ---- */
1835 2 : if (nSources != 1)
1836 0 : return CE_Failure;
1837 2 : if (GDALDataTypeIsComplex(eSrcType))
1838 0 : return CE_Failure;
1839 :
1840 : double power;
1841 2 : if (FetchDoubleArg(papszArgs, "power", &power) != CE_None)
1842 0 : return CE_Failure;
1843 :
1844 2 : double dfNoData{0};
1845 2 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1846 2 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1847 0 : return CE_Failure;
1848 :
1849 : /* ---- Set pixels ---- */
1850 2 : size_t ii = 0;
1851 23 : for (int iLine = 0; iLine < nYSize; ++iLine)
1852 : {
1853 423 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1854 : {
1855 402 : const double dfVal = GetSrcVal(papoSources[0], eSrcType, ii);
1856 :
1857 2 : const double dfPixVal = bHasNoData && IsNoData(dfVal, dfNoData)
1858 404 : ? dfNoData
1859 400 : : std::pow(dfVal, power);
1860 :
1861 402 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1862 : static_cast<GByte *>(pData) +
1863 402 : static_cast<GSpacing>(nLineSpace) * iLine +
1864 402 : iCol * nPixelSpace,
1865 : eBufType, nPixelSpace, 1);
1866 : }
1867 : }
1868 :
1869 : /* ---- Return success ---- */
1870 2 : return CE_None;
1871 : }
1872 :
1873 : // Given nt intervals spaced by dt and beginning at t0, return the index of
1874 : // the lower bound of the interval that should be used to
1875 : // interpolate/extrapolate a value for t.
1876 17 : static std::size_t intervalLeft(double t0, double dt, std::size_t nt, double t)
1877 : {
1878 17 : if (t < t0)
1879 : {
1880 4 : return 0;
1881 : }
1882 :
1883 13 : std::size_t n = static_cast<std::size_t>((t - t0) / dt);
1884 :
1885 13 : if (n >= nt - 1)
1886 : {
1887 3 : return nt - 2;
1888 : }
1889 :
1890 10 : return n;
1891 : }
1892 :
1893 17 : static double InterpolateLinear(double dfX0, double dfX1, double dfY0,
1894 : double dfY1, double dfX)
1895 : {
1896 17 : return dfY0 + (dfX - dfX0) * (dfY1 - dfY0) / (dfX1 - dfX0);
1897 : }
1898 :
1899 13 : static double InterpolateExponential(double dfX0, double dfX1, double dfY0,
1900 : double dfY1, double dfX)
1901 : {
1902 13 : const double r = std::log(dfY1 / dfY0) / (dfX1 - dfX0);
1903 13 : return dfY0 * std::exp(r * (dfX - dfX0));
1904 : }
1905 :
1906 : static const char pszInterpolatePixelFuncMetadata[] =
1907 : "<PixelFunctionArgumentsList>"
1908 : " <Argument name='t0' description='t0' type='double' mandatory='1' />"
1909 : " <Argument name='dt' description='dt' type='double' mandatory='1' />"
1910 : " <Argument name='t' description='t' type='double' mandatory='1' />"
1911 : " <Argument type='builtin' value='NoData' optional='true' />"
1912 : "</PixelFunctionArgumentsList>";
1913 :
1914 : template <decltype(InterpolateLinear) InterpolationFunction>
1915 17 : CPLErr InterpolatePixelFunc(void **papoSources, int nSources, void *pData,
1916 : int nXSize, int nYSize, GDALDataType eSrcType,
1917 : GDALDataType eBufType, int nPixelSpace,
1918 : int nLineSpace, CSLConstList papszArgs)
1919 : {
1920 : /* ---- Init ---- */
1921 17 : if (GDALDataTypeIsComplex(eSrcType))
1922 0 : return CE_Failure;
1923 :
1924 : double dfT0;
1925 17 : if (FetchDoubleArg(papszArgs, "t0", &dfT0) == CE_Failure)
1926 0 : return CE_Failure;
1927 :
1928 : double dfT;
1929 17 : if (FetchDoubleArg(papszArgs, "t", &dfT) == CE_Failure)
1930 0 : return CE_Failure;
1931 :
1932 : double dfDt;
1933 17 : if (FetchDoubleArg(papszArgs, "dt", &dfDt) == CE_Failure)
1934 0 : return CE_Failure;
1935 :
1936 17 : double dfNoData{0};
1937 17 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
1938 17 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
1939 0 : return CE_Failure;
1940 :
1941 17 : if (nSources < 2)
1942 : {
1943 0 : CPLError(CE_Failure, CPLE_AppDefined,
1944 : "At least two sources required for interpolation.");
1945 0 : return CE_Failure;
1946 : }
1947 :
1948 17 : if (dfT == 0 || !std::isfinite(dfT))
1949 : {
1950 0 : CPLError(CE_Failure, CPLE_AppDefined, "dt must be finite and non-zero");
1951 0 : return CE_Failure;
1952 : }
1953 :
1954 17 : const auto i0 = intervalLeft(dfT0, dfDt, nSources, dfT);
1955 17 : const auto i1 = i0 + 1;
1956 17 : const double dfX0 = dfT0 + static_cast<double>(i0) * dfDt;
1957 17 : const double dfX1 = dfT0 + static_cast<double>(i0 + 1) * dfDt;
1958 :
1959 : /* ---- Set pixels ---- */
1960 17 : size_t ii = 0;
1961 41 : for (int iLine = 0; iLine < nYSize; ++iLine)
1962 : {
1963 72 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1964 : {
1965 48 : const double dfY0 = GetSrcVal(papoSources[i0], eSrcType, ii);
1966 48 : const double dfY1 = GetSrcVal(papoSources[i1], eSrcType, ii);
1967 :
1968 48 : double dfPixVal = dfNoData;
1969 48 : if (dfT == dfX0)
1970 8 : dfPixVal = dfY0;
1971 40 : else if (dfT == dfX1)
1972 0 : dfPixVal = dfY1;
1973 52 : else if (!bHasNoData ||
1974 12 : (!IsNoData(dfY0, dfNoData) && !IsNoData(dfY1, dfNoData)))
1975 30 : dfPixVal = InterpolationFunction(dfX0, dfX1, dfY0, dfY1, dfT);
1976 :
1977 48 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1978 : static_cast<GByte *>(pData) +
1979 48 : static_cast<GSpacing>(nLineSpace) * iLine +
1980 48 : iCol * nPixelSpace,
1981 : eBufType, nPixelSpace, 1);
1982 : }
1983 : }
1984 :
1985 : /* ---- Return success ---- */
1986 17 : return CE_None;
1987 : }
1988 :
1989 : static const char pszReplaceNoDataPixelFuncMetadata[] =
1990 : "<PixelFunctionArgumentsList>"
1991 : " <Argument type='builtin' value='NoData' />"
1992 : " <Argument name='to' type='double' description='New NoData value to be "
1993 : "replaced' default='nan' />"
1994 : "</PixelFunctionArgumentsList>";
1995 :
1996 2 : static CPLErr ReplaceNoDataPixelFunc(void **papoSources, int nSources,
1997 : void *pData, int nXSize, int nYSize,
1998 : GDALDataType eSrcType,
1999 : GDALDataType eBufType, int nPixelSpace,
2000 : int nLineSpace, CSLConstList papszArgs)
2001 : {
2002 : /* ---- Init ---- */
2003 2 : if (nSources != 1)
2004 0 : return CE_Failure;
2005 2 : if (GDALDataTypeIsComplex(eSrcType))
2006 : {
2007 0 : CPLError(CE_Failure, CPLE_AppDefined,
2008 : "replace_nodata cannot convert complex data types");
2009 0 : return CE_Failure;
2010 : }
2011 :
2012 2 : double dfOldNoData, dfNewNoData = NAN;
2013 2 : if (FetchDoubleArg(papszArgs, "NoData", &dfOldNoData) != CE_None)
2014 0 : return CE_Failure;
2015 2 : if (FetchDoubleArg(papszArgs, "to", &dfNewNoData, &dfNewNoData) != CE_None)
2016 0 : return CE_Failure;
2017 :
2018 2 : if (!GDALDataTypeIsFloating(eBufType) && std::isnan(dfNewNoData))
2019 : {
2020 0 : CPLError(CE_Failure, CPLE_AppDefined,
2021 : "Using nan requires a floating point type output buffer");
2022 0 : return CE_Failure;
2023 : }
2024 :
2025 : /* ---- Set pixels ---- */
2026 2 : size_t ii = 0;
2027 102 : for (int iLine = 0; iLine < nYSize; ++iLine)
2028 : {
2029 5100 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
2030 : {
2031 5000 : double dfPixVal = GetSrcVal(papoSources[0], eSrcType, ii);
2032 5000 : if (dfPixVal == dfOldNoData || std::isnan(dfPixVal))
2033 3200 : dfPixVal = dfNewNoData;
2034 :
2035 5000 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
2036 : static_cast<GByte *>(pData) +
2037 5000 : static_cast<GSpacing>(nLineSpace) * iLine +
2038 5000 : iCol * nPixelSpace,
2039 : eBufType, nPixelSpace, 1);
2040 : }
2041 : }
2042 :
2043 : /* ---- Return success ---- */
2044 2 : return CE_None;
2045 : }
2046 :
2047 : static const char pszScalePixelFuncMetadata[] =
2048 : "<PixelFunctionArgumentsList>"
2049 : " <Argument type='builtin' value='offset' />"
2050 : " <Argument type='builtin' value='scale' />"
2051 : " <Argument type='builtin' value='NoData' optional='true' />"
2052 : "</PixelFunctionArgumentsList>";
2053 :
2054 2 : static CPLErr ScalePixelFunc(void **papoSources, int nSources, void *pData,
2055 : int nXSize, int nYSize, GDALDataType eSrcType,
2056 : GDALDataType eBufType, int nPixelSpace,
2057 : int nLineSpace, CSLConstList papszArgs)
2058 : {
2059 : /* ---- Init ---- */
2060 2 : if (nSources != 1)
2061 0 : return CE_Failure;
2062 2 : if (GDALDataTypeIsComplex(eSrcType))
2063 : {
2064 0 : CPLError(CE_Failure, CPLE_AppDefined,
2065 : "scale cannot by applied to complex data types");
2066 0 : return CE_Failure;
2067 : }
2068 :
2069 2 : double dfNoData{0};
2070 2 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
2071 2 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
2072 0 : return CE_Failure;
2073 :
2074 : double dfScale, dfOffset;
2075 2 : if (FetchDoubleArg(papszArgs, "scale", &dfScale) != CE_None)
2076 0 : return CE_Failure;
2077 2 : if (FetchDoubleArg(papszArgs, "offset", &dfOffset) != CE_None)
2078 0 : return CE_Failure;
2079 :
2080 : /* ---- Set pixels ---- */
2081 2 : size_t ii = 0;
2082 23 : for (int iLine = 0; iLine < nYSize; ++iLine)
2083 : {
2084 423 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
2085 : {
2086 402 : const double dfVal = GetSrcVal(papoSources[0], eSrcType, ii);
2087 :
2088 2 : const double dfPixVal = bHasNoData && IsNoData(dfVal, dfNoData)
2089 404 : ? dfNoData
2090 400 : : dfVal * dfScale + dfOffset;
2091 :
2092 402 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
2093 : static_cast<GByte *>(pData) +
2094 402 : static_cast<GSpacing>(nLineSpace) * iLine +
2095 402 : iCol * nPixelSpace,
2096 : eBufType, nPixelSpace, 1);
2097 : }
2098 : }
2099 :
2100 : /* ---- Return success ---- */
2101 2 : return CE_None;
2102 : }
2103 :
2104 : static const char pszNormDiffPixelFuncMetadata[] =
2105 : "<PixelFunctionArgumentsList>"
2106 : " <Argument type='builtin' value='NoData' optional='true' />"
2107 : "</PixelFunctionArgumentsList>";
2108 :
2109 3 : static CPLErr NormDiffPixelFunc(void **papoSources, int nSources, void *pData,
2110 : int nXSize, int nYSize, GDALDataType eSrcType,
2111 : GDALDataType eBufType, int nPixelSpace,
2112 : int nLineSpace, CSLConstList papszArgs)
2113 : {
2114 : /* ---- Init ---- */
2115 3 : if (nSources != 2)
2116 0 : return CE_Failure;
2117 :
2118 3 : if (GDALDataTypeIsComplex(eSrcType))
2119 : {
2120 0 : CPLError(CE_Failure, CPLE_AppDefined,
2121 : "norm_diff cannot by applied to complex data types");
2122 0 : return CE_Failure;
2123 : }
2124 :
2125 3 : double dfNoData{0};
2126 3 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
2127 3 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
2128 0 : return CE_Failure;
2129 :
2130 : /* ---- Set pixels ---- */
2131 3 : size_t ii = 0;
2132 11 : for (int iLine = 0; iLine < nYSize; ++iLine)
2133 : {
2134 42 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
2135 : {
2136 34 : const double dfLeftVal = GetSrcVal(papoSources[0], eSrcType, ii);
2137 34 : const double dfRightVal = GetSrcVal(papoSources[1], eSrcType, ii);
2138 :
2139 34 : double dfPixVal = dfNoData;
2140 :
2141 35 : if (!bHasNoData || (!IsNoData(dfLeftVal, dfNoData) &&
2142 1 : !IsNoData(dfRightVal, dfNoData)))
2143 : {
2144 30 : const double dfDenom = (dfLeftVal + dfRightVal);
2145 : // coverity[divide_by_zero]
2146 30 : dfPixVal = dfDenom == 0
2147 30 : ? std::numeric_limits<double>::infinity()
2148 30 : : (dfLeftVal - dfRightVal) / dfDenom;
2149 : }
2150 :
2151 34 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
2152 : static_cast<GByte *>(pData) +
2153 34 : static_cast<GSpacing>(nLineSpace) * iLine +
2154 34 : iCol * nPixelSpace,
2155 : eBufType, nPixelSpace, 1);
2156 : }
2157 : }
2158 :
2159 : /* ---- Return success ---- */
2160 3 : return CE_None;
2161 : } // NormDiffPixelFunc
2162 :
2163 : /************************************************************************/
2164 : /* pszMinMaxFuncMetadataNodata */
2165 : /************************************************************************/
2166 :
2167 : static const char pszMinMaxFuncMetadataNodata[] =
2168 : "<PixelFunctionArgumentsList>"
2169 : " <Argument type='builtin' value='NoData' optional='true' />"
2170 : " <Argument name='propagateNoData' description='Whether the output value "
2171 : "should be NoData as as soon as one source is NoData' type='boolean' "
2172 : "default='false' />"
2173 : "</PixelFunctionArgumentsList>";
2174 :
2175 : template <class Comparator>
2176 12 : static CPLErr MinOrMaxPixelFunc(void **papoSources, int nSources, void *pData,
2177 : int nXSize, int nYSize, GDALDataType eSrcType,
2178 : GDALDataType eBufType, int nPixelSpace,
2179 : int nLineSpace, CSLConstList papszArgs)
2180 : {
2181 : /* ---- Init ---- */
2182 12 : if (GDALDataTypeIsComplex(eSrcType))
2183 : {
2184 0 : CPLError(CE_Failure, CPLE_AppDefined,
2185 : "Complex data type not supported for min/max().");
2186 0 : return CE_Failure;
2187 : }
2188 :
2189 12 : double dfNoData = std::numeric_limits<double>::quiet_NaN();
2190 12 : if (FetchDoubleArg(papszArgs, "NoData", &dfNoData, &dfNoData) != CE_None)
2191 0 : return CE_Failure;
2192 12 : const bool bPropagateNoData = CPLTestBool(
2193 : CSLFetchNameValueDef(papszArgs, "propagateNoData", "false"));
2194 :
2195 : /* ---- Set pixels ---- */
2196 12 : size_t ii = 0;
2197 269 : for (int iLine = 0; iLine < nYSize; ++iLine)
2198 : {
2199 12773 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
2200 : {
2201 12516 : double dfRes = std::numeric_limits<double>::quiet_NaN();
2202 :
2203 43445 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
2204 : {
2205 32985 : const double dfVal = GetSrcVal(papoSources[iSrc], eSrcType, ii);
2206 :
2207 32985 : if (std::isnan(dfVal) || dfVal == dfNoData)
2208 : {
2209 14370 : if (bPropagateNoData)
2210 : {
2211 2056 : dfRes = dfNoData;
2212 2056 : break;
2213 : }
2214 : }
2215 18615 : else if (Comparator::compare(dfVal, dfRes))
2216 : {
2217 9961 : dfRes = dfVal;
2218 : }
2219 : }
2220 :
2221 12516 : if (!bPropagateNoData && std::isnan(dfRes))
2222 : {
2223 3203 : dfRes = dfNoData;
2224 : }
2225 :
2226 12516 : GDALCopyWords(&dfRes, GDT_Float64, 0,
2227 : static_cast<GByte *>(pData) +
2228 12516 : static_cast<GSpacing>(nLineSpace) * iLine +
2229 12516 : iCol * nPixelSpace,
2230 : eBufType, nPixelSpace, 1);
2231 : }
2232 : }
2233 :
2234 : /* ---- Return success ---- */
2235 12 : return CE_None;
2236 : } /* MinOrMaxPixelFunc */
2237 :
2238 7 : static CPLErr MinPixelFunc(void **papoSources, int nSources, void *pData,
2239 : int nXSize, int nYSize, GDALDataType eSrcType,
2240 : GDALDataType eBufType, int nPixelSpace,
2241 : int nLineSpace, CSLConstList papszArgs)
2242 : {
2243 : struct Comparator
2244 : {
2245 2712 : static bool compare(double x, double resVal)
2246 : {
2247 : // Written this way to deal with resVal being NaN
2248 2712 : return !(x >= resVal);
2249 : }
2250 : };
2251 :
2252 7 : return MinOrMaxPixelFunc<Comparator>(papoSources, nSources, pData, nXSize,
2253 : nYSize, eSrcType, eBufType,
2254 7 : nPixelSpace, nLineSpace, papszArgs);
2255 : }
2256 :
2257 5 : static CPLErr MaxPixelFunc(void **papoSources, int nSources, void *pData,
2258 : int nXSize, int nYSize, GDALDataType eSrcType,
2259 : GDALDataType eBufType, int nPixelSpace,
2260 : int nLineSpace, CSLConstList papszArgs)
2261 : {
2262 : struct Comparator
2263 : {
2264 15903 : static bool compare(double x, double resVal)
2265 : {
2266 : // Written this way to deal with resVal being NaN
2267 15903 : return !(x <= resVal);
2268 : }
2269 : };
2270 :
2271 5 : return MinOrMaxPixelFunc<Comparator>(papoSources, nSources, pData, nXSize,
2272 : nYSize, eSrcType, eBufType,
2273 5 : nPixelSpace, nLineSpace, papszArgs);
2274 : }
2275 :
2276 : static const char pszExprPixelFuncMetadata[] =
2277 : "<PixelFunctionArgumentsList>"
2278 : " <Argument name='expression' "
2279 : " description='Expression to be evaluated' "
2280 : " type='string'></Argument>"
2281 : " <Argument name='dialect' "
2282 : " description='Expression dialect' "
2283 : " type='string-select'"
2284 : " default='muparser'>"
2285 : " <Value>exprtk</Value>"
2286 : " <Value>muparser</Value>"
2287 : " </Argument>"
2288 : " <Argument type='builtin' value='source_names' />"
2289 : "</PixelFunctionArgumentsList>";
2290 :
2291 60 : static CPLErr ExprPixelFunc(void **papoSources, int nSources, void *pData,
2292 : int nXSize, int nYSize, GDALDataType eSrcType,
2293 : GDALDataType eBufType, int nPixelSpace,
2294 : int nLineSpace, CSLConstList papszArgs)
2295 : {
2296 : /* ---- Init ---- */
2297 :
2298 60 : if (GDALDataTypeIsComplex(eSrcType))
2299 : {
2300 0 : CPLError(CE_Failure, CPLE_AppDefined,
2301 : "expression cannot by applied to complex data types");
2302 0 : return CE_Failure;
2303 : }
2304 :
2305 60 : std::unique_ptr<gdal::MathExpression> poExpression;
2306 :
2307 60 : const char *pszExpression = CSLFetchNameValue(papszArgs, "expression");
2308 :
2309 60 : const char *pszSourceNames = CSLFetchNameValue(papszArgs, "source_names");
2310 : const CPLStringList aosSourceNames(
2311 120 : CSLTokenizeString2(pszSourceNames, "|", 0));
2312 :
2313 120 : std::vector<double> adfValuesForPixel(nSources);
2314 :
2315 60 : const char *pszDialect = CSLFetchNameValue(papszArgs, "dialect");
2316 60 : if (!pszDialect)
2317 : {
2318 0 : pszDialect = "muparser";
2319 : }
2320 :
2321 60 : poExpression = gdal::MathExpression::Create(pszExpression, pszDialect);
2322 :
2323 : // cppcheck-suppress knownConditionTrueFalse
2324 60 : if (!poExpression)
2325 : {
2326 0 : return CE_Failure;
2327 : }
2328 :
2329 : {
2330 60 : int iSource = 0;
2331 181 : for (const auto &osName : aosSourceNames)
2332 : {
2333 242 : poExpression->RegisterVariable(osName,
2334 121 : &adfValuesForPixel[iSource++]);
2335 : }
2336 : }
2337 :
2338 60 : if (strstr(pszExpression, "BANDS"))
2339 : {
2340 2 : poExpression->RegisterVector("BANDS", &adfValuesForPixel);
2341 : }
2342 :
2343 : std::unique_ptr<double, VSIFreeReleaser> padfResults(
2344 120 : static_cast<double *>(VSI_MALLOC2_VERBOSE(nXSize, sizeof(double))));
2345 60 : if (!padfResults)
2346 0 : return CE_Failure;
2347 :
2348 : /* ---- Set pixels ---- */
2349 60 : size_t ii = 0;
2350 598 : for (int iLine = 0; iLine < nYSize; ++iLine)
2351 : {
2352 19278 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
2353 : {
2354 67602 : for (int iSrc = 0; iSrc < nSources; iSrc++)
2355 : {
2356 : // cppcheck-suppress unreadVariable
2357 48862 : adfValuesForPixel[iSrc] =
2358 48862 : GetSrcVal(papoSources[iSrc], eSrcType, ii);
2359 : }
2360 :
2361 18740 : if (auto eErr = poExpression->Evaluate(); eErr != CE_None)
2362 : {
2363 12 : return CE_Failure;
2364 : }
2365 : else
2366 : {
2367 18728 : padfResults.get()[iCol] = poExpression->Results()[0];
2368 : }
2369 : }
2370 :
2371 538 : GDALCopyWords(padfResults.get(), GDT_Float64, sizeof(double),
2372 : static_cast<GByte *>(pData) +
2373 538 : static_cast<GSpacing>(nLineSpace) * iLine,
2374 : eBufType, nPixelSpace, nXSize);
2375 : }
2376 :
2377 : /* ---- Return success ---- */
2378 48 : return CE_None;
2379 : } // ExprPixelFunc
2380 :
2381 : static const char pszReclassifyPixelFuncMetadata[] =
2382 : "<PixelFunctionArgumentsList>"
2383 : " <Argument name='mapping' "
2384 : " description='Lookup table for mapping, in format "
2385 : "from=to,from=to' "
2386 : " type='string'></Argument>"
2387 : " <Argument type='builtin' value='NoData' optional='true' />"
2388 : "</PixelFunctionArgumentsList>";
2389 :
2390 33 : static CPLErr ReclassifyPixelFunc(void **papoSources, int nSources, void *pData,
2391 : int nXSize, int nYSize, GDALDataType eSrcType,
2392 : GDALDataType eBufType, int nPixelSpace,
2393 : int nLineSpace, CSLConstList papszArgs)
2394 : {
2395 33 : if (GDALDataTypeIsComplex(eSrcType))
2396 : {
2397 0 : CPLError(CE_Failure, CPLE_AppDefined,
2398 : "reclassify cannot by applied to complex data types");
2399 0 : return CE_Failure;
2400 : }
2401 :
2402 33 : if (nSources != 1)
2403 : {
2404 0 : CPLError(CE_Failure, CPLE_AppDefined,
2405 : "reclassify only be applied to a single source at a time");
2406 0 : return CE_Failure;
2407 : }
2408 33 : std::optional<double> noDataValue{};
2409 :
2410 33 : const char *pszNoData = CSLFetchNameValue(papszArgs, "NoData");
2411 33 : if (pszNoData != nullptr)
2412 : {
2413 10 : noDataValue = CPLAtof(pszNoData);
2414 : }
2415 :
2416 33 : const char *pszMappings = CSLFetchNameValue(papszArgs, "mapping");
2417 33 : if (pszMappings == nullptr)
2418 : {
2419 0 : CPLError(CE_Failure, CPLE_AppDefined,
2420 : "reclassify must be called with 'mapping' argument");
2421 0 : return CE_Failure;
2422 : }
2423 :
2424 66 : gdal::Reclassifier oReclassifier;
2425 33 : if (auto eErr = oReclassifier.Init(pszMappings, noDataValue, eBufType);
2426 : eErr != CE_None)
2427 : {
2428 14 : return eErr;
2429 : }
2430 :
2431 : std::unique_ptr<double, VSIFreeReleaser> padfResults(
2432 38 : static_cast<double *>(VSI_MALLOC2_VERBOSE(nXSize, sizeof(double))));
2433 19 : if (!padfResults)
2434 0 : return CE_Failure;
2435 :
2436 19 : size_t ii = 0;
2437 19 : bool bSuccess = false;
2438 435 : for (int iLine = 0; iLine < nYSize; ++iLine)
2439 : {
2440 20805 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
2441 : {
2442 20389 : double srcVal = GetSrcVal(papoSources[0], eSrcType, ii);
2443 40778 : padfResults.get()[iCol] =
2444 20389 : oReclassifier.Reclassify(srcVal, bSuccess);
2445 20389 : if (!bSuccess)
2446 : {
2447 2 : CPLError(CE_Failure, CPLE_AppDefined,
2448 : "Encountered value %g with no specified mapping",
2449 : srcVal);
2450 2 : return CE_Failure;
2451 : }
2452 : }
2453 :
2454 416 : GDALCopyWords(padfResults.get(), GDT_Float64, sizeof(double),
2455 : static_cast<GByte *>(pData) +
2456 416 : static_cast<GSpacing>(nLineSpace) * iLine,
2457 : eBufType, nPixelSpace, nXSize);
2458 : }
2459 :
2460 17 : return CE_None;
2461 : } // ReclassifyPixelFunc
2462 :
2463 : struct MeanKernel
2464 : {
2465 : static constexpr const char *pszName = "mean";
2466 :
2467 : double dfSum = 0;
2468 : int nValidSources = 0;
2469 :
2470 6 : void Reset()
2471 : {
2472 6 : dfSum = 0;
2473 6 : nValidSources = 0;
2474 6 : }
2475 :
2476 3 : static CPLErr ProcessArguments(CSLConstList)
2477 : {
2478 3 : return CE_None;
2479 : }
2480 :
2481 3 : void ProcessPixel(double dfVal)
2482 : {
2483 3 : dfSum += dfVal;
2484 3 : nValidSources++;
2485 3 : }
2486 :
2487 4 : bool HasValue() const
2488 : {
2489 4 : return nValidSources > 0;
2490 : }
2491 :
2492 1 : double GetValue() const
2493 : {
2494 1 : return dfSum / nValidSources;
2495 : }
2496 : };
2497 :
2498 : struct GeoMeanKernel
2499 : {
2500 : static constexpr const char *pszName = "geometric_mean";
2501 :
2502 : double dfProduct = 1;
2503 : int nValidSources = 0;
2504 :
2505 6 : void Reset()
2506 : {
2507 6 : dfProduct = 1;
2508 6 : nValidSources = 0;
2509 6 : }
2510 :
2511 3 : static CPLErr ProcessArguments(CSLConstList)
2512 : {
2513 3 : return CE_None;
2514 : }
2515 :
2516 3 : void ProcessPixel(double dfVal)
2517 : {
2518 3 : dfProduct *= dfVal;
2519 3 : nValidSources++;
2520 3 : }
2521 :
2522 4 : bool HasValue() const
2523 : {
2524 4 : return nValidSources > 0;
2525 : }
2526 :
2527 1 : double GetValue() const
2528 : {
2529 1 : return std::pow(dfProduct, 1.0 / nValidSources);
2530 : }
2531 : };
2532 :
2533 : struct HarmonicMeanKernel
2534 : {
2535 : static constexpr const char *pszName = "harmonic_mean";
2536 :
2537 : double dfDenom = 0;
2538 : int nValidSources = 0;
2539 : bool bValueIsZero = false;
2540 : bool bPropagateZero = false;
2541 :
2542 10 : void Reset()
2543 : {
2544 10 : dfDenom = 0;
2545 10 : nValidSources = 0;
2546 10 : bValueIsZero = false;
2547 10 : }
2548 :
2549 7 : void ProcessPixel(double dfVal)
2550 : {
2551 7 : if (dfVal == 0)
2552 : {
2553 2 : bValueIsZero = true;
2554 : }
2555 : else
2556 : {
2557 5 : dfDenom += 1 / dfVal;
2558 : }
2559 7 : nValidSources++;
2560 7 : }
2561 :
2562 5 : CPLErr ProcessArguments(CSLConstList papszArgs)
2563 : {
2564 5 : bPropagateZero =
2565 5 : CPLTestBool(CSLFetchNameValueDef(papszArgs, "propagateZero", "0"));
2566 5 : return CE_None;
2567 : }
2568 :
2569 8 : bool HasValue() const
2570 : {
2571 8 : return dfDenom > 0 && (bPropagateZero || !bValueIsZero);
2572 : }
2573 :
2574 2 : double GetValue() const
2575 : {
2576 2 : if (bPropagateZero && bValueIsZero)
2577 : {
2578 1 : return 0;
2579 : }
2580 1 : return static_cast<double>(nValidSources) / dfDenom;
2581 : }
2582 : };
2583 :
2584 : struct MedianKernel
2585 : {
2586 : static constexpr const char *pszName = "median";
2587 :
2588 : mutable std::vector<double> values{};
2589 :
2590 8 : void Reset()
2591 : {
2592 8 : values.clear();
2593 8 : }
2594 :
2595 4 : static CPLErr ProcessArguments(CSLConstList)
2596 : {
2597 4 : return CE_None;
2598 : }
2599 :
2600 6 : void ProcessPixel(double dfVal)
2601 : {
2602 6 : if (!std::isnan(dfVal))
2603 : {
2604 6 : values.push_back(dfVal);
2605 : }
2606 6 : }
2607 :
2608 6 : bool HasValue() const
2609 : {
2610 6 : return !values.empty();
2611 : }
2612 :
2613 2 : double GetValue() const
2614 : {
2615 2 : std::sort(values.begin(), values.end());
2616 2 : if (values.size() % 2 == 0)
2617 : {
2618 : return 0.5 *
2619 1 : (values[values.size() / 2 - 1] + values[values.size() / 2]);
2620 : }
2621 :
2622 1 : return values[values.size() / 2];
2623 : }
2624 : };
2625 :
2626 : struct ModeKernel
2627 : {
2628 : static constexpr const char *pszName = "mode";
2629 :
2630 : std::map<double, size_t> counts{};
2631 : std::size_t nanCount{0};
2632 : double dfMax = std::numeric_limits<double>::quiet_NaN();
2633 : decltype(counts.begin()) oMax = counts.end();
2634 :
2635 6 : void Reset()
2636 : {
2637 6 : nanCount = 0;
2638 6 : counts.clear();
2639 6 : oMax = counts.end();
2640 6 : }
2641 :
2642 3 : static CPLErr ProcessArguments(CSLConstList)
2643 : {
2644 3 : return CE_None;
2645 : }
2646 :
2647 8 : void ProcessPixel(double dfVal)
2648 : {
2649 8 : if (std::isnan(dfVal))
2650 : {
2651 2 : nanCount += 1;
2652 2 : return;
2653 : }
2654 :
2655 : // if dfVal is NaN, try_emplace will return an entry for a different key!
2656 6 : auto [it, inserted] = counts.try_emplace(dfVal, 0);
2657 :
2658 6 : it->second += 1;
2659 :
2660 6 : if (oMax == counts.end() || it->second > oMax->second)
2661 : {
2662 4 : oMax = it;
2663 : }
2664 : }
2665 :
2666 4 : bool HasValue() const
2667 : {
2668 4 : return nanCount > 0 || oMax != counts.end();
2669 : }
2670 :
2671 2 : double GetValue() const
2672 : {
2673 2 : double ret = std::numeric_limits<double>::quiet_NaN();
2674 2 : if (oMax != counts.end())
2675 : {
2676 2 : const size_t nCount = oMax->second;
2677 2 : if (nCount > nanCount)
2678 1 : ret = oMax->first;
2679 : }
2680 2 : return ret;
2681 : }
2682 : };
2683 :
2684 : static const char pszBasicPixelFuncMetadata[] =
2685 : "<PixelFunctionArgumentsList>"
2686 : " <Argument type='builtin' value='NoData' optional='true' />"
2687 : " <Argument name='propagateNoData' description='Whether the output value "
2688 : "should be NoData as as soon as one source is NoData' type='boolean' "
2689 : "default='false' />"
2690 : "</PixelFunctionArgumentsList>";
2691 :
2692 : template <typename T>
2693 18 : static CPLErr BasicPixelFunc(void **papoSources, int nSources, void *pData,
2694 : int nXSize, int nYSize, GDALDataType eSrcType,
2695 : GDALDataType eBufType, int nPixelSpace,
2696 : int nLineSpace, CSLConstList papszArgs)
2697 : {
2698 : /* ---- Init ---- */
2699 25 : T oKernel;
2700 :
2701 18 : if (GDALDataTypeIsComplex(eSrcType))
2702 : {
2703 0 : CPLError(CE_Failure, CPLE_AppDefined,
2704 : "Complex data types not supported by %s", oKernel.pszName);
2705 0 : return CE_Failure;
2706 : }
2707 :
2708 18 : double dfNoData{0};
2709 18 : const bool bHasNoData = CSLFindName(papszArgs, "NoData") != -1;
2710 18 : if (bHasNoData && FetchDoubleArg(papszArgs, "NoData", &dfNoData) != CE_None)
2711 0 : return CE_Failure;
2712 :
2713 18 : const bool bPropagateNoData = CPLTestBool(
2714 : CSLFetchNameValueDef(papszArgs, "propagateNoData", "false"));
2715 :
2716 18 : if (oKernel.ProcessArguments(papszArgs) == CE_Failure)
2717 : {
2718 0 : return CE_Failure;
2719 : }
2720 :
2721 : /* ---- Set pixels ---- */
2722 18 : size_t ii = 0;
2723 36 : for (int iLine = 0; iLine < nYSize; ++iLine)
2724 : {
2725 54 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
2726 : {
2727 36 : oKernel.Reset();
2728 36 : bool bWriteNoData = false;
2729 :
2730 127 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
2731 : {
2732 101 : const double dfVal = GetSrcVal(papoSources[iSrc], eSrcType, ii);
2733 :
2734 101 : if (bHasNoData && IsNoData(dfVal, dfNoData))
2735 : {
2736 74 : if (bPropagateNoData)
2737 : {
2738 10 : bWriteNoData = true;
2739 10 : break;
2740 : }
2741 : }
2742 : else
2743 : {
2744 27 : oKernel.ProcessPixel(dfVal);
2745 : }
2746 : }
2747 :
2748 36 : double dfPixVal{dfNoData};
2749 36 : if (!bWriteNoData && oKernel.HasValue())
2750 : {
2751 8 : dfPixVal = oKernel.GetValue();
2752 : }
2753 :
2754 36 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
2755 : static_cast<GByte *>(pData) +
2756 36 : static_cast<GSpacing>(nLineSpace) * iLine +
2757 36 : iCol * nPixelSpace,
2758 : eBufType, nPixelSpace, 1);
2759 : }
2760 : }
2761 :
2762 : /* ---- Return success ---- */
2763 18 : return CE_None;
2764 : } // BasicPixelFunc
2765 :
2766 : /************************************************************************/
2767 : /* GDALRegisterDefaultPixelFunc() */
2768 : /************************************************************************/
2769 :
2770 : /**
2771 : * This adds a default set of pixel functions to the global list of
2772 : * available pixel functions for derived bands:
2773 : *
2774 : * - "real": extract real part from a single raster band (just a copy if the
2775 : * input is non-complex)
2776 : * - "imag": extract imaginary part from a single raster band (0 for
2777 : * non-complex)
2778 : * - "complex": make a complex band merging two bands used as real and
2779 : * imag values
2780 : * - "polar": make a complex band using input bands for amplitude and
2781 : * phase values (b1 * exp( j * b2 ))
2782 : * - "mod": extract module from a single raster band (real or complex)
2783 : * - "phase": extract phase from a single raster band [-PI,PI] (0 or PI for
2784 : non-complex)
2785 : * - "conj": computes the complex conjugate of a single raster band (just a
2786 : * copy if the input is non-complex)
2787 : * - "sum": sum 2 or more raster bands
2788 : * - "diff": computes the difference between 2 raster bands (b1 - b2)
2789 : * - "mul": multiply 2 or more raster bands
2790 : * - "div": divide one raster band by another (b1 / b2).
2791 : * - "min": minimum value of 2 or more raster bands
2792 : * - "max": maximum value of 2 or more raster bands
2793 : * - "norm_diff": computes the normalized difference between two raster bands:
2794 : * ``(b1 - b2)/(b1 + b2)``.
2795 : * - "cmul": multiply the first band for the complex conjugate of the second
2796 : * - "inv": inverse (1./x).
2797 : * - "intensity": computes the intensity Re(x*conj(x)) of a single raster band
2798 : * (real or complex)
2799 : * - "sqrt": perform the square root of a single raster band (real only)
2800 : * - "log10": compute the logarithm (base 10) of the abs of a single raster
2801 : * band (real or complex): log10( abs( x ) )
2802 : * - "dB": perform conversion to dB of the abs of a single raster
2803 : * band (real or complex): 20. * log10( abs( x ) ).
2804 : * Note: the optional fact parameter can be set to 10. to get the
2805 : * alternative formula: 10. * log10( abs( x ) )
2806 : * - "exp": computes the exponential of each element in the input band ``x``
2807 : * (of real values): ``e ^ x``.
2808 : * The function also accepts two optional parameters: ``base`` and
2809 : ``fact``
2810 : * that allow to compute the generalized formula: ``base ^ ( fact *
2811 : x)``.
2812 : * Note: this function is the recommended one to perform conversion
2813 : * form logarithmic scale (dB): `` 10. ^ (x / 20.)``, in this case
2814 : * ``base = 10.`` and ``fact = 1./20``
2815 : * - "dB2amp": perform scale conversion from logarithmic to linear
2816 : * (amplitude) (i.e. 10 ^ ( x / 20 ) ) of a single raster
2817 : * band (real only).
2818 : * Deprecated in GDAL v3.5. Please use the ``exp`` pixel function
2819 : with
2820 : * ``base = 10.`` and ``fact = 0.05`` i.e. ``1./20``
2821 : * - "dB2pow": perform scale conversion from logarithmic to linear
2822 : * (power) (i.e. 10 ^ ( x / 10 ) ) of a single raster
2823 : * band (real only)
2824 : * Deprecated in GDAL v3.5. Please use the ``exp`` pixel function
2825 : with
2826 : * ``base = 10.`` and ``fact = 0.1`` i.e. ``1./10``
2827 : * - "pow": raise a single raster band to a constant power
2828 : * - "interpolate_linear": interpolate values between two raster bands
2829 : * using linear interpolation
2830 : * - "interpolate_exp": interpolate values between two raster bands using
2831 : * exponential interpolation
2832 : * - "scale": Apply the RasterBand metadata values of "offset" and "scale"
2833 : * - "reclassify": Reclassify values matching ranges in a table
2834 : * - "nan": Convert incoming NoData values to IEEE 754 nan
2835 : *
2836 : * @see GDALAddDerivedBandPixelFunc
2837 : *
2838 : * @return CE_None
2839 : */
2840 1407 : CPLErr GDALRegisterDefaultPixelFunc()
2841 : {
2842 1407 : GDALAddDerivedBandPixelFunc("real", RealPixelFunc);
2843 1407 : GDALAddDerivedBandPixelFunc("imag", ImagPixelFunc);
2844 1407 : GDALAddDerivedBandPixelFunc("complex", ComplexPixelFunc);
2845 1407 : GDALAddDerivedBandPixelFuncWithArgs("polar", PolarPixelFunc,
2846 : pszPolarPixelFuncMetadata);
2847 1407 : GDALAddDerivedBandPixelFunc("mod", ModulePixelFunc);
2848 1407 : GDALAddDerivedBandPixelFunc("phase", PhasePixelFunc);
2849 1407 : GDALAddDerivedBandPixelFunc("conj", ConjPixelFunc);
2850 1407 : GDALAddDerivedBandPixelFuncWithArgs("sum", SumPixelFunc,
2851 : pszSumPixelFuncMetadata);
2852 1407 : GDALAddDerivedBandPixelFuncWithArgs("diff", DiffPixelFunc,
2853 : pszDiffPixelFuncMetadata);
2854 1407 : GDALAddDerivedBandPixelFuncWithArgs("mul", MulPixelFunc,
2855 : pszMulPixelFuncMetadata);
2856 1407 : GDALAddDerivedBandPixelFuncWithArgs("div", DivPixelFunc,
2857 : pszDivPixelFuncMetadata);
2858 1407 : GDALAddDerivedBandPixelFunc("cmul", CMulPixelFunc);
2859 1407 : GDALAddDerivedBandPixelFuncWithArgs("inv", InvPixelFunc,
2860 : pszInvPixelFuncMetadata);
2861 1407 : GDALAddDerivedBandPixelFunc("intensity", IntensityPixelFunc);
2862 1407 : GDALAddDerivedBandPixelFuncWithArgs("sqrt", SqrtPixelFunc,
2863 : pszSqrtPixelFuncMetadata);
2864 1407 : GDALAddDerivedBandPixelFuncWithArgs("log10", Log10PixelFunc,
2865 : pszLog10PixelFuncMetadata);
2866 1407 : GDALAddDerivedBandPixelFuncWithArgs("dB", DBPixelFunc,
2867 : pszDBPixelFuncMetadata);
2868 1407 : GDALAddDerivedBandPixelFuncWithArgs("exp", ExpPixelFunc,
2869 : pszExpPixelFuncMetadata);
2870 1407 : GDALAddDerivedBandPixelFunc("dB2amp",
2871 : dB2AmpPixelFunc); // deprecated in v3.5
2872 1407 : GDALAddDerivedBandPixelFunc("dB2pow",
2873 : dB2PowPixelFunc); // deprecated in v3.5
2874 1407 : GDALAddDerivedBandPixelFuncWithArgs("pow", PowPixelFunc,
2875 : pszPowPixelFuncMetadata);
2876 1407 : GDALAddDerivedBandPixelFuncWithArgs("interpolate_linear",
2877 : InterpolatePixelFunc<InterpolateLinear>,
2878 : pszInterpolatePixelFuncMetadata);
2879 1407 : GDALAddDerivedBandPixelFuncWithArgs(
2880 : "interpolate_exp", InterpolatePixelFunc<InterpolateExponential>,
2881 : pszInterpolatePixelFuncMetadata);
2882 1407 : GDALAddDerivedBandPixelFuncWithArgs("replace_nodata",
2883 : ReplaceNoDataPixelFunc,
2884 : pszReplaceNoDataPixelFuncMetadata);
2885 1407 : GDALAddDerivedBandPixelFuncWithArgs("scale", ScalePixelFunc,
2886 : pszScalePixelFuncMetadata);
2887 1407 : GDALAddDerivedBandPixelFuncWithArgs("norm_diff", NormDiffPixelFunc,
2888 : pszNormDiffPixelFuncMetadata);
2889 1407 : GDALAddDerivedBandPixelFuncWithArgs("min", MinPixelFunc,
2890 : pszMinMaxFuncMetadataNodata);
2891 1407 : GDALAddDerivedBandPixelFuncWithArgs("max", MaxPixelFunc,
2892 : pszMinMaxFuncMetadataNodata);
2893 1407 : GDALAddDerivedBandPixelFuncWithArgs("expression", ExprPixelFunc,
2894 : pszExprPixelFuncMetadata);
2895 1407 : GDALAddDerivedBandPixelFuncWithArgs("reclassify", ReclassifyPixelFunc,
2896 : pszReclassifyPixelFuncMetadata);
2897 1407 : GDALAddDerivedBandPixelFuncWithArgs("mean", BasicPixelFunc<MeanKernel>,
2898 : pszBasicPixelFuncMetadata);
2899 1407 : GDALAddDerivedBandPixelFuncWithArgs("geometric_mean",
2900 : BasicPixelFunc<GeoMeanKernel>,
2901 : pszBasicPixelFuncMetadata);
2902 1407 : GDALAddDerivedBandPixelFuncWithArgs("harmonic_mean",
2903 : BasicPixelFunc<HarmonicMeanKernel>,
2904 : pszBasicPixelFuncMetadata);
2905 1407 : GDALAddDerivedBandPixelFuncWithArgs("median", BasicPixelFunc<MedianKernel>,
2906 : pszBasicPixelFuncMetadata);
2907 1407 : GDALAddDerivedBandPixelFuncWithArgs("mode", BasicPixelFunc<ModeKernel>,
2908 : pszBasicPixelFuncMetadata);
2909 1407 : return CE_None;
2910 : }
|
