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