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 :
19 : #include <limits>
20 :
21 : template <typename T>
22 26057400 : inline double GetSrcVal(const void *pSource, GDALDataType eSrcType, T ii)
23 : {
24 26057400 : switch (eSrcType)
25 : {
26 0 : case GDT_Unknown:
27 0 : return 0;
28 32950 : case GDT_Byte:
29 32950 : return static_cast<const GByte *>(pSource)[ii];
30 0 : case GDT_Int8:
31 0 : return static_cast<const GInt8 *>(pSource)[ii];
32 0 : case GDT_UInt16:
33 0 : return static_cast<const GUInt16 *>(pSource)[ii];
34 0 : case GDT_Int16:
35 0 : return static_cast<const GInt16 *>(pSource)[ii];
36 0 : case GDT_UInt32:
37 0 : return static_cast<const GUInt32 *>(pSource)[ii];
38 400 : case GDT_Int32:
39 400 : return static_cast<const GInt32 *>(pSource)[ii];
40 : // Precision loss currently for int64/uint64
41 0 : case GDT_UInt64:
42 : return static_cast<double>(
43 0 : static_cast<const uint64_t *>(pSource)[ii]);
44 0 : case GDT_Int64:
45 : return static_cast<double>(
46 0 : static_cast<const int64_t *>(pSource)[ii]);
47 0 : case GDT_Float16:
48 0 : return static_cast<const GFloat16 *>(pSource)[ii];
49 8636 : case GDT_Float32:
50 8636 : return static_cast<const float *>(pSource)[ii];
51 25977200 : case GDT_Float64:
52 25977200 : return static_cast<const double *>(pSource)[ii];
53 360 : case GDT_CInt16:
54 360 : return static_cast<const GInt16 *>(pSource)[2 * ii];
55 0 : case GDT_CInt32:
56 0 : return static_cast<const GInt32 *>(pSource)[2 * ii];
57 0 : case GDT_CFloat16:
58 0 : return static_cast<const GFloat16 *>(pSource)[2 * ii];
59 4660 : case GDT_CFloat32:
60 4660 : return static_cast<const float *>(pSource)[2 * ii];
61 33250 : case GDT_CFloat64:
62 33250 : return static_cast<const double *>(pSource)[2 * ii];
63 0 : case GDT_TypeCount:
64 0 : break;
65 : }
66 0 : return 0;
67 : }
68 :
69 57 : static CPLErr FetchDoubleArg(CSLConstList papszArgs, const char *pszName,
70 : double *pdfX, double *pdfDefault = nullptr)
71 : {
72 57 : const char *pszVal = CSLFetchNameValue(papszArgs, pszName);
73 :
74 57 : if (pszVal == nullptr)
75 : {
76 15 : if (pdfDefault == nullptr)
77 : {
78 0 : CPLError(CE_Failure, CPLE_AppDefined,
79 : "Missing pixel function argument: %s", pszName);
80 0 : return CE_Failure;
81 : }
82 : else
83 : {
84 15 : *pdfX = *pdfDefault;
85 15 : return CE_None;
86 : }
87 : }
88 :
89 42 : char *pszEnd = nullptr;
90 42 : *pdfX = std::strtod(pszVal, &pszEnd);
91 42 : if (pszEnd == pszVal)
92 : {
93 0 : CPLError(CE_Failure, CPLE_AppDefined,
94 : "Failed to parse pixel function argument: %s", pszName);
95 0 : return CE_Failure;
96 : }
97 :
98 42 : return CE_None;
99 : }
100 :
101 7 : static CPLErr RealPixelFunc(void **papoSources, int nSources, void *pData,
102 : int nXSize, int nYSize, GDALDataType eSrcType,
103 : GDALDataType eBufType, int nPixelSpace,
104 : int nLineSpace)
105 : {
106 : /* ---- Init ---- */
107 7 : if (nSources != 1)
108 1 : return CE_Failure;
109 :
110 6 : const int nPixelSpaceSrc = GDALGetDataTypeSizeBytes(eSrcType);
111 6 : const size_t nLineSpaceSrc = static_cast<size_t>(nPixelSpaceSrc) * nXSize;
112 :
113 : /* ---- Set pixels ---- */
114 98 : for (int iLine = 0; iLine < nYSize; ++iLine)
115 : {
116 92 : GDALCopyWords(static_cast<GByte *>(papoSources[0]) +
117 92 : nLineSpaceSrc * iLine,
118 : eSrcType, nPixelSpaceSrc,
119 : static_cast<GByte *>(pData) +
120 92 : static_cast<GSpacing>(nLineSpace) * iLine,
121 : eBufType, nPixelSpace, nXSize);
122 : }
123 :
124 : /* ---- Return success ---- */
125 6 : return CE_None;
126 : } // RealPixelFunc
127 :
128 8 : static CPLErr ImagPixelFunc(void **papoSources, int nSources, void *pData,
129 : int nXSize, int nYSize, GDALDataType eSrcType,
130 : GDALDataType eBufType, int nPixelSpace,
131 : int nLineSpace)
132 : {
133 : /* ---- Init ---- */
134 8 : if (nSources != 1)
135 1 : return CE_Failure;
136 :
137 7 : if (GDALDataTypeIsComplex(eSrcType))
138 : {
139 6 : const GDALDataType eSrcBaseType = GDALGetNonComplexDataType(eSrcType);
140 6 : const int nPixelSpaceSrc = GDALGetDataTypeSizeBytes(eSrcType);
141 6 : const size_t nLineSpaceSrc =
142 6 : static_cast<size_t>(nPixelSpaceSrc) * nXSize;
143 :
144 6 : const void *const pImag = static_cast<GByte *>(papoSources[0]) +
145 6 : GDALGetDataTypeSizeBytes(eSrcType) / 2;
146 :
147 : /* ---- Set pixels ---- */
148 56 : for (int iLine = 0; iLine < nYSize; ++iLine)
149 : {
150 50 : GDALCopyWords(static_cast<const GByte *>(pImag) +
151 50 : nLineSpaceSrc * iLine,
152 : eSrcBaseType, nPixelSpaceSrc,
153 : static_cast<GByte *>(pData) +
154 50 : static_cast<GSpacing>(nLineSpace) * iLine,
155 : eBufType, nPixelSpace, nXSize);
156 : }
157 : }
158 : else
159 : {
160 1 : const double dfImag = 0;
161 :
162 : /* ---- Set pixels ---- */
163 21 : for (int iLine = 0; iLine < nYSize; ++iLine)
164 : {
165 : // Always copy from the same location.
166 20 : GDALCopyWords(&dfImag, eSrcType, 0,
167 : static_cast<GByte *>(pData) +
168 20 : static_cast<GSpacing>(nLineSpace) * iLine,
169 : eBufType, nPixelSpace, nXSize);
170 : }
171 : }
172 :
173 : /* ---- Return success ---- */
174 7 : return CE_None;
175 : } // ImagPixelFunc
176 :
177 6 : static CPLErr ComplexPixelFunc(void **papoSources, int nSources, void *pData,
178 : int nXSize, int nYSize, GDALDataType eSrcType,
179 : GDALDataType eBufType, int nPixelSpace,
180 : int nLineSpace)
181 : {
182 : /* ---- Init ---- */
183 6 : if (nSources != 2)
184 1 : return CE_Failure;
185 :
186 5 : const void *const pReal = papoSources[0];
187 5 : const void *const pImag = papoSources[1];
188 :
189 : /* ---- Set pixels ---- */
190 5 : size_t ii = 0;
191 281 : for (int iLine = 0; iLine < nYSize; ++iLine)
192 : {
193 17060 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
194 : {
195 : // Source raster pixels may be obtained with GetSrcVal macro.
196 : const double adfPixVal[2] = {
197 16784 : GetSrcVal(pReal, eSrcType, ii), // re
198 33568 : GetSrcVal(pImag, eSrcType, ii) // im
199 16784 : };
200 :
201 16784 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
202 : static_cast<GByte *>(pData) +
203 16784 : static_cast<GSpacing>(nLineSpace) * iLine +
204 16784 : iCol * nPixelSpace,
205 : eBufType, nPixelSpace, 1);
206 : }
207 : }
208 :
209 : /* ---- Return success ---- */
210 5 : return CE_None;
211 : } // ComplexPixelFunc
212 :
213 : typedef enum
214 : {
215 : GAT_amplitude,
216 : GAT_intensity,
217 : GAT_dB
218 : } PolarAmplitudeType;
219 :
220 : static const char pszPolarPixelFuncMetadata[] =
221 : "<PixelFunctionArgumentsList>"
222 : " <Argument name='amplitude_type' description='Amplitude Type' "
223 : "type='string-select' default='AMPLITUDE'>"
224 : " <Value>INTENSITY</Value>"
225 : " <Value>dB</Value>"
226 : " <Value>AMPLITUDE</Value>"
227 : " </Argument>"
228 : "</PixelFunctionArgumentsList>";
229 :
230 4 : static CPLErr PolarPixelFunc(void **papoSources, int nSources, void *pData,
231 : int nXSize, int nYSize, GDALDataType eSrcType,
232 : GDALDataType eBufType, int nPixelSpace,
233 : int nLineSpace, CSLConstList papszArgs)
234 : {
235 : /* ---- Init ---- */
236 4 : if (nSources != 2)
237 0 : return CE_Failure;
238 :
239 4 : const char pszName[] = "amplitude_type";
240 4 : const char *pszVal = CSLFetchNameValue(papszArgs, pszName);
241 4 : PolarAmplitudeType amplitudeType = GAT_amplitude;
242 4 : if (pszVal != nullptr)
243 : {
244 3 : if (strcmp(pszVal, "INTENSITY") == 0)
245 1 : amplitudeType = GAT_intensity;
246 2 : else if (strcmp(pszVal, "dB") == 0)
247 1 : amplitudeType = GAT_dB;
248 1 : else if (strcmp(pszVal, "AMPLITUDE") != 0)
249 : {
250 0 : CPLError(CE_Failure, CPLE_AppDefined,
251 : "Invalid value for pixel function argument '%s': %s",
252 : pszName, pszVal);
253 0 : return CE_Failure;
254 : }
255 : }
256 :
257 4 : const void *const pAmp = papoSources[0];
258 4 : const void *const pPhase = papoSources[1];
259 :
260 : /* ---- Set pixels ---- */
261 4 : size_t ii = 0;
262 84 : for (int iLine = 0; iLine < nYSize; ++iLine)
263 : {
264 1680 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
265 : {
266 : // Source raster pixels may be obtained with GetSrcVal macro.
267 1600 : double dfAmp = GetSrcVal(pAmp, eSrcType, ii);
268 1600 : switch (amplitudeType)
269 : {
270 400 : case GAT_intensity:
271 : // clip to zero
272 400 : dfAmp = dfAmp <= 0 ? 0 : std::sqrt(dfAmp);
273 400 : break;
274 400 : case GAT_dB:
275 400 : dfAmp = dfAmp <= 0
276 400 : ? -std::numeric_limits<double>::infinity()
277 400 : : pow(10, dfAmp / 20.);
278 400 : break;
279 800 : case GAT_amplitude:
280 800 : break;
281 : }
282 1600 : const double dfPhase = GetSrcVal(pPhase, eSrcType, ii);
283 : const double adfPixVal[2] = {
284 1600 : dfAmp * std::cos(dfPhase), // re
285 1600 : dfAmp * std::sin(dfPhase) // im
286 1600 : };
287 :
288 1600 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
289 : static_cast<GByte *>(pData) +
290 1600 : static_cast<GSpacing>(nLineSpace) * iLine +
291 1600 : iCol * nPixelSpace,
292 : eBufType, nPixelSpace, 1);
293 : }
294 : }
295 :
296 : /* ---- Return success ---- */
297 4 : return CE_None;
298 : } // PolarPixelFunc
299 :
300 4 : static CPLErr ModulePixelFunc(void **papoSources, int nSources, void *pData,
301 : int nXSize, int nYSize, GDALDataType eSrcType,
302 : GDALDataType eBufType, int nPixelSpace,
303 : int nLineSpace)
304 : {
305 : /* ---- Init ---- */
306 4 : if (nSources != 1)
307 1 : return CE_Failure;
308 :
309 3 : if (GDALDataTypeIsComplex(eSrcType))
310 : {
311 2 : const void *pReal = papoSources[0];
312 2 : const void *pImag = static_cast<GByte *>(papoSources[0]) +
313 2 : GDALGetDataTypeSizeBytes(eSrcType) / 2;
314 :
315 : /* ---- Set pixels ---- */
316 2 : size_t ii = 0;
317 14 : for (int iLine = 0; iLine < nYSize; ++iLine)
318 : {
319 72 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
320 : {
321 : // Source raster pixels may be obtained with GetSrcVal macro.
322 60 : const double dfReal = GetSrcVal(pReal, eSrcType, ii);
323 60 : const double dfImag = GetSrcVal(pImag, eSrcType, ii);
324 :
325 60 : const double dfPixVal = sqrt(dfReal * dfReal + dfImag * dfImag);
326 :
327 60 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
328 : static_cast<GByte *>(pData) +
329 60 : static_cast<GSpacing>(nLineSpace) * iLine +
330 60 : iCol * nPixelSpace,
331 : eBufType, nPixelSpace, 1);
332 : }
333 : }
334 : }
335 : else
336 : {
337 : /* ---- Set pixels ---- */
338 1 : size_t ii = 0;
339 21 : for (int iLine = 0; iLine < nYSize; ++iLine)
340 : {
341 420 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
342 : {
343 : // Source raster pixels may be obtained with GetSrcVal macro.
344 : const double dfPixVal =
345 400 : fabs(GetSrcVal(papoSources[0], eSrcType, ii));
346 :
347 400 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
348 : static_cast<GByte *>(pData) +
349 400 : static_cast<GSpacing>(nLineSpace) * iLine +
350 400 : iCol * nPixelSpace,
351 : eBufType, nPixelSpace, 1);
352 : }
353 : }
354 : }
355 :
356 : /* ---- Return success ---- */
357 3 : return CE_None;
358 : } // ModulePixelFunc
359 :
360 5 : static CPLErr PhasePixelFunc(void **papoSources, int nSources, void *pData,
361 : int nXSize, int nYSize, GDALDataType eSrcType,
362 : GDALDataType eBufType, int nPixelSpace,
363 : int nLineSpace)
364 : {
365 : /* ---- Init ---- */
366 5 : if (nSources != 1)
367 1 : return CE_Failure;
368 :
369 4 : if (GDALDataTypeIsComplex(eSrcType))
370 : {
371 2 : const void *const pReal = papoSources[0];
372 2 : const void *const pImag = static_cast<GByte *>(papoSources[0]) +
373 2 : GDALGetDataTypeSizeBytes(eSrcType) / 2;
374 :
375 : /* ---- Set pixels ---- */
376 2 : size_t ii = 0;
377 14 : for (int iLine = 0; iLine < nYSize; ++iLine)
378 : {
379 72 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
380 : {
381 : // Source raster pixels may be obtained with GetSrcVal macro.
382 60 : const double dfReal = GetSrcVal(pReal, eSrcType, ii);
383 60 : const double dfImag = GetSrcVal(pImag, eSrcType, ii);
384 :
385 60 : const double dfPixVal = atan2(dfImag, dfReal);
386 :
387 60 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
388 : static_cast<GByte *>(pData) +
389 60 : static_cast<GSpacing>(nLineSpace) * iLine +
390 60 : iCol * nPixelSpace,
391 : eBufType, nPixelSpace, 1);
392 : }
393 : }
394 : }
395 2 : else if (GDALDataTypeIsInteger(eSrcType) && !GDALDataTypeIsSigned(eSrcType))
396 : {
397 1 : constexpr double dfZero = 0;
398 7 : for (int iLine = 0; iLine < nYSize; ++iLine)
399 : {
400 6 : GDALCopyWords(&dfZero, GDT_Float64, 0,
401 : static_cast<GByte *>(pData) +
402 6 : static_cast<GSpacing>(nLineSpace) * iLine,
403 : eBufType, nPixelSpace, nXSize);
404 : }
405 : }
406 : else
407 : {
408 : /* ---- Set pixels ---- */
409 1 : size_t ii = 0;
410 7 : for (int iLine = 0; iLine < nYSize; ++iLine)
411 : {
412 36 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
413 : {
414 30 : const void *const pReal = papoSources[0];
415 :
416 : // Source raster pixels may be obtained with GetSrcVal macro.
417 30 : const double dfReal = GetSrcVal(pReal, eSrcType, ii);
418 30 : const double dfPixVal = (dfReal < 0) ? M_PI : 0.0;
419 :
420 30 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
421 : static_cast<GByte *>(pData) +
422 30 : static_cast<GSpacing>(nLineSpace) * iLine +
423 30 : iCol * nPixelSpace,
424 : eBufType, nPixelSpace, 1);
425 : }
426 : }
427 : }
428 :
429 : /* ---- Return success ---- */
430 4 : return CE_None;
431 : } // PhasePixelFunc
432 :
433 4 : static CPLErr ConjPixelFunc(void **papoSources, int nSources, void *pData,
434 : int nXSize, int nYSize, GDALDataType eSrcType,
435 : GDALDataType eBufType, int nPixelSpace,
436 : int nLineSpace)
437 : {
438 : /* ---- Init ---- */
439 4 : if (nSources != 1)
440 1 : return CE_Failure;
441 :
442 3 : if (GDALDataTypeIsComplex(eSrcType) && GDALDataTypeIsComplex(eBufType))
443 : {
444 2 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
445 2 : const void *const pReal = papoSources[0];
446 2 : const void *const pImag =
447 2 : static_cast<GByte *>(papoSources[0]) + nOffset;
448 :
449 : /* ---- Set pixels ---- */
450 2 : size_t ii = 0;
451 14 : for (int iLine = 0; iLine < nYSize; ++iLine)
452 : {
453 72 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
454 : {
455 : // Source raster pixels may be obtained with GetSrcVal macro.
456 : const double adfPixVal[2] = {
457 60 : +GetSrcVal(pReal, eSrcType, ii), // re
458 120 : -GetSrcVal(pImag, eSrcType, ii) // im
459 60 : };
460 :
461 60 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
462 : static_cast<GByte *>(pData) +
463 60 : static_cast<GSpacing>(nLineSpace) * iLine +
464 60 : iCol * nPixelSpace,
465 : eBufType, nPixelSpace, 1);
466 : }
467 : }
468 : }
469 : else
470 : {
471 : // No complex data type.
472 1 : return RealPixelFunc(papoSources, nSources, pData, nXSize, nYSize,
473 1 : eSrcType, eBufType, nPixelSpace, nLineSpace);
474 : }
475 :
476 : /* ---- Return success ---- */
477 2 : return CE_None;
478 : } // ConjPixelFunc
479 :
480 : static const char pszSumPixelFuncMetadata[] =
481 : "<PixelFunctionArgumentsList>"
482 : " <Argument name='k' description='Optional constant term' type='double' "
483 : "default='0.0' />"
484 : "</PixelFunctionArgumentsList>";
485 :
486 4 : static CPLErr SumPixelFunc(void **papoSources, int nSources, void *pData,
487 : int nXSize, int nYSize, GDALDataType eSrcType,
488 : GDALDataType eBufType, int nPixelSpace,
489 : int nLineSpace, CSLConstList papszArgs)
490 : {
491 : /* ---- Init ---- */
492 4 : if (nSources < 2)
493 1 : return CE_Failure;
494 :
495 3 : double dfK = 0.0;
496 3 : if (FetchDoubleArg(papszArgs, "k", &dfK, &dfK) != CE_None)
497 0 : return CE_Failure;
498 :
499 : /* ---- Set pixels ---- */
500 3 : if (GDALDataTypeIsComplex(eSrcType))
501 : {
502 1 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
503 :
504 : /* ---- Set pixels ---- */
505 1 : size_t ii = 0;
506 7 : for (int iLine = 0; iLine < nYSize; ++iLine)
507 : {
508 36 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
509 : {
510 30 : double adfSum[2] = {dfK, 0.0};
511 :
512 120 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
513 : {
514 90 : const void *const pReal = papoSources[iSrc];
515 90 : const void *const pImag =
516 90 : static_cast<const GByte *>(pReal) + nOffset;
517 :
518 : // Source raster pixels may be obtained with GetSrcVal
519 : // macro.
520 90 : adfSum[0] += GetSrcVal(pReal, eSrcType, ii);
521 90 : adfSum[1] += GetSrcVal(pImag, eSrcType, ii);
522 : }
523 :
524 30 : GDALCopyWords(adfSum, GDT_CFloat64, 0,
525 : static_cast<GByte *>(pData) +
526 30 : static_cast<GSpacing>(nLineSpace) * iLine +
527 30 : iCol * nPixelSpace,
528 : eBufType, nPixelSpace, 1);
529 : }
530 : }
531 : }
532 : else
533 : {
534 : /* ---- Set pixels ---- */
535 2 : size_t ii = 0;
536 42 : for (int iLine = 0; iLine < nYSize; ++iLine)
537 : {
538 840 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
539 : {
540 800 : double dfSum = dfK; // Not complex.
541 :
542 3200 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
543 : {
544 : // Source raster pixels may be obtained with GetSrcVal
545 : // macro.
546 2400 : dfSum += GetSrcVal(papoSources[iSrc], eSrcType, ii);
547 : }
548 :
549 800 : GDALCopyWords(&dfSum, GDT_Float64, 0,
550 : static_cast<GByte *>(pData) +
551 800 : static_cast<GSpacing>(nLineSpace) * iLine +
552 800 : iCol * nPixelSpace,
553 : eBufType, nPixelSpace, 1);
554 : }
555 : }
556 : }
557 :
558 : /* ---- Return success ---- */
559 3 : return CE_None;
560 : } /* SumPixelFunc */
561 :
562 3 : static CPLErr DiffPixelFunc(void **papoSources, int nSources, void *pData,
563 : int nXSize, int nYSize, GDALDataType eSrcType,
564 : GDALDataType eBufType, int nPixelSpace,
565 : int nLineSpace)
566 : {
567 : /* ---- Init ---- */
568 3 : if (nSources != 2)
569 1 : return CE_Failure;
570 :
571 2 : if (GDALDataTypeIsComplex(eSrcType))
572 : {
573 1 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
574 1 : const void *const pReal0 = papoSources[0];
575 1 : const void *const pImag0 =
576 1 : static_cast<GByte *>(papoSources[0]) + nOffset;
577 1 : const void *const pReal1 = papoSources[1];
578 1 : const void *const pImag1 =
579 1 : static_cast<GByte *>(papoSources[1]) + nOffset;
580 :
581 : /* ---- Set pixels ---- */
582 1 : size_t ii = 0;
583 7 : for (int iLine = 0; iLine < nYSize; ++iLine)
584 : {
585 36 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
586 : {
587 : // Source raster pixels may be obtained with GetSrcVal macro.
588 30 : double adfPixVal[2] = {GetSrcVal(pReal0, eSrcType, ii) -
589 30 : GetSrcVal(pReal1, eSrcType, ii),
590 90 : GetSrcVal(pImag0, eSrcType, ii) -
591 30 : GetSrcVal(pImag1, eSrcType, ii)};
592 :
593 30 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
594 : static_cast<GByte *>(pData) +
595 30 : static_cast<GSpacing>(nLineSpace) * iLine +
596 30 : iCol * nPixelSpace,
597 : eBufType, nPixelSpace, 1);
598 : }
599 : }
600 : }
601 : else
602 : {
603 : /* ---- Set pixels ---- */
604 1 : size_t ii = 0;
605 7 : for (int iLine = 0; iLine < nYSize; ++iLine)
606 : {
607 36 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
608 : {
609 : // Source raster pixels may be obtained with GetSrcVal macro.
610 : // Not complex.
611 : const double dfPixVal =
612 30 : GetSrcVal(papoSources[0], eSrcType, ii) -
613 30 : GetSrcVal(papoSources[1], eSrcType, ii);
614 :
615 30 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
616 : static_cast<GByte *>(pData) +
617 30 : static_cast<GSpacing>(nLineSpace) * iLine +
618 30 : iCol * nPixelSpace,
619 : eBufType, nPixelSpace, 1);
620 : }
621 : }
622 : }
623 :
624 : /* ---- Return success ---- */
625 2 : return CE_None;
626 : } // DiffPixelFunc
627 :
628 : static const char pszMulPixelFuncMetadata[] =
629 : "<PixelFunctionArgumentsList>"
630 : " <Argument name='k' description='Optional constant factor' "
631 : "type='double' default='1.0' />"
632 : "</PixelFunctionArgumentsList>";
633 :
634 4 : static CPLErr MulPixelFunc(void **papoSources, int nSources, void *pData,
635 : int nXSize, int nYSize, GDALDataType eSrcType,
636 : GDALDataType eBufType, int nPixelSpace,
637 : int nLineSpace, CSLConstList papszArgs)
638 : {
639 : /* ---- Init ---- */
640 4 : if (nSources < 2)
641 1 : return CE_Failure;
642 :
643 3 : double dfK = 1.0;
644 3 : if (FetchDoubleArg(papszArgs, "k", &dfK, &dfK) != CE_None)
645 0 : return CE_Failure;
646 :
647 : /* ---- Set pixels ---- */
648 3 : if (GDALDataTypeIsComplex(eSrcType))
649 : {
650 1 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
651 :
652 : /* ---- Set pixels ---- */
653 1 : size_t ii = 0;
654 7 : for (int iLine = 0; iLine < nYSize; ++iLine)
655 : {
656 36 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
657 : {
658 30 : double adfPixVal[2] = {dfK, 0.0};
659 :
660 90 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
661 : {
662 60 : const void *const pReal = papoSources[iSrc];
663 60 : const void *const pImag =
664 60 : static_cast<const GByte *>(pReal) + nOffset;
665 :
666 60 : const double dfOldR = adfPixVal[0];
667 60 : const double dfOldI = adfPixVal[1];
668 :
669 : // Source raster pixels may be obtained with GetSrcVal
670 : // macro.
671 60 : const double dfNewR = GetSrcVal(pReal, eSrcType, ii);
672 60 : const double dfNewI = GetSrcVal(pImag, eSrcType, ii);
673 :
674 60 : adfPixVal[0] = dfOldR * dfNewR - dfOldI * dfNewI;
675 60 : adfPixVal[1] = dfOldR * dfNewI + dfOldI * dfNewR;
676 : }
677 :
678 30 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
679 : static_cast<GByte *>(pData) +
680 30 : static_cast<GSpacing>(nLineSpace) * iLine +
681 30 : iCol * nPixelSpace,
682 : eBufType, nPixelSpace, 1);
683 : }
684 : }
685 : }
686 : else
687 : {
688 : /* ---- Set pixels ---- */
689 2 : size_t ii = 0;
690 42 : for (int iLine = 0; iLine < nYSize; ++iLine)
691 : {
692 840 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
693 : {
694 800 : double dfPixVal = dfK; // Not complex.
695 :
696 3200 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
697 : {
698 : // Source raster pixels may be obtained with GetSrcVal
699 : // macro.
700 2400 : dfPixVal *= GetSrcVal(papoSources[iSrc], eSrcType, ii);
701 : }
702 :
703 800 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
704 : static_cast<GByte *>(pData) +
705 800 : static_cast<GSpacing>(nLineSpace) * iLine +
706 800 : iCol * nPixelSpace,
707 : eBufType, nPixelSpace, 1);
708 : }
709 : }
710 : }
711 :
712 : /* ---- Return success ---- */
713 3 : return CE_None;
714 : } // MulPixelFunc
715 :
716 2 : static CPLErr DivPixelFunc(void **papoSources, int nSources, void *pData,
717 : int nXSize, int nYSize, GDALDataType eSrcType,
718 : GDALDataType eBufType, int nPixelSpace,
719 : int nLineSpace)
720 : {
721 : /* ---- Init ---- */
722 2 : if (nSources != 2)
723 0 : return CE_Failure;
724 :
725 : /* ---- Set pixels ---- */
726 2 : if (GDALDataTypeIsComplex(eSrcType))
727 : {
728 1 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
729 1 : const void *const pReal0 = papoSources[0];
730 1 : const void *const pImag0 =
731 1 : static_cast<GByte *>(papoSources[0]) + nOffset;
732 1 : const void *const pReal1 = papoSources[1];
733 1 : const void *const pImag1 =
734 1 : static_cast<GByte *>(papoSources[1]) + nOffset;
735 :
736 : /* ---- Set pixels ---- */
737 1 : size_t ii = 0;
738 7 : for (int iLine = 0; iLine < nYSize; ++iLine)
739 : {
740 36 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
741 : {
742 : // Source raster pixels may be obtained with GetSrcVal macro.
743 30 : const double dfReal0 = GetSrcVal(pReal0, eSrcType, ii);
744 30 : const double dfReal1 = GetSrcVal(pReal1, eSrcType, ii);
745 30 : const double dfImag0 = GetSrcVal(pImag0, eSrcType, ii);
746 30 : const double dfImag1 = GetSrcVal(pImag1, eSrcType, ii);
747 30 : const double dfAux = dfReal1 * dfReal1 + dfImag1 * dfImag1;
748 :
749 : const double adfPixVal[2] = {
750 : dfAux == 0
751 30 : ? std::numeric_limits<double>::infinity()
752 30 : : dfReal0 * dfReal1 / dfAux + dfImag0 * dfImag1 / dfAux,
753 0 : dfAux == 0 ? std::numeric_limits<double>::infinity()
754 30 : : dfReal1 / dfAux * dfImag0 -
755 30 : dfReal0 * dfImag1 / dfAux};
756 :
757 30 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
758 : static_cast<GByte *>(pData) +
759 30 : static_cast<GSpacing>(nLineSpace) * iLine +
760 30 : iCol * nPixelSpace,
761 : eBufType, nPixelSpace, 1);
762 : }
763 : }
764 : }
765 : else
766 : {
767 : /* ---- Set pixels ---- */
768 1 : size_t ii = 0;
769 7 : for (int iLine = 0; iLine < nYSize; ++iLine)
770 : {
771 36 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
772 : {
773 30 : const double dfVal = GetSrcVal(papoSources[1], eSrcType, ii);
774 : // coverity[divide_by_zero]
775 : double dfPixVal =
776 : dfVal == 0
777 30 : ? std::numeric_limits<double>::infinity()
778 30 : : GetSrcVal(papoSources[0], eSrcType, ii) / dfVal;
779 :
780 30 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
781 : static_cast<GByte *>(pData) +
782 30 : static_cast<GSpacing>(nLineSpace) * iLine +
783 30 : iCol * nPixelSpace,
784 : eBufType, nPixelSpace, 1);
785 : }
786 : }
787 : }
788 :
789 : /* ---- Return success ---- */
790 2 : return CE_None;
791 : } // DivPixelFunc
792 :
793 3 : static CPLErr CMulPixelFunc(void **papoSources, int nSources, void *pData,
794 : int nXSize, int nYSize, GDALDataType eSrcType,
795 : GDALDataType eBufType, int nPixelSpace,
796 : int nLineSpace)
797 : {
798 : /* ---- Init ---- */
799 3 : if (nSources != 2)
800 1 : return CE_Failure;
801 :
802 : /* ---- Set pixels ---- */
803 2 : if (GDALDataTypeIsComplex(eSrcType))
804 : {
805 1 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
806 1 : const void *const pReal0 = papoSources[0];
807 1 : const void *const pImag0 =
808 1 : static_cast<GByte *>(papoSources[0]) + nOffset;
809 1 : const void *const pReal1 = papoSources[1];
810 1 : const void *const pImag1 =
811 1 : static_cast<GByte *>(papoSources[1]) + nOffset;
812 :
813 1 : size_t ii = 0;
814 7 : for (int iLine = 0; iLine < nYSize; ++iLine)
815 : {
816 36 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
817 : {
818 : // Source raster pixels may be obtained with GetSrcVal macro.
819 30 : const double dfReal0 = GetSrcVal(pReal0, eSrcType, ii);
820 30 : const double dfReal1 = GetSrcVal(pReal1, eSrcType, ii);
821 30 : const double dfImag0 = GetSrcVal(pImag0, eSrcType, ii);
822 30 : const double dfImag1 = GetSrcVal(pImag1, eSrcType, ii);
823 : const double adfPixVal[2] = {
824 30 : dfReal0 * dfReal1 + dfImag0 * dfImag1,
825 30 : dfReal1 * dfImag0 - dfReal0 * dfImag1};
826 :
827 30 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
828 : static_cast<GByte *>(pData) +
829 30 : static_cast<GSpacing>(nLineSpace) * iLine +
830 30 : iCol * nPixelSpace,
831 : eBufType, nPixelSpace, 1);
832 : }
833 : }
834 : }
835 : else
836 : {
837 1 : size_t ii = 0;
838 21 : for (int iLine = 0; iLine < nYSize; ++iLine)
839 : {
840 420 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
841 : {
842 : // Source raster pixels may be obtained with GetSrcVal macro.
843 : // Not complex.
844 400 : const double adfPixVal[2] = {
845 400 : GetSrcVal(papoSources[0], eSrcType, ii) *
846 400 : GetSrcVal(papoSources[1], eSrcType, ii),
847 400 : 0.0};
848 :
849 400 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
850 : static_cast<GByte *>(pData) +
851 400 : static_cast<GSpacing>(nLineSpace) * iLine +
852 400 : iCol * nPixelSpace,
853 : eBufType, nPixelSpace, 1);
854 : }
855 : }
856 : }
857 :
858 : /* ---- Return success ---- */
859 2 : return CE_None;
860 : } // CMulPixelFunc
861 :
862 : static const char pszInvPixelFuncMetadata[] =
863 : "<PixelFunctionArgumentsList>"
864 : " <Argument name='k' description='Optional constant factor' "
865 : "type='double' default='1.0' />"
866 : "</PixelFunctionArgumentsList>";
867 :
868 6 : static CPLErr InvPixelFunc(void **papoSources, int nSources, void *pData,
869 : int nXSize, int nYSize, GDALDataType eSrcType,
870 : GDALDataType eBufType, int nPixelSpace,
871 : int nLineSpace, CSLConstList papszArgs)
872 : {
873 : /* ---- Init ---- */
874 6 : if (nSources != 1)
875 1 : return CE_Failure;
876 :
877 5 : double dfK = 1.0;
878 5 : if (FetchDoubleArg(papszArgs, "k", &dfK, &dfK) != CE_None)
879 0 : return CE_Failure;
880 :
881 : /* ---- Set pixels ---- */
882 5 : if (GDALDataTypeIsComplex(eSrcType))
883 : {
884 2 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
885 2 : const void *const pReal = papoSources[0];
886 2 : const void *const pImag =
887 2 : static_cast<GByte *>(papoSources[0]) + nOffset;
888 :
889 2 : size_t ii = 0;
890 9 : for (int iLine = 0; iLine < nYSize; ++iLine)
891 : {
892 38 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
893 : {
894 : // Source raster pixels may be obtained with GetSrcVal macro.
895 31 : const double dfReal = GetSrcVal(pReal, eSrcType, ii);
896 31 : const double dfImag = GetSrcVal(pImag, eSrcType, ii);
897 31 : const double dfAux = dfReal * dfReal + dfImag * dfImag;
898 : const double adfPixVal[2] = {
899 31 : dfAux == 0 ? std::numeric_limits<double>::infinity()
900 30 : : dfK * dfReal / dfAux,
901 1 : dfAux == 0 ? std::numeric_limits<double>::infinity()
902 31 : : -dfK * dfImag / dfAux};
903 :
904 31 : GDALCopyWords(adfPixVal, GDT_CFloat64, 0,
905 : static_cast<GByte *>(pData) +
906 31 : static_cast<GSpacing>(nLineSpace) * iLine +
907 31 : iCol * nPixelSpace,
908 : eBufType, nPixelSpace, 1);
909 : }
910 : }
911 : }
912 : else
913 : {
914 : /* ---- Set pixels ---- */
915 3 : size_t ii = 0;
916 44 : for (int iLine = 0; iLine < nYSize; ++iLine)
917 : {
918 842 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
919 : {
920 : // Source raster pixels may be obtained with GetSrcVal macro.
921 : // Not complex.
922 801 : const double dfVal = GetSrcVal(papoSources[0], eSrcType, ii);
923 : // coverity[divide_by_zero]
924 : const double dfPixVal =
925 801 : dfVal == 0 ? std::numeric_limits<double>::infinity()
926 800 : : dfK / dfVal;
927 :
928 801 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
929 : static_cast<GByte *>(pData) +
930 801 : static_cast<GSpacing>(nLineSpace) * iLine +
931 801 : iCol * nPixelSpace,
932 : eBufType, nPixelSpace, 1);
933 : }
934 : }
935 : }
936 :
937 : /* ---- Return success ---- */
938 5 : return CE_None;
939 : } // InvPixelFunc
940 :
941 4 : static CPLErr IntensityPixelFunc(void **papoSources, int nSources, void *pData,
942 : int nXSize, int nYSize, GDALDataType eSrcType,
943 : GDALDataType eBufType, int nPixelSpace,
944 : int nLineSpace)
945 : {
946 : /* ---- Init ---- */
947 4 : if (nSources != 1)
948 1 : return CE_Failure;
949 :
950 3 : if (GDALDataTypeIsComplex(eSrcType))
951 : {
952 2 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
953 2 : const void *const pReal = papoSources[0];
954 2 : const void *const pImag =
955 2 : static_cast<GByte *>(papoSources[0]) + nOffset;
956 :
957 : /* ---- Set pixels ---- */
958 2 : size_t ii = 0;
959 14 : for (int iLine = 0; iLine < nYSize; ++iLine)
960 : {
961 72 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
962 : {
963 : // Source raster pixels may be obtained with GetSrcVal macro.
964 60 : const double dfReal = GetSrcVal(pReal, eSrcType, ii);
965 60 : const double dfImag = GetSrcVal(pImag, eSrcType, ii);
966 :
967 60 : const double dfPixVal = dfReal * dfReal + dfImag * dfImag;
968 :
969 60 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
970 : static_cast<GByte *>(pData) +
971 60 : static_cast<GSpacing>(nLineSpace) * iLine +
972 60 : iCol * nPixelSpace,
973 : eBufType, nPixelSpace, 1);
974 : }
975 : }
976 : }
977 : else
978 : {
979 : /* ---- Set pixels ---- */
980 1 : size_t ii = 0;
981 21 : for (int iLine = 0; iLine < nYSize; ++iLine)
982 : {
983 420 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
984 : {
985 : // Source raster pixels may be obtained with GetSrcVal macro.
986 400 : double dfPixVal = GetSrcVal(papoSources[0], eSrcType, ii);
987 400 : dfPixVal *= dfPixVal;
988 :
989 400 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
990 : static_cast<GByte *>(pData) +
991 400 : static_cast<GSpacing>(nLineSpace) * iLine +
992 400 : iCol * nPixelSpace,
993 : eBufType, nPixelSpace, 1);
994 : }
995 : }
996 : }
997 :
998 : /* ---- Return success ---- */
999 3 : return CE_None;
1000 : } // IntensityPixelFunc
1001 :
1002 2 : static CPLErr SqrtPixelFunc(void **papoSources, int nSources, void *pData,
1003 : int nXSize, int nYSize, GDALDataType eSrcType,
1004 : GDALDataType eBufType, int nPixelSpace,
1005 : int nLineSpace)
1006 : {
1007 : /* ---- Init ---- */
1008 2 : if (nSources != 1)
1009 1 : return CE_Failure;
1010 1 : if (GDALDataTypeIsComplex(eSrcType))
1011 0 : return CE_Failure;
1012 :
1013 : /* ---- Set pixels ---- */
1014 1 : size_t ii = 0;
1015 21 : for (int iLine = 0; iLine < nYSize; ++iLine)
1016 : {
1017 420 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1018 : {
1019 : // Source raster pixels may be obtained with GetSrcVal macro.
1020 : const double dfPixVal =
1021 400 : sqrt(GetSrcVal(papoSources[0], eSrcType, ii));
1022 :
1023 400 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1024 : static_cast<GByte *>(pData) +
1025 400 : static_cast<GSpacing>(nLineSpace) * iLine +
1026 400 : iCol * nPixelSpace,
1027 : eBufType, nPixelSpace, 1);
1028 : }
1029 : }
1030 :
1031 : /* ---- Return success ---- */
1032 1 : return CE_None;
1033 : } // SqrtPixelFunc
1034 :
1035 11 : static CPLErr Log10PixelFuncHelper(void **papoSources, int nSources,
1036 : void *pData, int nXSize, int nYSize,
1037 : GDALDataType eSrcType, GDALDataType eBufType,
1038 : int nPixelSpace, int nLineSpace, double fact)
1039 : {
1040 : /* ---- Init ---- */
1041 11 : if (nSources != 1)
1042 2 : return CE_Failure;
1043 :
1044 9 : if (GDALDataTypeIsComplex(eSrcType))
1045 : {
1046 : // Complex input datatype.
1047 5 : const int nOffset = GDALGetDataTypeSizeBytes(eSrcType) / 2;
1048 5 : const void *const pReal = papoSources[0];
1049 5 : const void *const pImag =
1050 5 : static_cast<GByte *>(papoSources[0]) + nOffset;
1051 :
1052 : /* We should compute fact * log10( sqrt( dfReal * dfReal + dfImag *
1053 : * dfImag ) ) */
1054 : /* Given that log10(sqrt(x)) = 0.5 * log10(x) */
1055 : /* we can remove the sqrt() by multiplying fact by 0.5 */
1056 5 : fact *= 0.5;
1057 :
1058 : /* ---- Set pixels ---- */
1059 5 : size_t ii = 0;
1060 35 : for (int iLine = 0; iLine < nYSize; ++iLine)
1061 : {
1062 180 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1063 : {
1064 : // Source raster pixels may be obtained with GetSrcVal macro.
1065 150 : const double dfReal = GetSrcVal(pReal, eSrcType, ii);
1066 150 : const double dfImag = GetSrcVal(pImag, eSrcType, ii);
1067 :
1068 : const double dfPixVal =
1069 150 : fact * log10(dfReal * dfReal + dfImag * dfImag);
1070 :
1071 150 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1072 : static_cast<GByte *>(pData) +
1073 150 : static_cast<GSpacing>(nLineSpace) * iLine +
1074 150 : iCol * nPixelSpace,
1075 : eBufType, nPixelSpace, 1);
1076 : }
1077 : }
1078 : }
1079 : else
1080 : {
1081 : /* ---- Set pixels ---- */
1082 4 : size_t ii = 0;
1083 84 : for (int iLine = 0; iLine < nYSize; ++iLine)
1084 : {
1085 1680 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1086 : {
1087 : // Source raster pixels may be obtained with GetSrcVal macro.
1088 : const double dfPixVal =
1089 1600 : fact * log10(fabs(GetSrcVal(papoSources[0], eSrcType, ii)));
1090 :
1091 1600 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1092 : static_cast<GByte *>(pData) +
1093 1600 : static_cast<GSpacing>(nLineSpace) * iLine +
1094 1600 : iCol * nPixelSpace,
1095 : eBufType, nPixelSpace, 1);
1096 : }
1097 : }
1098 : }
1099 :
1100 : /* ---- Return success ---- */
1101 9 : return CE_None;
1102 : } // Log10PixelFuncHelper
1103 :
1104 4 : static CPLErr Log10PixelFunc(void **papoSources, int nSources, void *pData,
1105 : int nXSize, int nYSize, GDALDataType eSrcType,
1106 : GDALDataType eBufType, int nPixelSpace,
1107 : int nLineSpace)
1108 : {
1109 4 : return Log10PixelFuncHelper(papoSources, nSources, pData, nXSize, nYSize,
1110 : eSrcType, eBufType, nPixelSpace, nLineSpace,
1111 4 : 1.0);
1112 : } // Log10PixelFunc
1113 :
1114 : static const char pszDBPixelFuncMetadata[] =
1115 : "<PixelFunctionArgumentsList>"
1116 : " <Argument name='fact' description='Factor' type='double' "
1117 : "default='20.0' />"
1118 : "</PixelFunctionArgumentsList>";
1119 :
1120 7 : static CPLErr DBPixelFunc(void **papoSources, int nSources, void *pData,
1121 : int nXSize, int nYSize, GDALDataType eSrcType,
1122 : GDALDataType eBufType, int nPixelSpace,
1123 : int nLineSpace, CSLConstList papszArgs)
1124 : {
1125 7 : double dfFact = 20.;
1126 7 : if (FetchDoubleArg(papszArgs, "fact", &dfFact, &dfFact) != CE_None)
1127 0 : return CE_Failure;
1128 :
1129 7 : return Log10PixelFuncHelper(papoSources, nSources, pData, nXSize, nYSize,
1130 : eSrcType, eBufType, nPixelSpace, nLineSpace,
1131 7 : dfFact);
1132 : } // DBPixelFunc
1133 :
1134 7 : static CPLErr ExpPixelFuncHelper(void **papoSources, int nSources, void *pData,
1135 : int nXSize, int nYSize, GDALDataType eSrcType,
1136 : GDALDataType eBufType, int nPixelSpace,
1137 : int nLineSpace, double base, double fact)
1138 : {
1139 : /* ---- Init ---- */
1140 7 : if (nSources != 1)
1141 2 : return CE_Failure;
1142 5 : if (GDALDataTypeIsComplex(eSrcType))
1143 0 : return CE_Failure;
1144 :
1145 : /* ---- Set pixels ---- */
1146 5 : size_t ii = 0;
1147 105 : for (int iLine = 0; iLine < nYSize; ++iLine)
1148 : {
1149 2100 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1150 : {
1151 : // Source raster pixels may be obtained with GetSrcVal macro.
1152 : const double dfPixVal =
1153 2000 : pow(base, GetSrcVal(papoSources[0], eSrcType, ii) * fact);
1154 :
1155 2000 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1156 : static_cast<GByte *>(pData) +
1157 2000 : static_cast<GSpacing>(nLineSpace) * iLine +
1158 2000 : iCol * nPixelSpace,
1159 : eBufType, nPixelSpace, 1);
1160 : }
1161 : }
1162 :
1163 : /* ---- Return success ---- */
1164 5 : return CE_None;
1165 : } // ExpPixelFuncHelper
1166 :
1167 : static const char pszExpPixelFuncMetadata[] =
1168 : "<PixelFunctionArgumentsList>"
1169 : " <Argument name='base' description='Base' type='double' "
1170 : "default='2.7182818284590452353602874713526624' />"
1171 : " <Argument name='fact' description='Factor' type='double' default='1' />"
1172 : "</PixelFunctionArgumentsList>";
1173 :
1174 3 : static CPLErr ExpPixelFunc(void **papoSources, int nSources, void *pData,
1175 : int nXSize, int nYSize, GDALDataType eSrcType,
1176 : GDALDataType eBufType, int nPixelSpace,
1177 : int nLineSpace, CSLConstList papszArgs)
1178 : {
1179 3 : double dfBase = 2.7182818284590452353602874713526624;
1180 3 : double dfFact = 1.;
1181 :
1182 3 : if (FetchDoubleArg(papszArgs, "base", &dfBase, &dfBase) != CE_None)
1183 0 : return CE_Failure;
1184 :
1185 3 : if (FetchDoubleArg(papszArgs, "fact", &dfFact, &dfFact) != CE_None)
1186 0 : return CE_Failure;
1187 :
1188 3 : return ExpPixelFuncHelper(papoSources, nSources, pData, nXSize, nYSize,
1189 : eSrcType, eBufType, nPixelSpace, nLineSpace,
1190 3 : dfBase, dfFact);
1191 : } // ExpPixelFunc
1192 :
1193 2 : static CPLErr dB2AmpPixelFunc(void **papoSources, int nSources, void *pData,
1194 : int nXSize, int nYSize, GDALDataType eSrcType,
1195 : GDALDataType eBufType, int nPixelSpace,
1196 : int nLineSpace)
1197 : {
1198 2 : return ExpPixelFuncHelper(papoSources, nSources, pData, nXSize, nYSize,
1199 : eSrcType, eBufType, nPixelSpace, nLineSpace, 10.0,
1200 2 : 1. / 20);
1201 : } // dB2AmpPixelFunc
1202 :
1203 2 : static CPLErr dB2PowPixelFunc(void **papoSources, int nSources, void *pData,
1204 : int nXSize, int nYSize, GDALDataType eSrcType,
1205 : GDALDataType eBufType, int nPixelSpace,
1206 : int nLineSpace)
1207 : {
1208 2 : return ExpPixelFuncHelper(papoSources, nSources, pData, nXSize, nYSize,
1209 : eSrcType, eBufType, nPixelSpace, nLineSpace, 10.0,
1210 2 : 1. / 10);
1211 : } // dB2PowPixelFunc
1212 :
1213 : static const char pszPowPixelFuncMetadata[] =
1214 : "<PixelFunctionArgumentsList>"
1215 : " <Argument name='power' description='Exponent' type='double' "
1216 : "mandatory='1' />"
1217 : "</PixelFunctionArgumentsList>";
1218 :
1219 1 : static CPLErr PowPixelFunc(void **papoSources, int nSources, void *pData,
1220 : int nXSize, int nYSize, GDALDataType eSrcType,
1221 : GDALDataType eBufType, int nPixelSpace,
1222 : int nLineSpace, CSLConstList papszArgs)
1223 : {
1224 : /* ---- Init ---- */
1225 1 : if (nSources != 1)
1226 0 : return CE_Failure;
1227 1 : if (GDALDataTypeIsComplex(eSrcType))
1228 0 : return CE_Failure;
1229 :
1230 : double power;
1231 1 : if (FetchDoubleArg(papszArgs, "power", &power) != CE_None)
1232 0 : return CE_Failure;
1233 :
1234 : /* ---- Set pixels ---- */
1235 1 : size_t ii = 0;
1236 21 : for (int iLine = 0; iLine < nYSize; ++iLine)
1237 : {
1238 420 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1239 : {
1240 : const double dfPixVal =
1241 400 : std::pow(GetSrcVal(papoSources[0], eSrcType, ii), power);
1242 :
1243 400 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1244 : static_cast<GByte *>(pData) +
1245 400 : static_cast<GSpacing>(nLineSpace) * iLine +
1246 400 : iCol * nPixelSpace,
1247 : eBufType, nPixelSpace, 1);
1248 : }
1249 : }
1250 :
1251 : /* ---- Return success ---- */
1252 1 : return CE_None;
1253 : }
1254 :
1255 : // Given nt intervals spaced by dt and beginning at t0, return the index of
1256 : // the lower bound of the interval that should be used to
1257 : // interpolate/extrapolate a value for t.
1258 7 : static std::size_t intervalLeft(double t0, double dt, std::size_t nt, double t)
1259 : {
1260 7 : if (t < t0)
1261 : {
1262 2 : return 0;
1263 : }
1264 :
1265 5 : std::size_t n = static_cast<std::size_t>((t - t0) / dt);
1266 :
1267 5 : if (n >= nt - 1)
1268 : {
1269 3 : return nt - 2;
1270 : }
1271 :
1272 2 : return n;
1273 : }
1274 :
1275 16 : static double InterpolateLinear(double dfX0, double dfX1, double dfY0,
1276 : double dfY1, double dfX)
1277 : {
1278 16 : return dfY0 + (dfX - dfX0) * (dfY1 - dfY0) / (dfX1 - dfX0);
1279 : }
1280 :
1281 12 : static double InterpolateExponential(double dfX0, double dfX1, double dfY0,
1282 : double dfY1, double dfX)
1283 : {
1284 12 : const double r = std::log(dfY1 / dfY0) / (dfX1 - dfX0);
1285 12 : return dfY0 * std::exp(r * (dfX - dfX0));
1286 : }
1287 :
1288 : static const char pszInterpolatePixelFuncMetadata[] =
1289 : "<PixelFunctionArgumentsList>"
1290 : " <Argument name='t0' description='t0' type='double' mandatory='1' />"
1291 : " <Argument name='dt' description='dt' type='double' mandatory='1' />"
1292 : " <Argument name='t' description='t' type='double' mandatory='1' />"
1293 : "</PixelFunctionArgumentsList>";
1294 :
1295 : template <decltype(InterpolateLinear) InterpolationFunction>
1296 7 : CPLErr InterpolatePixelFunc(void **papoSources, int nSources, void *pData,
1297 : int nXSize, int nYSize, GDALDataType eSrcType,
1298 : GDALDataType eBufType, int nPixelSpace,
1299 : int nLineSpace, CSLConstList papszArgs)
1300 : {
1301 : /* ---- Init ---- */
1302 7 : if (GDALDataTypeIsComplex(eSrcType))
1303 0 : return CE_Failure;
1304 :
1305 : double dfT0;
1306 7 : if (FetchDoubleArg(papszArgs, "t0", &dfT0) == CE_Failure)
1307 0 : return CE_Failure;
1308 :
1309 : double dfT;
1310 7 : if (FetchDoubleArg(papszArgs, "t", &dfT) == CE_Failure)
1311 0 : return CE_Failure;
1312 :
1313 : double dfDt;
1314 7 : if (FetchDoubleArg(papszArgs, "dt", &dfDt) == CE_Failure)
1315 0 : return CE_Failure;
1316 :
1317 7 : if (nSources < 2)
1318 : {
1319 0 : CPLError(CE_Failure, CPLE_AppDefined,
1320 : "At least two sources required for interpolation.");
1321 0 : return CE_Failure;
1322 : }
1323 :
1324 7 : if (dfT == 0 || !std::isfinite(dfT))
1325 : {
1326 0 : CPLError(CE_Failure, CPLE_AppDefined, "dt must be finite and non-zero");
1327 0 : return CE_Failure;
1328 : }
1329 :
1330 7 : const auto i0 = intervalLeft(dfT0, dfDt, nSources, dfT);
1331 7 : const auto i1 = i0 + 1;
1332 7 : dfT0 = dfT0 + static_cast<double>(i0) * dfDt;
1333 7 : double dfX1 = dfT0 + dfDt;
1334 :
1335 : /* ---- Set pixels ---- */
1336 7 : size_t ii = 0;
1337 21 : for (int iLine = 0; iLine < nYSize; ++iLine)
1338 : {
1339 42 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1340 : {
1341 28 : const double dfY0 = GetSrcVal(papoSources[i0], eSrcType, ii);
1342 28 : const double dfY1 = GetSrcVal(papoSources[i1], eSrcType, ii);
1343 :
1344 28 : const double dfPixVal =
1345 28 : InterpolationFunction(dfT0, dfX1, dfY0, dfY1, dfT);
1346 :
1347 28 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1348 : static_cast<GByte *>(pData) +
1349 28 : static_cast<GSpacing>(nLineSpace) * iLine +
1350 28 : iCol * nPixelSpace,
1351 : eBufType, nPixelSpace, 1);
1352 : }
1353 : }
1354 :
1355 : /* ---- Return success ---- */
1356 7 : return CE_None;
1357 : }
1358 :
1359 : static const char pszReplaceNoDataPixelFuncMetadata[] =
1360 : "<PixelFunctionArgumentsList>"
1361 : " <Argument type='builtin' value='NoData' />"
1362 : " <Argument name='to' type='double' description='New NoData value to be "
1363 : "replaced' default='nan' />"
1364 : "</PixelFunctionArgumentsList>";
1365 :
1366 2 : static CPLErr ReplaceNoDataPixelFunc(void **papoSources, int nSources,
1367 : void *pData, int nXSize, int nYSize,
1368 : GDALDataType eSrcType,
1369 : GDALDataType eBufType, int nPixelSpace,
1370 : int nLineSpace, CSLConstList papszArgs)
1371 : {
1372 : /* ---- Init ---- */
1373 2 : if (nSources != 1)
1374 0 : return CE_Failure;
1375 2 : if (GDALDataTypeIsComplex(eSrcType))
1376 : {
1377 0 : CPLError(CE_Failure, CPLE_AppDefined,
1378 : "replace_nodata cannot convert complex data types");
1379 0 : return CE_Failure;
1380 : }
1381 :
1382 2 : double dfOldNoData, dfNewNoData = NAN;
1383 2 : if (FetchDoubleArg(papszArgs, "NoData", &dfOldNoData) != CE_None)
1384 0 : return CE_Failure;
1385 2 : if (FetchDoubleArg(papszArgs, "to", &dfNewNoData, &dfNewNoData) != CE_None)
1386 0 : return CE_Failure;
1387 :
1388 2 : if (!GDALDataTypeIsFloating(eBufType) && std::isnan(dfNewNoData))
1389 : {
1390 0 : CPLError(CE_Failure, CPLE_AppDefined,
1391 : "Using nan requires a floating point type output buffer");
1392 0 : return CE_Failure;
1393 : }
1394 :
1395 : /* ---- Set pixels ---- */
1396 2 : size_t ii = 0;
1397 102 : for (int iLine = 0; iLine < nYSize; ++iLine)
1398 : {
1399 5100 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1400 : {
1401 5000 : double dfPixVal = GetSrcVal(papoSources[0], eSrcType, ii);
1402 5000 : if (dfPixVal == dfOldNoData || std::isnan(dfPixVal))
1403 3200 : dfPixVal = dfNewNoData;
1404 :
1405 5000 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1406 : static_cast<GByte *>(pData) +
1407 5000 : static_cast<GSpacing>(nLineSpace) * iLine +
1408 5000 : iCol * nPixelSpace,
1409 : eBufType, nPixelSpace, 1);
1410 : }
1411 : }
1412 :
1413 : /* ---- Return success ---- */
1414 2 : return CE_None;
1415 : }
1416 :
1417 : static const char pszScalePixelFuncMetadata[] =
1418 : "<PixelFunctionArgumentsList>"
1419 : " <Argument type='builtin' value='offset' />"
1420 : " <Argument type='builtin' value='scale' />"
1421 : "</PixelFunctionArgumentsList>";
1422 :
1423 1 : static CPLErr ScalePixelFunc(void **papoSources, int nSources, void *pData,
1424 : int nXSize, int nYSize, GDALDataType eSrcType,
1425 : GDALDataType eBufType, int nPixelSpace,
1426 : int nLineSpace, CSLConstList papszArgs)
1427 : {
1428 : /* ---- Init ---- */
1429 1 : if (nSources != 1)
1430 0 : return CE_Failure;
1431 1 : if (GDALDataTypeIsComplex(eSrcType))
1432 : {
1433 0 : CPLError(CE_Failure, CPLE_AppDefined,
1434 : "scale cannot by applied to complex data types");
1435 0 : return CE_Failure;
1436 : }
1437 :
1438 : double dfScale, dfOffset;
1439 1 : if (FetchDoubleArg(papszArgs, "scale", &dfScale) != CE_None)
1440 0 : return CE_Failure;
1441 1 : if (FetchDoubleArg(papszArgs, "offset", &dfOffset) != CE_None)
1442 0 : return CE_Failure;
1443 :
1444 : /* ---- Set pixels ---- */
1445 1 : size_t ii = 0;
1446 21 : for (int iLine = 0; iLine < nYSize; ++iLine)
1447 : {
1448 420 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1449 : {
1450 : const double dfPixVal =
1451 400 : GetSrcVal(papoSources[0], eSrcType, ii) * dfScale + dfOffset;
1452 :
1453 400 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1454 : static_cast<GByte *>(pData) +
1455 400 : static_cast<GSpacing>(nLineSpace) * iLine +
1456 400 : iCol * nPixelSpace,
1457 : eBufType, nPixelSpace, 1);
1458 : }
1459 : }
1460 :
1461 : /* ---- Return success ---- */
1462 1 : return CE_None;
1463 : }
1464 :
1465 1 : static CPLErr NormDiffPixelFunc(void **papoSources, int nSources, void *pData,
1466 : int nXSize, int nYSize, GDALDataType eSrcType,
1467 : GDALDataType eBufType, int nPixelSpace,
1468 : int nLineSpace)
1469 : {
1470 : /* ---- Init ---- */
1471 1 : if (nSources != 2)
1472 0 : return CE_Failure;
1473 :
1474 1 : if (GDALDataTypeIsComplex(eSrcType))
1475 : {
1476 0 : CPLError(CE_Failure, CPLE_AppDefined,
1477 : "norm_diff cannot by applied to complex data types");
1478 0 : return CE_Failure;
1479 : }
1480 :
1481 : /* ---- Set pixels ---- */
1482 1 : size_t ii = 0;
1483 7 : for (int iLine = 0; iLine < nYSize; ++iLine)
1484 : {
1485 36 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1486 : {
1487 30 : const double dfLeftVal = GetSrcVal(papoSources[0], eSrcType, ii);
1488 30 : const double dfRightVal = GetSrcVal(papoSources[1], eSrcType, ii);
1489 :
1490 30 : const double dfDenom = (dfLeftVal + dfRightVal);
1491 :
1492 : // coverity[divide_by_zero]
1493 : const double dfPixVal =
1494 30 : dfDenom == 0 ? std::numeric_limits<double>::infinity()
1495 30 : : (dfLeftVal - dfRightVal) / dfDenom;
1496 :
1497 30 : GDALCopyWords(&dfPixVal, GDT_Float64, 0,
1498 : static_cast<GByte *>(pData) +
1499 30 : static_cast<GSpacing>(nLineSpace) * iLine +
1500 30 : iCol * nPixelSpace,
1501 : eBufType, nPixelSpace, 1);
1502 : }
1503 : }
1504 :
1505 : /* ---- Return success ---- */
1506 1 : return CE_None;
1507 : } // NormDiffPixelFunc
1508 :
1509 : /************************************************************************/
1510 : /* pszMinMaxFuncMetadataNodata */
1511 : /************************************************************************/
1512 :
1513 : static const char pszMinMaxFuncMetadataNodata[] =
1514 : "<PixelFunctionArgumentsList>"
1515 : " <Argument type='builtin' value='NoData' optional='true' />"
1516 : " <Argument name='propagateNoData' description='Whether the output value "
1517 : "should be NoData as as soon as one source is NoData' type='boolean' "
1518 : "default='false' />"
1519 : "</PixelFunctionArgumentsList>";
1520 :
1521 : template <class Comparator>
1522 5 : static CPLErr MinOrMaxPixelFunc(void **papoSources, int nSources, void *pData,
1523 : int nXSize, int nYSize, GDALDataType eSrcType,
1524 : GDALDataType eBufType, int nPixelSpace,
1525 : int nLineSpace, CSLConstList papszArgs)
1526 : {
1527 : /* ---- Init ---- */
1528 5 : if (nSources < 2)
1529 0 : return CE_Failure;
1530 :
1531 5 : if (GDALDataTypeIsComplex(eSrcType))
1532 : {
1533 0 : CPLError(CE_Failure, CPLE_AppDefined,
1534 : "Complex data type not supported for min/max().");
1535 0 : return CE_Failure;
1536 : }
1537 :
1538 5 : double dfNoData = std::numeric_limits<double>::quiet_NaN();
1539 5 : if (FetchDoubleArg(papszArgs, "NoData", &dfNoData, &dfNoData) != CE_None)
1540 0 : return CE_Failure;
1541 5 : const bool bPropagateNoData = CPLTestBool(
1542 : CSLFetchNameValueDef(papszArgs, "propagateNoData", "false"));
1543 :
1544 : /* ---- Set pixels ---- */
1545 5 : size_t ii = 0;
1546 255 : for (int iLine = 0; iLine < nYSize; ++iLine)
1547 : {
1548 12750 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1549 : {
1550 12500 : double dfRes = std::numeric_limits<double>::quiet_NaN();
1551 :
1552 43400 : for (int iSrc = 0; iSrc < nSources; ++iSrc)
1553 : {
1554 32950 : const double dfVal = GetSrcVal(papoSources[iSrc], eSrcType, ii);
1555 :
1556 32950 : if (std::isnan(dfVal) || dfVal == dfNoData)
1557 : {
1558 14350 : if (bPropagateNoData)
1559 : {
1560 2050 : dfRes = dfNoData;
1561 2050 : break;
1562 : }
1563 : }
1564 18600 : else if (Comparator::compare(dfVal, dfRes))
1565 : {
1566 9950 : dfRes = dfVal;
1567 : }
1568 : }
1569 :
1570 12500 : if (!bPropagateNoData && std::isnan(dfRes))
1571 : {
1572 3200 : dfRes = dfNoData;
1573 : }
1574 :
1575 12500 : GDALCopyWords(&dfRes, GDT_Float64, 0,
1576 : static_cast<GByte *>(pData) +
1577 12500 : static_cast<GSpacing>(nLineSpace) * iLine +
1578 12500 : iCol * nPixelSpace,
1579 : eBufType, nPixelSpace, 1);
1580 : }
1581 : }
1582 :
1583 : /* ---- Return success ---- */
1584 5 : return CE_None;
1585 : } /* MinOrMaxPixelFunc */
1586 :
1587 2 : static CPLErr MinPixelFunc(void **papoSources, int nSources, void *pData,
1588 : int nXSize, int nYSize, GDALDataType eSrcType,
1589 : GDALDataType eBufType, int nPixelSpace,
1590 : int nLineSpace, CSLConstList papszArgs)
1591 : {
1592 : struct Comparator
1593 : {
1594 2700 : static bool compare(double x, double resVal)
1595 : {
1596 : // Written this way to deal with resVal being NaN
1597 2700 : return !(x >= resVal);
1598 : }
1599 : };
1600 :
1601 2 : return MinOrMaxPixelFunc<Comparator>(papoSources, nSources, pData, nXSize,
1602 : nYSize, eSrcType, eBufType,
1603 2 : nPixelSpace, nLineSpace, papszArgs);
1604 : }
1605 :
1606 3 : static CPLErr MaxPixelFunc(void **papoSources, int nSources, void *pData,
1607 : int nXSize, int nYSize, GDALDataType eSrcType,
1608 : GDALDataType eBufType, int nPixelSpace,
1609 : int nLineSpace, CSLConstList papszArgs)
1610 : {
1611 : struct Comparator
1612 : {
1613 15900 : static bool compare(double x, double resVal)
1614 : {
1615 : // Written this way to deal with resVal being NaN
1616 15900 : return !(x <= resVal);
1617 : }
1618 : };
1619 :
1620 3 : return MinOrMaxPixelFunc<Comparator>(papoSources, nSources, pData, nXSize,
1621 : nYSize, eSrcType, eBufType,
1622 3 : nPixelSpace, nLineSpace, papszArgs);
1623 : }
1624 :
1625 : static const char pszExprPixelFuncMetadata[] =
1626 : "<PixelFunctionArgumentsList>"
1627 : " <Argument name='expression' "
1628 : " description='Expression to be evaluated' "
1629 : " type='string'></Argument>"
1630 : " <Argument name='dialect' "
1631 : " description='Expression dialect' "
1632 : " type='string-select'"
1633 : " default='exprtk'>"
1634 : " <Value>exprtk</Value>"
1635 : " <Value>muparser</Value>"
1636 : " </Argument>"
1637 : "</PixelFunctionArgumentsList>";
1638 :
1639 70 : static CPLErr ExprPixelFunc(void **papoSources, int nSources, void *pData,
1640 : int nXSize, int nYSize, GDALDataType eSrcType,
1641 : GDALDataType eBufType, int nPixelSpace,
1642 : int nLineSpace, CSLConstList papszArgs)
1643 : {
1644 : /* ---- Init ---- */
1645 :
1646 70 : if (GDALDataTypeIsComplex(eSrcType))
1647 : {
1648 0 : CPLError(CE_Failure, CPLE_AppDefined,
1649 : "expression cannot by applied to complex data types");
1650 0 : return CE_Failure;
1651 : }
1652 :
1653 70 : std::unique_ptr<gdal::MathExpression> poExpression;
1654 :
1655 70 : const char *pszExpression = CSLFetchNameValue(papszArgs, "expression");
1656 :
1657 70 : const char *pszSourceNames = CSLFetchNameValue(papszArgs, "SOURCE_NAMES");
1658 : const CPLStringList aosSourceNames(
1659 140 : CSLTokenizeString2(pszSourceNames, "|", 0));
1660 :
1661 140 : std::vector<double> adfValuesForPixel(nSources);
1662 :
1663 70 : const char *pszDialect = CSLFetchNameValue(papszArgs, "dialect");
1664 70 : if (!pszDialect)
1665 : {
1666 0 : pszDialect = "muparser";
1667 : }
1668 :
1669 70 : poExpression = gdal::MathExpression::Create(pszExpression, pszDialect);
1670 :
1671 : // cppcheck-suppress knownConditionTrueFalse
1672 70 : if (!poExpression)
1673 : {
1674 0 : return CE_Failure;
1675 : }
1676 :
1677 : {
1678 70 : int iSource = 0;
1679 216 : for (const auto &osName : aosSourceNames)
1680 : {
1681 292 : poExpression->RegisterVariable(osName,
1682 146 : &adfValuesForPixel[iSource++]);
1683 : }
1684 : }
1685 :
1686 70 : if (strstr(pszExpression, "BANDS"))
1687 : {
1688 1 : poExpression->RegisterVector("BANDS", &adfValuesForPixel);
1689 : }
1690 :
1691 : std::unique_ptr<double, VSIFreeReleaser> padfResults(
1692 140 : static_cast<double *>(VSI_MALLOC2_VERBOSE(nXSize, sizeof(double))));
1693 70 : if (!padfResults)
1694 0 : return CE_Failure;
1695 :
1696 : /* ---- Set pixels ---- */
1697 70 : size_t ii = 0;
1698 4203 : for (int iLine = 0; iLine < nYSize; ++iLine)
1699 : {
1700 12982900 : for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
1701 : {
1702 38947700 : for (int iSrc = 0; iSrc < nSources; iSrc++)
1703 : {
1704 : // cppcheck-suppress unreadVariable
1705 25968900 : adfValuesForPixel[iSrc] =
1706 25968900 : GetSrcVal(papoSources[iSrc], eSrcType, ii);
1707 : }
1708 :
1709 12978700 : if (auto eErr = poExpression->Evaluate(); eErr != CE_None)
1710 : {
1711 2 : return CE_Failure;
1712 : }
1713 : else
1714 : {
1715 12978700 : padfResults.get()[iCol] = poExpression->Results()[0];
1716 : }
1717 : }
1718 :
1719 4133 : GDALCopyWords(padfResults.get(), GDT_Float64, sizeof(double),
1720 : static_cast<GByte *>(pData) +
1721 4133 : static_cast<GSpacing>(nLineSpace) * iLine,
1722 : eBufType, nPixelSpace, nXSize);
1723 : }
1724 :
1725 : /* ---- Return success ---- */
1726 68 : return CE_None;
1727 : } // ExprPixelFunc
1728 :
1729 : /************************************************************************/
1730 : /* GDALRegisterDefaultPixelFunc() */
1731 : /************************************************************************/
1732 :
1733 : /**
1734 : * This adds a default set of pixel functions to the global list of
1735 : * available pixel functions for derived bands:
1736 : *
1737 : * - "real": extract real part from a single raster band (just a copy if the
1738 : * input is non-complex)
1739 : * - "imag": extract imaginary part from a single raster band (0 for
1740 : * non-complex)
1741 : * - "complex": make a complex band merging two bands used as real and
1742 : * imag values
1743 : * - "polar": make a complex band using input bands for amplitude and
1744 : * phase values (b1 * exp( j * b2 ))
1745 : * - "mod": extract module from a single raster band (real or complex)
1746 : * - "phase": extract phase from a single raster band [-PI,PI] (0 or PI for
1747 : non-complex)
1748 : * - "conj": computes the complex conjugate of a single raster band (just a
1749 : * copy if the input is non-complex)
1750 : * - "sum": sum 2 or more raster bands
1751 : * - "diff": computes the difference between 2 raster bands (b1 - b2)
1752 : * - "mul": multiply 2 or more raster bands
1753 : * - "div": divide one raster band by another (b1 / b2).
1754 : * - "min": minimum value of 2 or more raster bands
1755 : * - "max": maximum value of 2 or more raster bands
1756 : * - "norm_diff": computes the normalized difference between two raster bands:
1757 : * ``(b1 - b2)/(b1 + b2)``.
1758 : * - "cmul": multiply the first band for the complex conjugate of the second
1759 : * - "inv": inverse (1./x).
1760 : * - "intensity": computes the intensity Re(x*conj(x)) of a single raster band
1761 : * (real or complex)
1762 : * - "sqrt": perform the square root of a single raster band (real only)
1763 : * - "log10": compute the logarithm (base 10) of the abs of a single raster
1764 : * band (real or complex): log10( abs( x ) )
1765 : * - "dB": perform conversion to dB of the abs of a single raster
1766 : * band (real or complex): 20. * log10( abs( x ) ).
1767 : * Note: the optional fact parameter can be set to 10. to get the
1768 : * alternative formula: 10. * log10( abs( x ) )
1769 : * - "exp": computes the exponential of each element in the input band ``x``
1770 : * (of real values): ``e ^ x``.
1771 : * The function also accepts two optional parameters: ``base`` and
1772 : ``fact``
1773 : * that allow to compute the generalized formula: ``base ^ ( fact *
1774 : x)``.
1775 : * Note: this function is the recommended one to perform conversion
1776 : * form logarithmic scale (dB): `` 10. ^ (x / 20.)``, in this case
1777 : * ``base = 10.`` and ``fact = 1./20``
1778 : * - "dB2amp": perform scale conversion from logarithmic to linear
1779 : * (amplitude) (i.e. 10 ^ ( x / 20 ) ) of a single raster
1780 : * band (real only).
1781 : * Deprecated in GDAL v3.5. Please use the ``exp`` pixel function
1782 : with
1783 : * ``base = 10.`` and ``fact = 0.05`` i.e. ``1./20``
1784 : * - "dB2pow": perform scale conversion from logarithmic to linear
1785 : * (power) (i.e. 10 ^ ( x / 10 ) ) of a single raster
1786 : * band (real only)
1787 : * Deprecated in GDAL v3.5. Please use the ``exp`` pixel function
1788 : with
1789 : * ``base = 10.`` and ``fact = 0.1`` i.e. ``1./10``
1790 : * - "pow": raise a single raster band to a constant power
1791 : * - "interpolate_linear": interpolate values between two raster bands
1792 : * using linear interpolation
1793 : * - "interpolate_exp": interpolate values between two raster bands using
1794 : * exponential interpolation
1795 : * - "scale": Apply the RasterBand metadata values of "offset" and "scale"
1796 : * - "nan": Convert incoming NoData values to IEEE 754 nan
1797 : *
1798 : * @see GDALAddDerivedBandPixelFunc
1799 : *
1800 : * @return CE_None
1801 : */
1802 1384 : CPLErr GDALRegisterDefaultPixelFunc()
1803 : {
1804 1384 : GDALAddDerivedBandPixelFunc("real", RealPixelFunc);
1805 1384 : GDALAddDerivedBandPixelFunc("imag", ImagPixelFunc);
1806 1384 : GDALAddDerivedBandPixelFunc("complex", ComplexPixelFunc);
1807 1384 : GDALAddDerivedBandPixelFuncWithArgs("polar", PolarPixelFunc,
1808 : pszPolarPixelFuncMetadata);
1809 1384 : GDALAddDerivedBandPixelFunc("mod", ModulePixelFunc);
1810 1384 : GDALAddDerivedBandPixelFunc("phase", PhasePixelFunc);
1811 1384 : GDALAddDerivedBandPixelFunc("conj", ConjPixelFunc);
1812 1384 : GDALAddDerivedBandPixelFuncWithArgs("sum", SumPixelFunc,
1813 : pszSumPixelFuncMetadata);
1814 1384 : GDALAddDerivedBandPixelFunc("diff", DiffPixelFunc);
1815 1384 : GDALAddDerivedBandPixelFuncWithArgs("mul", MulPixelFunc,
1816 : pszMulPixelFuncMetadata);
1817 1384 : GDALAddDerivedBandPixelFunc("div", DivPixelFunc);
1818 1384 : GDALAddDerivedBandPixelFunc("cmul", CMulPixelFunc);
1819 1384 : GDALAddDerivedBandPixelFuncWithArgs("inv", InvPixelFunc,
1820 : pszInvPixelFuncMetadata);
1821 1384 : GDALAddDerivedBandPixelFunc("intensity", IntensityPixelFunc);
1822 1384 : GDALAddDerivedBandPixelFunc("sqrt", SqrtPixelFunc);
1823 1384 : GDALAddDerivedBandPixelFunc("log10", Log10PixelFunc);
1824 1384 : GDALAddDerivedBandPixelFuncWithArgs("dB", DBPixelFunc,
1825 : pszDBPixelFuncMetadata);
1826 1384 : GDALAddDerivedBandPixelFuncWithArgs("exp", ExpPixelFunc,
1827 : pszExpPixelFuncMetadata);
1828 1384 : GDALAddDerivedBandPixelFunc("dB2amp",
1829 : dB2AmpPixelFunc); // deprecated in v3.5
1830 1384 : GDALAddDerivedBandPixelFunc("dB2pow",
1831 : dB2PowPixelFunc); // deprecated in v3.5
1832 1384 : GDALAddDerivedBandPixelFuncWithArgs("pow", PowPixelFunc,
1833 : pszPowPixelFuncMetadata);
1834 1384 : GDALAddDerivedBandPixelFuncWithArgs("interpolate_linear",
1835 : InterpolatePixelFunc<InterpolateLinear>,
1836 : pszInterpolatePixelFuncMetadata);
1837 1384 : GDALAddDerivedBandPixelFuncWithArgs(
1838 : "interpolate_exp", InterpolatePixelFunc<InterpolateExponential>,
1839 : pszInterpolatePixelFuncMetadata);
1840 1384 : GDALAddDerivedBandPixelFuncWithArgs("replace_nodata",
1841 : ReplaceNoDataPixelFunc,
1842 : pszReplaceNoDataPixelFuncMetadata);
1843 1384 : GDALAddDerivedBandPixelFuncWithArgs("scale", ScalePixelFunc,
1844 : pszScalePixelFuncMetadata);
1845 1384 : GDALAddDerivedBandPixelFunc("norm_diff", NormDiffPixelFunc);
1846 1384 : GDALAddDerivedBandPixelFuncWithArgs("min", MinPixelFunc,
1847 : pszMinMaxFuncMetadataNodata);
1848 1384 : GDALAddDerivedBandPixelFuncWithArgs("max", MaxPixelFunc,
1849 : pszMinMaxFuncMetadataNodata);
1850 1384 : GDALAddDerivedBandPixelFuncWithArgs("expression", ExprPixelFunc,
1851 : pszExprPixelFuncMetadata);
1852 1384 : return CE_None;
1853 : }
|
