Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Virtual GDAL Datasets
4 : * Purpose: Implementation of VRTProcessedDataset processing functions
5 : * Author: Even Rouault <even.rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2024, Even Rouault <even.rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "cpl_minixml.h"
14 : #include "cpl_string.h"
15 : #include "vrtdataset.h"
16 : #include "vrtexpression.h"
17 :
18 : #include <algorithm>
19 : #include <functional>
20 : #include <limits>
21 : #include <map>
22 : #include <optional>
23 : #include <set>
24 : #include <vector>
25 :
26 : /************************************************************************/
27 : /* GetDstValue() */
28 : /************************************************************************/
29 :
30 : /** Return a destination value given an initial value, the destination no data
31 : * value and its replacement value
32 : */
33 3178 : static inline double GetDstValue(double dfVal, double dfDstNoData,
34 : double dfReplacementDstNodata,
35 : GDALDataType eIntendedDstDT,
36 : bool bDstIntendedDTIsInteger)
37 : {
38 3178 : if (bDstIntendedDTIsInteger && std::round(dfVal) == dfDstNoData)
39 : {
40 1 : return dfReplacementDstNodata;
41 : }
42 3177 : else if (eIntendedDstDT == GDT_Float16 &&
43 3177 : static_cast<GFloat16>(dfVal) == static_cast<GFloat16>(dfDstNoData))
44 : {
45 0 : return dfReplacementDstNodata;
46 : }
47 3177 : else if (eIntendedDstDT == GDT_Float32 &&
48 0 : static_cast<float>(dfVal) == static_cast<float>(dfDstNoData))
49 : {
50 0 : return dfReplacementDstNodata;
51 : }
52 3177 : else if (eIntendedDstDT == GDT_Float64 && dfVal == dfDstNoData)
53 : {
54 1 : return dfReplacementDstNodata;
55 : }
56 : else
57 : {
58 3176 : return dfVal;
59 : }
60 : }
61 :
62 : /************************************************************************/
63 : /* BandAffineCombinationData */
64 : /************************************************************************/
65 :
66 : namespace
67 : {
68 : /** Working structure for 'BandAffineCombination' builtin function. */
69 : struct BandAffineCombinationData
70 : {
71 : static constexpr const char *const EXPECTED_SIGNATURE =
72 : "BandAffineCombination";
73 : //! Signature (to make sure callback functions are called with the right argument)
74 : const std::string m_osSignature = EXPECTED_SIGNATURE;
75 :
76 : /** Replacement nodata value */
77 : std::vector<double> m_adfReplacementDstNodata{};
78 :
79 : /** Intended destination data type. */
80 : GDALDataType m_eIntendedDstDT = GDT_Float64;
81 :
82 : /** Affine transformation coefficients.
83 : * m_aadfCoefficients[i][0] is the constant term for the i(th) dst band
84 : * m_aadfCoefficients[i][j] is the weight of the j(th) src band for the
85 : * i(th) dst vand.
86 : * Said otherwise dst[i] = m_aadfCoefficients[i][0] +
87 : * sum(m_aadfCoefficients[i][j + 1] * src[j] for j in 0...nSrcBands-1)
88 : */
89 : std::vector<std::vector<double>> m_aadfCoefficients{};
90 :
91 : //! Minimum clamping value.
92 : double m_dfClampMin = std::numeric_limits<double>::quiet_NaN();
93 :
94 : //! Maximum clamping value.
95 : double m_dfClampMax = std::numeric_limits<double>::quiet_NaN();
96 : };
97 : } // namespace
98 :
99 : /************************************************************************/
100 : /* SetOutputValuesForInNoDataAndOutNoData() */
101 : /************************************************************************/
102 :
103 44 : static std::vector<double> SetOutputValuesForInNoDataAndOutNoData(
104 : int nInBands, double *padfInNoData, int *pnOutBands,
105 : double **ppadfOutNoData, bool bSrcNodataSpecified, double dfSrcNoData,
106 : bool bDstNodataSpecified, double dfDstNoData, bool bIsFinalStep)
107 : {
108 44 : if (bSrcNodataSpecified)
109 : {
110 3 : std::vector<double> adfNoData(nInBands, dfSrcNoData);
111 3 : memcpy(padfInNoData, adfNoData.data(),
112 3 : adfNoData.size() * sizeof(double));
113 : }
114 :
115 44 : std::vector<double> adfDstNoData;
116 44 : if (bDstNodataSpecified)
117 : {
118 3 : adfDstNoData.resize(*pnOutBands, dfDstNoData);
119 : }
120 41 : else if (bIsFinalStep)
121 : {
122 : adfDstNoData =
123 35 : std::vector<double>(*ppadfOutNoData, *ppadfOutNoData + *pnOutBands);
124 : }
125 : else
126 : {
127 : adfDstNoData =
128 6 : std::vector<double>(padfInNoData, padfInNoData + nInBands);
129 6 : adfDstNoData.resize(*pnOutBands, *padfInNoData);
130 : }
131 :
132 44 : if (*ppadfOutNoData == nullptr)
133 : {
134 6 : *ppadfOutNoData =
135 6 : static_cast<double *>(CPLMalloc(*pnOutBands * sizeof(double)));
136 : }
137 44 : memcpy(*ppadfOutNoData, adfDstNoData.data(), *pnOutBands * sizeof(double));
138 :
139 44 : return adfDstNoData;
140 : }
141 :
142 : /************************************************************************/
143 : /* BandAffineCombinationInit() */
144 : /************************************************************************/
145 :
146 : /** Init function for 'BandAffineCombination' builtin function. */
147 38 : static CPLErr BandAffineCombinationInit(
148 : const char * /*pszFuncName*/, void * /*pUserData*/,
149 : CSLConstList papszFunctionArgs, int nInBands, GDALDataType eInDT,
150 : double *padfInNoData, int *pnOutBands, GDALDataType *peOutDT,
151 : double **ppadfOutNoData, const char * /* pszVRTPath */,
152 : VRTPDWorkingDataPtr *ppWorkingData)
153 : {
154 38 : CPLAssert(eInDT == GDT_Float64);
155 :
156 38 : *peOutDT = eInDT;
157 38 : *ppWorkingData = nullptr;
158 :
159 76 : auto data = std::make_unique<BandAffineCombinationData>();
160 :
161 76 : std::map<int, std::vector<double>> oMapCoefficients{};
162 38 : double dfSrcNoData = std::numeric_limits<double>::quiet_NaN();
163 38 : bool bSrcNodataSpecified = false;
164 38 : double dfDstNoData = std::numeric_limits<double>::quiet_NaN();
165 38 : bool bDstNodataSpecified = false;
166 38 : double dfReplacementDstNodata = std::numeric_limits<double>::quiet_NaN();
167 38 : bool bReplacementDstNodataSpecified = false;
168 :
169 195 : for (const auto &[pszKey, pszValue] :
170 232 : cpl::IterateNameValue(papszFunctionArgs))
171 : {
172 98 : if (EQUAL(pszKey, "src_nodata"))
173 : {
174 2 : bSrcNodataSpecified = true;
175 2 : dfSrcNoData = CPLAtof(pszValue);
176 : }
177 96 : else if (EQUAL(pszKey, "dst_nodata"))
178 : {
179 2 : bDstNodataSpecified = true;
180 2 : dfDstNoData = CPLAtof(pszValue);
181 : }
182 94 : else if (EQUAL(pszKey, "replacement_nodata"))
183 : {
184 1 : bReplacementDstNodataSpecified = true;
185 1 : dfReplacementDstNodata = CPLAtof(pszValue);
186 : }
187 93 : else if (EQUAL(pszKey, "dst_intended_datatype"))
188 : {
189 1 : for (GDALDataType eDT = GDT_Byte; eDT < GDT_TypeCount;
190 0 : eDT = static_cast<GDALDataType>(eDT + 1))
191 : {
192 1 : if (EQUAL(GDALGetDataTypeName(eDT), pszValue))
193 : {
194 1 : data->m_eIntendedDstDT = eDT;
195 1 : break;
196 : }
197 : }
198 : }
199 92 : else if (STARTS_WITH_CI(pszKey, "coefficients_"))
200 : {
201 88 : const int nTargetBand = atoi(pszKey + strlen("coefficients_"));
202 88 : if (nTargetBand <= 0 || nTargetBand > 65536)
203 : {
204 0 : CPLError(CE_Failure, CPLE_AppDefined,
205 : "Invalid band in argument '%s'", pszKey);
206 1 : return CE_Failure;
207 : }
208 88 : const CPLStringList aosTokens(CSLTokenizeString2(pszValue, ",", 0));
209 88 : if (aosTokens.size() != 1 + nInBands)
210 : {
211 1 : CPLError(CE_Failure, CPLE_AppDefined,
212 : "Argument %s has %d values, whereas %d are expected",
213 : pszKey, aosTokens.size(), 1 + nInBands);
214 1 : return CE_Failure;
215 : }
216 87 : std::vector<double> adfValues;
217 401 : for (int i = 0; i < aosTokens.size(); ++i)
218 : {
219 314 : adfValues.push_back(CPLAtof(aosTokens[i]));
220 : }
221 87 : oMapCoefficients[nTargetBand - 1] = std::move(adfValues);
222 : }
223 4 : else if (EQUAL(pszKey, "min"))
224 : {
225 2 : data->m_dfClampMin = CPLAtof(pszValue);
226 : }
227 2 : else if (EQUAL(pszKey, "max"))
228 : {
229 2 : data->m_dfClampMax = CPLAtof(pszValue);
230 : }
231 : else
232 : {
233 0 : CPLError(CE_Warning, CPLE_AppDefined,
234 : "Unrecognized argument name %s. Ignored", pszKey);
235 : }
236 : }
237 :
238 37 : const bool bIsFinalStep = *pnOutBands != 0;
239 37 : if (bIsFinalStep)
240 : {
241 31 : if (*pnOutBands != static_cast<int>(oMapCoefficients.size()))
242 : {
243 2 : CPLError(CE_Failure, CPLE_AppDefined,
244 : "Final step expect %d bands, but only %d coefficient_XX "
245 : "are provided",
246 2 : *pnOutBands, static_cast<int>(oMapCoefficients.size()));
247 2 : return CE_Failure;
248 : }
249 : }
250 : else
251 : {
252 6 : *pnOutBands = static_cast<int>(oMapCoefficients.size());
253 : }
254 :
255 : const std::vector<double> adfDstNoData =
256 : SetOutputValuesForInNoDataAndOutNoData(
257 : nInBands, padfInNoData, pnOutBands, ppadfOutNoData,
258 : bSrcNodataSpecified, dfSrcNoData, bDstNodataSpecified, dfDstNoData,
259 70 : bIsFinalStep);
260 :
261 35 : if (bReplacementDstNodataSpecified)
262 : {
263 1 : data->m_adfReplacementDstNodata.resize(*pnOutBands,
264 : dfReplacementDstNodata);
265 : }
266 : else
267 : {
268 116 : for (double dfVal : adfDstNoData)
269 : {
270 82 : data->m_adfReplacementDstNodata.emplace_back(
271 82 : GDALGetNoDataReplacementValue(data->m_eIntendedDstDT, dfVal));
272 : }
273 : }
274 :
275 : // Check we have a set of coefficient for all output bands and
276 : // convert the map to a vector
277 117 : for (auto &oIter : oMapCoefficients)
278 : {
279 84 : const int iExpected = static_cast<int>(data->m_aadfCoefficients.size());
280 84 : if (oIter.first != iExpected)
281 : {
282 2 : CPLError(CE_Failure, CPLE_AppDefined,
283 : "Argument coefficients_%d is missing", iExpected + 1);
284 2 : return CE_Failure;
285 : }
286 82 : data->m_aadfCoefficients.emplace_back(std::move(oIter.second));
287 : }
288 33 : *ppWorkingData = data.release();
289 33 : return CE_None;
290 : }
291 :
292 : /************************************************************************/
293 : /* BandAffineCombinationFree() */
294 : /************************************************************************/
295 :
296 : /** Free function for 'BandAffineCombination' builtin function. */
297 33 : static void BandAffineCombinationFree(const char * /*pszFuncName*/,
298 : void * /*pUserData*/,
299 : VRTPDWorkingDataPtr pWorkingData)
300 : {
301 33 : BandAffineCombinationData *data =
302 : static_cast<BandAffineCombinationData *>(pWorkingData);
303 33 : CPLAssert(data->m_osSignature ==
304 : BandAffineCombinationData::EXPECTED_SIGNATURE);
305 33 : CPL_IGNORE_RET_VAL(data->m_osSignature);
306 33 : delete data;
307 33 : }
308 :
309 : /************************************************************************/
310 : /* BandAffineCombinationProcess() */
311 : /************************************************************************/
312 :
313 : /** Processing function for 'BandAffineCombination' builtin function. */
314 41 : static CPLErr BandAffineCombinationProcess(
315 : const char * /*pszFuncName*/, void * /*pUserData*/,
316 : VRTPDWorkingDataPtr pWorkingData, CSLConstList /* papszFunctionArgs*/,
317 : int nBufXSize, int nBufYSize, const void *pInBuffer, size_t nInBufferSize,
318 : GDALDataType eInDT, int nInBands, const double *CPL_RESTRICT padfInNoData,
319 : void *pOutBuffer, size_t nOutBufferSize, GDALDataType eOutDT, int nOutBands,
320 : const double *CPL_RESTRICT padfOutNoData, double /*dfSrcXOff*/,
321 : double /*dfSrcYOff*/, double /*dfSrcXSize*/, double /*dfSrcYSize*/,
322 : const double /*adfSrcGT*/[], const char * /* pszVRTPath */,
323 : CSLConstList /*papszExtra*/)
324 : {
325 41 : const size_t nElts = static_cast<size_t>(nBufXSize) * nBufYSize;
326 :
327 41 : CPL_IGNORE_RET_VAL(eInDT);
328 41 : CPLAssert(eInDT == GDT_Float64);
329 41 : CPL_IGNORE_RET_VAL(eOutDT);
330 41 : CPLAssert(eOutDT == GDT_Float64);
331 41 : CPL_IGNORE_RET_VAL(nInBufferSize);
332 41 : CPLAssert(nInBufferSize == nElts * nInBands * sizeof(double));
333 41 : CPL_IGNORE_RET_VAL(nOutBufferSize);
334 41 : CPLAssert(nOutBufferSize == nElts * nOutBands * sizeof(double));
335 :
336 41 : const BandAffineCombinationData *data =
337 : static_cast<BandAffineCombinationData *>(pWorkingData);
338 41 : CPLAssert(data->m_osSignature ==
339 : BandAffineCombinationData::EXPECTED_SIGNATURE);
340 41 : const double *CPL_RESTRICT padfSrc = static_cast<const double *>(pInBuffer);
341 41 : double *CPL_RESTRICT padfDst = static_cast<double *>(pOutBuffer);
342 : const bool bDstIntendedDTIsInteger =
343 41 : GDALDataTypeIsInteger(data->m_eIntendedDstDT);
344 41 : const double dfClampMin = data->m_dfClampMin;
345 41 : const double dfClampMax = data->m_dfClampMax;
346 1919 : for (size_t i = 0; i < nElts; ++i)
347 : {
348 5068 : for (int iDst = 0; iDst < nOutBands; ++iDst)
349 : {
350 3190 : const auto &adfCoefficients = data->m_aadfCoefficients[iDst];
351 3190 : double dfVal = adfCoefficients[0];
352 3190 : bool bSetNoData = false;
353 7940 : for (int iSrc = 0; iSrc < nInBands; ++iSrc)
354 : {
355 : // written this way to work with a NaN value
356 4762 : if (!(padfSrc[iSrc] != padfInNoData[iSrc]))
357 : {
358 12 : bSetNoData = true;
359 12 : break;
360 : }
361 4750 : dfVal += adfCoefficients[iSrc + 1] * padfSrc[iSrc];
362 : }
363 3190 : if (bSetNoData)
364 : {
365 12 : *padfDst = padfOutNoData[iDst];
366 : }
367 : else
368 : {
369 9534 : double dfDstVal = GetDstValue(
370 3178 : dfVal, padfOutNoData[iDst],
371 3178 : data->m_adfReplacementDstNodata[iDst],
372 3178 : data->m_eIntendedDstDT, bDstIntendedDTIsInteger);
373 3178 : if (dfDstVal < dfClampMin)
374 2 : dfDstVal = dfClampMin;
375 3178 : if (dfDstVal > dfClampMax)
376 2 : dfDstVal = dfClampMax;
377 3178 : *padfDst = dfDstVal;
378 : }
379 3190 : ++padfDst;
380 : }
381 1878 : padfSrc += nInBands;
382 : }
383 :
384 41 : return CE_None;
385 : }
386 :
387 : /************************************************************************/
388 : /* LUTData */
389 : /************************************************************************/
390 :
391 : namespace
392 : {
393 : /** Working structure for 'LUT' builtin function. */
394 : struct LUTData
395 : {
396 : static constexpr const char *const EXPECTED_SIGNATURE = "LUT";
397 : //! Signature (to make sure callback functions are called with the right argument)
398 : const std::string m_osSignature = EXPECTED_SIGNATURE;
399 :
400 : //! m_aadfLUTInputs[i][j] is the j(th) input value for that LUT of band i.
401 : std::vector<std::vector<double>> m_aadfLUTInputs{};
402 :
403 : //! m_aadfLUTOutputs[i][j] is the j(th) output value for that LUT of band i.
404 : std::vector<std::vector<double>> m_aadfLUTOutputs{};
405 :
406 : /************************************************************************/
407 : /* LookupValue() */
408 : /************************************************************************/
409 :
410 18 : double LookupValue(int iBand, double dfInput) const
411 : {
412 18 : const auto &adfInput = m_aadfLUTInputs[iBand];
413 18 : const auto &afdOutput = m_aadfLUTOutputs[iBand];
414 :
415 : // Find the index of the first element in the LUT input array that
416 : // is not smaller than the input value.
417 : int i = static_cast<int>(
418 18 : std::lower_bound(adfInput.data(), adfInput.data() + adfInput.size(),
419 18 : dfInput) -
420 18 : adfInput.data());
421 :
422 18 : if (i == 0)
423 6 : return afdOutput[0];
424 :
425 : // If the index is beyond the end of the LUT input array, the input
426 : // value is larger than all the values in the array.
427 12 : if (i == static_cast<int>(adfInput.size()))
428 6 : return afdOutput.back();
429 :
430 6 : if (adfInput[i] == dfInput)
431 0 : return afdOutput[i];
432 :
433 : // Otherwise, interpolate.
434 6 : return afdOutput[i - 1] + (dfInput - adfInput[i - 1]) *
435 6 : ((afdOutput[i] - afdOutput[i - 1]) /
436 6 : (adfInput[i] - adfInput[i - 1]));
437 : }
438 : };
439 : } // namespace
440 :
441 : /************************************************************************/
442 : /* LUTInit() */
443 : /************************************************************************/
444 :
445 : /** Init function for 'LUT' builtin function. */
446 6 : static CPLErr LUTInit(const char * /*pszFuncName*/, void * /*pUserData*/,
447 : CSLConstList papszFunctionArgs, int nInBands,
448 : GDALDataType eInDT, double *padfInNoData, int *pnOutBands,
449 : GDALDataType *peOutDT, double **ppadfOutNoData,
450 : const char * /* pszVRTPath */,
451 : VRTPDWorkingDataPtr *ppWorkingData)
452 : {
453 6 : CPLAssert(eInDT == GDT_Float64);
454 :
455 6 : const bool bIsFinalStep = *pnOutBands != 0;
456 6 : *peOutDT = eInDT;
457 6 : *ppWorkingData = nullptr;
458 :
459 6 : if (!bIsFinalStep)
460 : {
461 0 : *pnOutBands = nInBands;
462 : }
463 :
464 12 : auto data = std::make_unique<LUTData>();
465 :
466 6 : double dfSrcNoData = std::numeric_limits<double>::quiet_NaN();
467 6 : bool bSrcNodataSpecified = false;
468 6 : double dfDstNoData = std::numeric_limits<double>::quiet_NaN();
469 6 : bool bDstNodataSpecified = false;
470 :
471 12 : std::map<int, std::pair<std::vector<double>, std::vector<double>>> oMap{};
472 :
473 22 : for (const auto &[pszKey, pszValue] :
474 26 : cpl::IterateNameValue(papszFunctionArgs))
475 : {
476 12 : if (EQUAL(pszKey, "src_nodata"))
477 : {
478 1 : bSrcNodataSpecified = true;
479 1 : dfSrcNoData = CPLAtof(pszValue);
480 : }
481 11 : else if (EQUAL(pszKey, "dst_nodata"))
482 : {
483 1 : bDstNodataSpecified = true;
484 1 : dfDstNoData = CPLAtof(pszValue);
485 : }
486 10 : else if (STARTS_WITH_CI(pszKey, "lut_"))
487 : {
488 10 : const int nBand = atoi(pszKey + strlen("lut_"));
489 10 : if (nBand <= 0 || nBand > nInBands)
490 : {
491 1 : CPLError(CE_Failure, CPLE_AppDefined,
492 : "Invalid band in argument '%s'", pszKey);
493 2 : return CE_Failure;
494 : }
495 9 : const CPLStringList aosTokens(CSLTokenizeString2(pszValue, ",", 0));
496 9 : std::vector<double> adfInputValues;
497 9 : std::vector<double> adfOutputValues;
498 26 : for (int i = 0; i < aosTokens.size(); ++i)
499 : {
500 : const CPLStringList aosTokens2(
501 18 : CSLTokenizeString2(aosTokens[i], ":", 0));
502 18 : if (aosTokens2.size() != 2)
503 : {
504 1 : CPLError(CE_Failure, CPLE_AppDefined,
505 : "Invalid value for argument '%s'", pszKey);
506 1 : return CE_Failure;
507 : }
508 17 : adfInputValues.push_back(CPLAtof(aosTokens2[0]));
509 17 : adfOutputValues.push_back(CPLAtof(aosTokens2[1]));
510 : }
511 16 : oMap[nBand - 1] = std::pair(std::move(adfInputValues),
512 16 : std::move(adfOutputValues));
513 : }
514 : else
515 : {
516 0 : CPLError(CE_Warning, CPLE_AppDefined,
517 : "Unrecognized argument name %s. Ignored", pszKey);
518 : }
519 : }
520 :
521 4 : SetOutputValuesForInNoDataAndOutNoData(
522 : nInBands, padfInNoData, pnOutBands, ppadfOutNoData, bSrcNodataSpecified,
523 : dfSrcNoData, bDstNodataSpecified, dfDstNoData, bIsFinalStep);
524 :
525 4 : int iExpected = 0;
526 : // Check we have values for all bands and convert to vector
527 11 : for (auto &oIter : oMap)
528 : {
529 7 : if (oIter.first != iExpected)
530 : {
531 0 : CPLError(CE_Failure, CPLE_AppDefined, "Argument lut_%d is missing",
532 : iExpected + 1);
533 0 : return CE_Failure;
534 : }
535 7 : ++iExpected;
536 7 : data->m_aadfLUTInputs.emplace_back(std::move(oIter.second.first));
537 7 : data->m_aadfLUTOutputs.emplace_back(std::move(oIter.second.second));
538 : }
539 :
540 4 : if (static_cast<int>(oMap.size()) < *pnOutBands)
541 : {
542 1 : CPLError(CE_Failure, CPLE_AppDefined, "Missing lut_XX element(s)");
543 1 : return CE_Failure;
544 : }
545 :
546 3 : *ppWorkingData = data.release();
547 3 : return CE_None;
548 : }
549 :
550 : /************************************************************************/
551 : /* LUTFree() */
552 : /************************************************************************/
553 :
554 : /** Free function for 'LUT' builtin function. */
555 3 : static void LUTFree(const char * /*pszFuncName*/, void * /*pUserData*/,
556 : VRTPDWorkingDataPtr pWorkingData)
557 : {
558 3 : LUTData *data = static_cast<LUTData *>(pWorkingData);
559 3 : CPLAssert(data->m_osSignature == LUTData::EXPECTED_SIGNATURE);
560 3 : CPL_IGNORE_RET_VAL(data->m_osSignature);
561 3 : delete data;
562 3 : }
563 :
564 : /************************************************************************/
565 : /* LUTProcess() */
566 : /************************************************************************/
567 :
568 : /** Processing function for 'LUT' builtin function. */
569 : static CPLErr
570 3 : LUTProcess(const char * /*pszFuncName*/, void * /*pUserData*/,
571 : VRTPDWorkingDataPtr pWorkingData,
572 : CSLConstList /* papszFunctionArgs*/, int nBufXSize, int nBufYSize,
573 : const void *pInBuffer, size_t nInBufferSize, GDALDataType eInDT,
574 : int nInBands, const double *CPL_RESTRICT padfInNoData,
575 : void *pOutBuffer, size_t nOutBufferSize, GDALDataType eOutDT,
576 : int nOutBands, const double *CPL_RESTRICT padfOutNoData,
577 : double /*dfSrcXOff*/, double /*dfSrcYOff*/, double /*dfSrcXSize*/,
578 : double /*dfSrcYSize*/, const double /*adfSrcGT*/[],
579 : const char * /* pszVRTPath */, CSLConstList /*papszExtra*/)
580 : {
581 3 : const size_t nElts = static_cast<size_t>(nBufXSize) * nBufYSize;
582 :
583 3 : CPL_IGNORE_RET_VAL(eInDT);
584 3 : CPLAssert(eInDT == GDT_Float64);
585 3 : CPL_IGNORE_RET_VAL(eOutDT);
586 3 : CPLAssert(eOutDT == GDT_Float64);
587 3 : CPL_IGNORE_RET_VAL(nInBufferSize);
588 3 : CPLAssert(nInBufferSize == nElts * nInBands * sizeof(double));
589 3 : CPL_IGNORE_RET_VAL(nOutBufferSize);
590 3 : CPLAssert(nOutBufferSize == nElts * nOutBands * sizeof(double));
591 3 : CPLAssert(nInBands == nOutBands);
592 3 : CPL_IGNORE_RET_VAL(nOutBands);
593 :
594 3 : const LUTData *data = static_cast<LUTData *>(pWorkingData);
595 3 : CPLAssert(data->m_osSignature == LUTData::EXPECTED_SIGNATURE);
596 3 : const double *CPL_RESTRICT padfSrc = static_cast<const double *>(pInBuffer);
597 3 : double *CPL_RESTRICT padfDst = static_cast<double *>(pOutBuffer);
598 14 : for (size_t i = 0; i < nElts; ++i)
599 : {
600 33 : for (int iBand = 0; iBand < nInBands; ++iBand)
601 : {
602 : // written this way to work with a NaN value
603 22 : if (!(*padfSrc != padfInNoData[iBand]))
604 4 : *padfDst = padfOutNoData[iBand];
605 : else
606 18 : *padfDst = data->LookupValue(iBand, *padfSrc);
607 22 : ++padfSrc;
608 22 : ++padfDst;
609 : }
610 : }
611 :
612 3 : return CE_None;
613 : }
614 :
615 : /************************************************************************/
616 : /* LocalScaleOffsetData */
617 : /************************************************************************/
618 :
619 : namespace
620 : {
621 : /** Working structure for 'LocalScaleOffset' builtin function. */
622 : struct LocalScaleOffsetData
623 : {
624 : static constexpr const char *const EXPECTED_SIGNATURE = "LocalScaleOffset";
625 : //! Signature (to make sure callback functions are called with the right argument)
626 : const std::string m_osSignature = EXPECTED_SIGNATURE;
627 :
628 : //! Nodata value for gain dataset(s)
629 : double m_dfGainNodata = std::numeric_limits<double>::quiet_NaN();
630 :
631 : //! Nodata value for offset dataset(s)
632 : double m_dfOffsetNodata = std::numeric_limits<double>::quiet_NaN();
633 :
634 : //! Minimum clamping value.
635 : double m_dfClampMin = std::numeric_limits<double>::quiet_NaN();
636 :
637 : //! Maximum clamping value.
638 : double m_dfClampMax = std::numeric_limits<double>::quiet_NaN();
639 :
640 : //! Map from gain/offset dataset name to datasets
641 : std::map<std::string, std::unique_ptr<GDALDataset>> m_oDatasetMap{};
642 :
643 : //! Vector of size nInBands that point to the raster band from which to read gains.
644 : std::vector<GDALRasterBand *> m_oGainBands{};
645 :
646 : //! Vector of size nInBands that point to the raster band from which to read offsets.
647 : std::vector<GDALRasterBand *> m_oOffsetBands{};
648 :
649 : //! Working buffer that contain gain values.
650 : std::vector<VRTProcessedDataset::NoInitByte> m_abyGainBuffer{};
651 :
652 : //! Working buffer that contain offset values.
653 : std::vector<VRTProcessedDataset::NoInitByte> m_abyOffsetBuffer{};
654 : };
655 : } // namespace
656 :
657 : /************************************************************************/
658 : /* CheckAllBands() */
659 : /************************************************************************/
660 :
661 : /** Return true if the key of oMap is the sequence of all integers between
662 : * 0 and nExpectedBandCount-1.
663 : */
664 : template <class T>
665 28 : static bool CheckAllBands(const std::map<int, T> &oMap, int nExpectedBandCount)
666 : {
667 28 : int iExpected = 0;
668 60 : for (const auto &kv : oMap)
669 : {
670 32 : if (kv.first != iExpected)
671 0 : return false;
672 32 : ++iExpected;
673 : }
674 28 : return iExpected == nExpectedBandCount;
675 : }
676 :
677 : /************************************************************************/
678 : /* LocalScaleOffsetInit() */
679 : /************************************************************************/
680 :
681 : /** Init function for 'LocalScaleOffset' builtin function. */
682 : static CPLErr
683 11 : LocalScaleOffsetInit(const char * /*pszFuncName*/, void * /*pUserData*/,
684 : CSLConstList papszFunctionArgs, int nInBands,
685 : GDALDataType eInDT, double *padfInNoData, int *pnOutBands,
686 : GDALDataType *peOutDT, double **ppadfOutNoData,
687 : const char *pszVRTPath, VRTPDWorkingDataPtr *ppWorkingData)
688 : {
689 11 : CPLAssert(eInDT == GDT_Float64);
690 :
691 11 : const bool bIsFinalStep = *pnOutBands != 0;
692 11 : *peOutDT = eInDT;
693 11 : *ppWorkingData = nullptr;
694 :
695 11 : if (!bIsFinalStep)
696 : {
697 0 : *pnOutBands = nInBands;
698 : }
699 :
700 22 : auto data = std::make_unique<LocalScaleOffsetData>();
701 :
702 11 : bool bNodataSpecified = false;
703 11 : double dfNoData = std::numeric_limits<double>::quiet_NaN();
704 :
705 11 : bool bGainNodataSpecified = false;
706 11 : bool bOffsetNodataSpecified = false;
707 :
708 22 : std::map<int, std::string> oGainDatasetNameMap;
709 22 : std::map<int, int> oGainDatasetBandMap;
710 :
711 22 : std::map<int, std::string> oOffsetDatasetNameMap;
712 22 : std::map<int, int> oOffsetDatasetBandMap;
713 :
714 11 : bool bRelativeToVRT = false;
715 :
716 84 : for (const auto &[pszKey, pszValue] :
717 91 : cpl::IterateNameValue(papszFunctionArgs))
718 : {
719 44 : if (EQUAL(pszKey, "relativeToVRT"))
720 : {
721 0 : bRelativeToVRT = CPLTestBool(pszValue);
722 : }
723 44 : else if (EQUAL(pszKey, "nodata"))
724 : {
725 0 : bNodataSpecified = true;
726 0 : dfNoData = CPLAtof(pszValue);
727 : }
728 44 : else if (EQUAL(pszKey, "gain_nodata"))
729 : {
730 0 : bGainNodataSpecified = true;
731 0 : data->m_dfGainNodata = CPLAtof(pszValue);
732 : }
733 44 : else if (EQUAL(pszKey, "offset_nodata"))
734 : {
735 0 : bOffsetNodataSpecified = true;
736 0 : data->m_dfOffsetNodata = CPLAtof(pszValue);
737 : }
738 44 : else if (STARTS_WITH_CI(pszKey, "gain_dataset_filename_"))
739 : {
740 12 : const int nBand = atoi(pszKey + strlen("gain_dataset_filename_"));
741 12 : if (nBand <= 0 || nBand > nInBands)
742 : {
743 1 : CPLError(CE_Failure, CPLE_AppDefined,
744 : "Invalid band in argument '%s'", pszKey);
745 4 : return CE_Failure;
746 : }
747 11 : oGainDatasetNameMap[nBand - 1] = pszValue;
748 : }
749 32 : else if (STARTS_WITH_CI(pszKey, "gain_dataset_band_"))
750 : {
751 11 : const int nBand = atoi(pszKey + strlen("gain_dataset_band_"));
752 11 : if (nBand <= 0 || nBand > nInBands)
753 : {
754 1 : CPLError(CE_Failure, CPLE_AppDefined,
755 : "Invalid band in argument '%s'", pszKey);
756 1 : return CE_Failure;
757 : }
758 10 : oGainDatasetBandMap[nBand - 1] = atoi(pszValue);
759 : }
760 21 : else if (STARTS_WITH_CI(pszKey, "offset_dataset_filename_"))
761 : {
762 10 : const int nBand = atoi(pszKey + strlen("offset_dataset_filename_"));
763 10 : if (nBand <= 0 || nBand > nInBands)
764 : {
765 1 : CPLError(CE_Failure, CPLE_AppDefined,
766 : "Invalid band in argument '%s'", pszKey);
767 1 : return CE_Failure;
768 : }
769 9 : oOffsetDatasetNameMap[nBand - 1] = pszValue;
770 : }
771 11 : else if (STARTS_WITH_CI(pszKey, "offset_dataset_band_"))
772 : {
773 9 : const int nBand = atoi(pszKey + strlen("offset_dataset_band_"));
774 9 : if (nBand <= 0 || nBand > nInBands)
775 : {
776 1 : CPLError(CE_Failure, CPLE_AppDefined,
777 : "Invalid band in argument '%s'", pszKey);
778 1 : return CE_Failure;
779 : }
780 8 : oOffsetDatasetBandMap[nBand - 1] = atoi(pszValue);
781 : }
782 2 : else if (EQUAL(pszKey, "min"))
783 : {
784 1 : data->m_dfClampMin = CPLAtof(pszValue);
785 : }
786 1 : else if (EQUAL(pszKey, "max"))
787 : {
788 1 : data->m_dfClampMax = CPLAtof(pszValue);
789 : }
790 : else
791 : {
792 0 : CPLError(CE_Warning, CPLE_AppDefined,
793 : "Unrecognized argument name %s. Ignored", pszKey);
794 : }
795 : }
796 :
797 7 : if (!CheckAllBands(oGainDatasetNameMap, nInBands))
798 : {
799 0 : CPLError(CE_Failure, CPLE_AppDefined,
800 : "Missing gain_dataset_filename_XX element(s)");
801 0 : return CE_Failure;
802 : }
803 7 : if (!CheckAllBands(oGainDatasetBandMap, nInBands))
804 : {
805 0 : CPLError(CE_Failure, CPLE_AppDefined,
806 : "Missing gain_dataset_band_XX element(s)");
807 0 : return CE_Failure;
808 : }
809 7 : if (!CheckAllBands(oOffsetDatasetNameMap, nInBands))
810 : {
811 0 : CPLError(CE_Failure, CPLE_AppDefined,
812 : "Missing offset_dataset_filename_XX element(s)");
813 0 : return CE_Failure;
814 : }
815 7 : if (!CheckAllBands(oOffsetDatasetBandMap, nInBands))
816 : {
817 0 : CPLError(CE_Failure, CPLE_AppDefined,
818 : "Missing offset_dataset_band_XX element(s)");
819 0 : return CE_Failure;
820 : }
821 :
822 7 : data->m_oGainBands.resize(nInBands);
823 7 : data->m_oOffsetBands.resize(nInBands);
824 :
825 7 : constexpr int IDX_GAIN = 0;
826 7 : constexpr int IDX_OFFSET = 1;
827 15 : for (int i : {IDX_GAIN, IDX_OFFSET})
828 : {
829 11 : const auto &oMapNames =
830 : (i == IDX_GAIN) ? oGainDatasetNameMap : oOffsetDatasetNameMap;
831 11 : const auto &oMapBands =
832 : (i == IDX_GAIN) ? oGainDatasetBandMap : oOffsetDatasetBandMap;
833 21 : for (const auto &kv : oMapNames)
834 : {
835 13 : const int nInBandIdx = kv.first;
836 : const auto osFilename = GDALDataset::BuildFilename(
837 13 : kv.second.c_str(), pszVRTPath, bRelativeToVRT);
838 13 : auto oIter = data->m_oDatasetMap.find(osFilename);
839 13 : if (oIter == data->m_oDatasetMap.end())
840 : {
841 : auto poDS = std::unique_ptr<GDALDataset>(GDALDataset::Open(
842 : osFilename.c_str(), GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR,
843 11 : nullptr, nullptr, nullptr));
844 11 : if (!poDS)
845 1 : return CE_Failure;
846 : double adfAuxGT[6];
847 10 : if (poDS->GetGeoTransform(adfAuxGT) != CE_None)
848 : {
849 1 : CPLError(CE_Failure, CPLE_AppDefined,
850 : "%s lacks a geotransform", osFilename.c_str());
851 1 : return CE_Failure;
852 : }
853 9 : oIter = data->m_oDatasetMap
854 9 : .insert(std::pair(osFilename, std::move(poDS)))
855 : .first;
856 : }
857 11 : auto poDS = oIter->second.get();
858 11 : const auto oIterBand = oMapBands.find(nInBandIdx);
859 11 : CPLAssert(oIterBand != oMapBands.end());
860 11 : const int nAuxBand = oIterBand->second;
861 11 : if (nAuxBand <= 0 || nAuxBand > poDS->GetRasterCount())
862 : {
863 1 : CPLError(CE_Failure, CPLE_AppDefined,
864 : "Invalid band number (%d) for a %s dataset", nAuxBand,
865 : (i == IDX_GAIN) ? "gain" : "offset");
866 1 : return CE_Failure;
867 : }
868 10 : auto poAuxBand = poDS->GetRasterBand(nAuxBand);
869 10 : int bAuxBandHasNoData = false;
870 : const double dfAuxNoData =
871 10 : poAuxBand->GetNoDataValue(&bAuxBandHasNoData);
872 10 : if (i == IDX_GAIN)
873 : {
874 5 : data->m_oGainBands[nInBandIdx] = poAuxBand;
875 5 : if (!bGainNodataSpecified && bAuxBandHasNoData)
876 2 : data->m_dfGainNodata = dfAuxNoData;
877 : }
878 : else
879 : {
880 5 : data->m_oOffsetBands[nInBandIdx] = poAuxBand;
881 5 : if (!bOffsetNodataSpecified && bAuxBandHasNoData)
882 2 : data->m_dfOffsetNodata = dfAuxNoData;
883 : }
884 : }
885 : }
886 :
887 4 : SetOutputValuesForInNoDataAndOutNoData(
888 : nInBands, padfInNoData, pnOutBands, ppadfOutNoData, bNodataSpecified,
889 : dfNoData, bNodataSpecified, dfNoData, bIsFinalStep);
890 :
891 4 : *ppWorkingData = data.release();
892 4 : return CE_None;
893 : }
894 :
895 : /************************************************************************/
896 : /* LocalScaleOffsetFree() */
897 : /************************************************************************/
898 :
899 : /** Free function for 'LocalScaleOffset' builtin function. */
900 4 : static void LocalScaleOffsetFree(const char * /*pszFuncName*/,
901 : void * /*pUserData*/,
902 : VRTPDWorkingDataPtr pWorkingData)
903 : {
904 4 : LocalScaleOffsetData *data =
905 : static_cast<LocalScaleOffsetData *>(pWorkingData);
906 4 : CPLAssert(data->m_osSignature == LocalScaleOffsetData::EXPECTED_SIGNATURE);
907 4 : CPL_IGNORE_RET_VAL(data->m_osSignature);
908 4 : delete data;
909 4 : }
910 :
911 : /************************************************************************/
912 : /* LoadAuxData() */
913 : /************************************************************************/
914 :
915 : // Load auxiliary corresponding offset, gain or trimming data.
916 17 : static bool LoadAuxData(double dfULX, double dfULY, double dfLRX, double dfLRY,
917 : size_t nElts, int nBufXSize, int nBufYSize,
918 : const char *pszAuxType, GDALRasterBand *poAuxBand,
919 : std::vector<VRTProcessedDataset::NoInitByte> &abyBuffer)
920 : {
921 : double adfAuxGT[6];
922 : double adfAuxInvGT[6];
923 :
924 : // Compute pixel/line coordinates from the georeferenced extent
925 34 : CPL_IGNORE_RET_VAL(poAuxBand->GetDataset()->GetGeoTransform(
926 17 : adfAuxGT)); // return code already tested
927 17 : CPL_IGNORE_RET_VAL(GDALInvGeoTransform(adfAuxGT, adfAuxInvGT));
928 17 : const double dfULPixel =
929 17 : adfAuxInvGT[0] + adfAuxInvGT[1] * dfULX + adfAuxInvGT[2] * dfULY;
930 17 : const double dfULLine =
931 17 : adfAuxInvGT[3] + adfAuxInvGT[4] * dfULX + adfAuxInvGT[5] * dfULY;
932 17 : const double dfLRPixel =
933 17 : adfAuxInvGT[0] + adfAuxInvGT[1] * dfLRX + adfAuxInvGT[2] * dfLRY;
934 17 : const double dfLRLine =
935 17 : adfAuxInvGT[3] + adfAuxInvGT[4] * dfLRX + adfAuxInvGT[5] * dfLRY;
936 17 : if (dfULPixel >= dfLRPixel || dfULLine >= dfLRLine)
937 : {
938 0 : CPLError(CE_Failure, CPLE_AppDefined,
939 : "Unexpected computed %s pixel/line", pszAuxType);
940 0 : return false;
941 : }
942 17 : if (dfULPixel < -1 || dfULLine < -1)
943 : {
944 0 : CPLError(CE_Failure, CPLE_AppDefined,
945 : "Unexpected computed %s upper left (pixel,line)=(%f,%f)",
946 : pszAuxType, dfULPixel, dfULLine);
947 0 : return false;
948 : }
949 34 : if (dfLRPixel > poAuxBand->GetXSize() + 1 ||
950 17 : dfLRLine > poAuxBand->GetYSize() + 1)
951 : {
952 0 : CPLError(CE_Failure, CPLE_AppDefined,
953 : "Unexpected computed %s lower right (pixel,line)=(%f,%f)",
954 : pszAuxType, dfLRPixel, dfLRLine);
955 0 : return false;
956 : }
957 :
958 34 : const int nAuxXOff = std::clamp(static_cast<int>(std::round(dfULPixel)), 0,
959 17 : poAuxBand->GetXSize() - 1);
960 34 : const int nAuxYOff = std::clamp(static_cast<int>(std::round(dfULLine)), 0,
961 17 : poAuxBand->GetYSize() - 1);
962 51 : const int nAuxX2Off = std::min(poAuxBand->GetXSize(),
963 17 : static_cast<int>(std::round(dfLRPixel)));
964 : const int nAuxY2Off =
965 17 : std::min(poAuxBand->GetYSize(), static_cast<int>(std::round(dfLRLine)));
966 :
967 : try
968 : {
969 17 : abyBuffer.resize(nElts * sizeof(float));
970 : }
971 0 : catch (const std::bad_alloc &)
972 : {
973 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
974 : "Out of memory allocating working buffer");
975 0 : return false;
976 : }
977 : GDALRasterIOExtraArg sExtraArg;
978 17 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
979 17 : sExtraArg.bFloatingPointWindowValidity = true;
980 17 : CPL_IGNORE_RET_VAL(sExtraArg.eResampleAlg);
981 17 : sExtraArg.eResampleAlg = GRIORA_Bilinear;
982 17 : sExtraArg.dfXOff = std::max(0.0, dfULPixel);
983 17 : sExtraArg.dfYOff = std::max(0.0, dfULLine);
984 17 : sExtraArg.dfXSize = std::min<double>(poAuxBand->GetXSize(), dfLRPixel) -
985 17 : std::max(0.0, dfULPixel);
986 17 : sExtraArg.dfYSize = std::min<double>(poAuxBand->GetYSize(), dfLRLine) -
987 17 : std::max(0.0, dfULLine);
988 17 : return (poAuxBand->RasterIO(
989 17 : GF_Read, nAuxXOff, nAuxYOff, std::max(1, nAuxX2Off - nAuxXOff),
990 17 : std::max(1, nAuxY2Off - nAuxYOff), abyBuffer.data(), nBufXSize,
991 17 : nBufYSize, GDT_Float32, 0, 0, &sExtraArg) == CE_None);
992 : }
993 :
994 : /************************************************************************/
995 : /* LocalScaleOffsetProcess() */
996 : /************************************************************************/
997 :
998 : /** Processing function for 'LocalScaleOffset' builtin function. */
999 7 : static CPLErr LocalScaleOffsetProcess(
1000 : const char * /*pszFuncName*/, void * /*pUserData*/,
1001 : VRTPDWorkingDataPtr pWorkingData, CSLConstList /* papszFunctionArgs*/,
1002 : int nBufXSize, int nBufYSize, const void *pInBuffer, size_t nInBufferSize,
1003 : GDALDataType eInDT, int nInBands, const double *CPL_RESTRICT padfInNoData,
1004 : void *pOutBuffer, size_t nOutBufferSize, GDALDataType eOutDT, int nOutBands,
1005 : const double *CPL_RESTRICT padfOutNoData, double dfSrcXOff,
1006 : double dfSrcYOff, double dfSrcXSize, double dfSrcYSize,
1007 : const double adfSrcGT[], const char * /* pszVRTPath */,
1008 : CSLConstList /*papszExtra*/)
1009 : {
1010 7 : const size_t nElts = static_cast<size_t>(nBufXSize) * nBufYSize;
1011 :
1012 7 : CPL_IGNORE_RET_VAL(eInDT);
1013 7 : CPLAssert(eInDT == GDT_Float64);
1014 7 : CPL_IGNORE_RET_VAL(eOutDT);
1015 7 : CPLAssert(eOutDT == GDT_Float64);
1016 7 : CPL_IGNORE_RET_VAL(nInBufferSize);
1017 7 : CPLAssert(nInBufferSize == nElts * nInBands * sizeof(double));
1018 7 : CPL_IGNORE_RET_VAL(nOutBufferSize);
1019 7 : CPLAssert(nOutBufferSize == nElts * nOutBands * sizeof(double));
1020 7 : CPLAssert(nInBands == nOutBands);
1021 7 : CPL_IGNORE_RET_VAL(nOutBands);
1022 :
1023 7 : LocalScaleOffsetData *data =
1024 : static_cast<LocalScaleOffsetData *>(pWorkingData);
1025 7 : CPLAssert(data->m_osSignature == LocalScaleOffsetData::EXPECTED_SIGNATURE);
1026 7 : const double *CPL_RESTRICT padfSrc = static_cast<const double *>(pInBuffer);
1027 7 : double *CPL_RESTRICT padfDst = static_cast<double *>(pOutBuffer);
1028 :
1029 : // Compute georeferenced extent of input region
1030 7 : const double dfULX =
1031 7 : adfSrcGT[0] + adfSrcGT[1] * dfSrcXOff + adfSrcGT[2] * dfSrcYOff;
1032 7 : const double dfULY =
1033 7 : adfSrcGT[3] + adfSrcGT[4] * dfSrcXOff + adfSrcGT[5] * dfSrcYOff;
1034 7 : const double dfLRX = adfSrcGT[0] + adfSrcGT[1] * (dfSrcXOff + dfSrcXSize) +
1035 7 : adfSrcGT[2] * (dfSrcYOff + dfSrcYSize);
1036 7 : const double dfLRY = adfSrcGT[3] + adfSrcGT[4] * (dfSrcXOff + dfSrcXSize) +
1037 7 : adfSrcGT[5] * (dfSrcYOff + dfSrcYSize);
1038 :
1039 7 : auto &abyOffsetBuffer = data->m_abyGainBuffer;
1040 7 : auto &abyGainBuffer = data->m_abyOffsetBuffer;
1041 :
1042 15 : for (int iBand = 0; iBand < nInBands; ++iBand)
1043 : {
1044 8 : if (!LoadAuxData(dfULX, dfULY, dfLRX, dfLRY, nElts, nBufXSize,
1045 8 : nBufYSize, "gain", data->m_oGainBands[iBand],
1046 16 : abyGainBuffer) ||
1047 8 : !LoadAuxData(dfULX, dfULY, dfLRX, dfLRY, nElts, nBufXSize,
1048 8 : nBufYSize, "offset", data->m_oOffsetBands[iBand],
1049 : abyOffsetBuffer))
1050 : {
1051 0 : return CE_Failure;
1052 : }
1053 :
1054 8 : const double *CPL_RESTRICT padfSrcThisBand = padfSrc + iBand;
1055 8 : double *CPL_RESTRICT padfDstThisBand = padfDst + iBand;
1056 : const float *pafGain =
1057 8 : reinterpret_cast<const float *>(abyGainBuffer.data());
1058 : const float *pafOffset =
1059 8 : reinterpret_cast<const float *>(abyOffsetBuffer.data());
1060 8 : const double dfSrcNodata = padfInNoData[iBand];
1061 8 : const double dfDstNodata = padfOutNoData[iBand];
1062 8 : const double dfGainNodata = data->m_dfGainNodata;
1063 8 : const double dfOffsetNodata = data->m_dfOffsetNodata;
1064 8 : const double dfClampMin = data->m_dfClampMin;
1065 8 : const double dfClampMax = data->m_dfClampMax;
1066 66084 : for (size_t i = 0; i < nElts; ++i)
1067 : {
1068 66076 : const double dfSrcVal = *padfSrcThisBand;
1069 : // written this way to work with a NaN value
1070 66076 : if (!(dfSrcVal != dfSrcNodata))
1071 : {
1072 2 : *padfDstThisBand = dfDstNodata;
1073 : }
1074 : else
1075 : {
1076 66074 : const double dfGain = pafGain[i];
1077 66074 : const double dfOffset = pafOffset[i];
1078 66074 : if (!(dfGain != dfGainNodata) || !(dfOffset != dfOffsetNodata))
1079 : {
1080 4 : *padfDstThisBand = dfDstNodata;
1081 : }
1082 : else
1083 : {
1084 66070 : double dfUnscaled = dfSrcVal * dfGain - dfOffset;
1085 66070 : if (dfUnscaled < dfClampMin)
1086 2 : dfUnscaled = dfClampMin;
1087 66070 : if (dfUnscaled > dfClampMax)
1088 1 : dfUnscaled = dfClampMax;
1089 :
1090 66070 : *padfDstThisBand = dfUnscaled;
1091 : }
1092 : }
1093 66076 : padfSrcThisBand += nInBands;
1094 66076 : padfDstThisBand += nInBands;
1095 : }
1096 : }
1097 :
1098 7 : return CE_None;
1099 : }
1100 :
1101 : /************************************************************************/
1102 : /* TrimmingData */
1103 : /************************************************************************/
1104 :
1105 : namespace
1106 : {
1107 : /** Working structure for 'Trimming' builtin function. */
1108 : struct TrimmingData
1109 : {
1110 : static constexpr const char *const EXPECTED_SIGNATURE = "Trimming";
1111 : //! Signature (to make sure callback functions are called with the right argument)
1112 : const std::string m_osSignature = EXPECTED_SIGNATURE;
1113 :
1114 : //! Nodata value for trimming dataset
1115 : double m_dfTrimmingNodata = std::numeric_limits<double>::quiet_NaN();
1116 :
1117 : //! Maximum saturating RGB output value.
1118 : double m_dfTopRGB = 0;
1119 :
1120 : //! Maximum threshold beyond which we give up saturation
1121 : double m_dfToneCeil = 0;
1122 :
1123 : //! Margin to allow for dynamics in brighest areas (in [0,1] range)
1124 : double m_dfTopMargin = 0;
1125 :
1126 : //! Index (zero-based) of input/output red band.
1127 : int m_nRedBand = 1 - 1;
1128 :
1129 : //! Index (zero-based) of input/output green band.
1130 : int m_nGreenBand = 2 - 1;
1131 :
1132 : //! Index (zero-based) of input/output blue band.
1133 : int m_nBlueBand = 3 - 1;
1134 :
1135 : //! Trimming dataset
1136 : std::unique_ptr<GDALDataset> m_poTrimmingDS{};
1137 :
1138 : //! Trimming raster band.
1139 : GDALRasterBand *m_poTrimmingBand = nullptr;
1140 :
1141 : //! Working buffer that contain trimming values.
1142 : std::vector<VRTProcessedDataset::NoInitByte> m_abyTrimmingBuffer{};
1143 : };
1144 : } // namespace
1145 :
1146 : /************************************************************************/
1147 : /* TrimmingInit() */
1148 : /************************************************************************/
1149 :
1150 : /** Init function for 'Trimming' builtin function. */
1151 12 : static CPLErr TrimmingInit(const char * /*pszFuncName*/, void * /*pUserData*/,
1152 : CSLConstList papszFunctionArgs, int nInBands,
1153 : GDALDataType eInDT, double *padfInNoData,
1154 : int *pnOutBands, GDALDataType *peOutDT,
1155 : double **ppadfOutNoData, const char *pszVRTPath,
1156 : VRTPDWorkingDataPtr *ppWorkingData)
1157 : {
1158 12 : CPLAssert(eInDT == GDT_Float64);
1159 :
1160 12 : const bool bIsFinalStep = *pnOutBands != 0;
1161 12 : *peOutDT = eInDT;
1162 12 : *ppWorkingData = nullptr;
1163 :
1164 12 : if (!bIsFinalStep)
1165 : {
1166 0 : *pnOutBands = nInBands;
1167 : }
1168 :
1169 24 : auto data = std::make_unique<TrimmingData>();
1170 :
1171 12 : bool bNodataSpecified = false;
1172 12 : double dfNoData = std::numeric_limits<double>::quiet_NaN();
1173 24 : std::string osTrimmingFilename;
1174 12 : bool bTrimmingNodataSpecified = false;
1175 12 : bool bRelativeToVRT = false;
1176 :
1177 84 : for (const auto &[pszKey, pszValue] :
1178 90 : cpl::IterateNameValue(papszFunctionArgs))
1179 : {
1180 45 : if (EQUAL(pszKey, "relativeToVRT"))
1181 : {
1182 0 : bRelativeToVRT = CPLTestBool(pszValue);
1183 : }
1184 45 : else if (EQUAL(pszKey, "nodata"))
1185 : {
1186 0 : bNodataSpecified = true;
1187 0 : dfNoData = CPLAtof(pszValue);
1188 : }
1189 45 : else if (EQUAL(pszKey, "trimming_nodata"))
1190 : {
1191 0 : bTrimmingNodataSpecified = true;
1192 0 : data->m_dfTrimmingNodata = CPLAtof(pszValue);
1193 : }
1194 45 : else if (EQUAL(pszKey, "trimming_dataset_filename"))
1195 : {
1196 12 : osTrimmingFilename = pszValue;
1197 : }
1198 33 : else if (EQUAL(pszKey, "red_band"))
1199 : {
1200 5 : const int nBand = atoi(pszValue) - 1;
1201 5 : if (nBand < 0 || nBand >= nInBands)
1202 : {
1203 2 : CPLError(CE_Failure, CPLE_AppDefined,
1204 : "Invalid band in argument '%s'", pszKey);
1205 6 : return CE_Failure;
1206 : }
1207 3 : data->m_nRedBand = nBand;
1208 : }
1209 28 : else if (EQUAL(pszKey, "green_band"))
1210 : {
1211 5 : const int nBand = atoi(pszValue) - 1;
1212 5 : if (nBand < 0 || nBand >= nInBands)
1213 : {
1214 2 : CPLError(CE_Failure, CPLE_AppDefined,
1215 : "Invalid band in argument '%s'", pszKey);
1216 2 : return CE_Failure;
1217 : }
1218 3 : data->m_nGreenBand = nBand;
1219 : }
1220 23 : else if (EQUAL(pszKey, "blue_band"))
1221 : {
1222 5 : const int nBand = atoi(pszValue) - 1;
1223 5 : if (nBand < 0 || nBand >= nInBands)
1224 : {
1225 2 : CPLError(CE_Failure, CPLE_AppDefined,
1226 : "Invalid band in argument '%s'", pszKey);
1227 2 : return CE_Failure;
1228 : }
1229 3 : data->m_nBlueBand = nBand;
1230 : }
1231 18 : else if (EQUAL(pszKey, "top_rgb"))
1232 : {
1233 6 : data->m_dfTopRGB = CPLAtof(pszValue);
1234 : }
1235 12 : else if (EQUAL(pszKey, "tone_ceil"))
1236 : {
1237 6 : data->m_dfToneCeil = CPLAtof(pszValue);
1238 : }
1239 6 : else if (EQUAL(pszKey, "top_margin"))
1240 : {
1241 6 : data->m_dfTopMargin = CPLAtof(pszValue);
1242 : }
1243 : else
1244 : {
1245 0 : CPLError(CE_Warning, CPLE_AppDefined,
1246 : "Unrecognized argument name %s. Ignored", pszKey);
1247 : }
1248 : }
1249 :
1250 6 : if (data->m_nRedBand == data->m_nGreenBand ||
1251 10 : data->m_nRedBand == data->m_nBlueBand ||
1252 4 : data->m_nGreenBand == data->m_nBlueBand)
1253 : {
1254 3 : CPLError(
1255 : CE_Failure, CPLE_NotSupported,
1256 : "red_band, green_band and blue_band must have distinct values");
1257 3 : return CE_Failure;
1258 : }
1259 :
1260 : const auto osFilename = GDALDataset::BuildFilename(
1261 6 : osTrimmingFilename.c_str(), pszVRTPath, bRelativeToVRT);
1262 3 : data->m_poTrimmingDS.reset(GDALDataset::Open(
1263 : osFilename.c_str(), GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, nullptr,
1264 : nullptr, nullptr));
1265 3 : if (!data->m_poTrimmingDS)
1266 1 : return CE_Failure;
1267 2 : if (data->m_poTrimmingDS->GetRasterCount() != 1)
1268 : {
1269 1 : CPLError(CE_Failure, CPLE_NotSupported,
1270 : "Trimming dataset should have a single band");
1271 1 : return CE_Failure;
1272 : }
1273 1 : data->m_poTrimmingBand = data->m_poTrimmingDS->GetRasterBand(1);
1274 :
1275 : double adfAuxGT[6];
1276 1 : if (data->m_poTrimmingDS->GetGeoTransform(adfAuxGT) != CE_None)
1277 : {
1278 0 : CPLError(CE_Failure, CPLE_AppDefined, "%s lacks a geotransform",
1279 : osFilename.c_str());
1280 0 : return CE_Failure;
1281 : }
1282 1 : int bAuxBandHasNoData = false;
1283 : const double dfAuxNoData =
1284 1 : data->m_poTrimmingBand->GetNoDataValue(&bAuxBandHasNoData);
1285 1 : if (!bTrimmingNodataSpecified && bAuxBandHasNoData)
1286 0 : data->m_dfTrimmingNodata = dfAuxNoData;
1287 :
1288 1 : SetOutputValuesForInNoDataAndOutNoData(
1289 : nInBands, padfInNoData, pnOutBands, ppadfOutNoData, bNodataSpecified,
1290 : dfNoData, bNodataSpecified, dfNoData, bIsFinalStep);
1291 :
1292 1 : *ppWorkingData = data.release();
1293 1 : return CE_None;
1294 : }
1295 :
1296 : /************************************************************************/
1297 : /* TrimmingFree() */
1298 : /************************************************************************/
1299 :
1300 : /** Free function for 'Trimming' builtin function. */
1301 1 : static void TrimmingFree(const char * /*pszFuncName*/, void * /*pUserData*/,
1302 : VRTPDWorkingDataPtr pWorkingData)
1303 : {
1304 1 : TrimmingData *data = static_cast<TrimmingData *>(pWorkingData);
1305 1 : CPLAssert(data->m_osSignature == TrimmingData::EXPECTED_SIGNATURE);
1306 1 : CPL_IGNORE_RET_VAL(data->m_osSignature);
1307 1 : delete data;
1308 1 : }
1309 :
1310 : /************************************************************************/
1311 : /* TrimmingProcess() */
1312 : /************************************************************************/
1313 :
1314 : /** Processing function for 'Trimming' builtin function. */
1315 1 : static CPLErr TrimmingProcess(
1316 : const char * /*pszFuncName*/, void * /*pUserData*/,
1317 : VRTPDWorkingDataPtr pWorkingData, CSLConstList /* papszFunctionArgs*/,
1318 : int nBufXSize, int nBufYSize, const void *pInBuffer, size_t nInBufferSize,
1319 : GDALDataType eInDT, int nInBands, const double *CPL_RESTRICT padfInNoData,
1320 : void *pOutBuffer, size_t nOutBufferSize, GDALDataType eOutDT, int nOutBands,
1321 : const double *CPL_RESTRICT padfOutNoData, double dfSrcXOff,
1322 : double dfSrcYOff, double dfSrcXSize, double dfSrcYSize,
1323 : const double adfSrcGT[], const char * /* pszVRTPath */,
1324 : CSLConstList /*papszExtra*/)
1325 : {
1326 1 : const size_t nElts = static_cast<size_t>(nBufXSize) * nBufYSize;
1327 :
1328 1 : CPL_IGNORE_RET_VAL(eInDT);
1329 1 : CPLAssert(eInDT == GDT_Float64);
1330 1 : CPL_IGNORE_RET_VAL(eOutDT);
1331 1 : CPLAssert(eOutDT == GDT_Float64);
1332 1 : CPL_IGNORE_RET_VAL(nInBufferSize);
1333 1 : CPLAssert(nInBufferSize == nElts * nInBands * sizeof(double));
1334 1 : CPL_IGNORE_RET_VAL(nOutBufferSize);
1335 1 : CPLAssert(nOutBufferSize == nElts * nOutBands * sizeof(double));
1336 1 : CPLAssert(nInBands == nOutBands);
1337 1 : CPL_IGNORE_RET_VAL(nOutBands);
1338 :
1339 1 : TrimmingData *data = static_cast<TrimmingData *>(pWorkingData);
1340 1 : CPLAssert(data->m_osSignature == TrimmingData::EXPECTED_SIGNATURE);
1341 1 : const double *CPL_RESTRICT padfSrc = static_cast<const double *>(pInBuffer);
1342 1 : double *CPL_RESTRICT padfDst = static_cast<double *>(pOutBuffer);
1343 :
1344 : // Compute georeferenced extent of input region
1345 1 : const double dfULX =
1346 1 : adfSrcGT[0] + adfSrcGT[1] * dfSrcXOff + adfSrcGT[2] * dfSrcYOff;
1347 1 : const double dfULY =
1348 1 : adfSrcGT[3] + adfSrcGT[4] * dfSrcXOff + adfSrcGT[5] * dfSrcYOff;
1349 1 : const double dfLRX = adfSrcGT[0] + adfSrcGT[1] * (dfSrcXOff + dfSrcXSize) +
1350 1 : adfSrcGT[2] * (dfSrcYOff + dfSrcYSize);
1351 1 : const double dfLRY = adfSrcGT[3] + adfSrcGT[4] * (dfSrcXOff + dfSrcXSize) +
1352 1 : adfSrcGT[5] * (dfSrcYOff + dfSrcYSize);
1353 :
1354 1 : if (!LoadAuxData(dfULX, dfULY, dfLRX, dfLRY, nElts, nBufXSize, nBufYSize,
1355 : "trimming", data->m_poTrimmingBand,
1356 1 : data->m_abyTrimmingBuffer))
1357 : {
1358 0 : return CE_Failure;
1359 : }
1360 :
1361 : const float *pafTrimming =
1362 1 : reinterpret_cast<const float *>(data->m_abyTrimmingBuffer.data());
1363 1 : const int nRedBand = data->m_nRedBand;
1364 1 : const int nGreenBand = data->m_nGreenBand;
1365 1 : const int nBlueBand = data->m_nBlueBand;
1366 1 : const double dfTopMargin = data->m_dfTopMargin;
1367 1 : const double dfTopRGB = data->m_dfTopRGB;
1368 1 : const double dfToneCeil = data->m_dfToneCeil;
1369 : #if !defined(trimming_non_optimized_version)
1370 1 : const double dfInvToneCeil = 1.0 / dfToneCeil;
1371 : #endif
1372 : const bool bRGBBandsAreFirst =
1373 1 : std::max(std::max(nRedBand, nGreenBand), nBlueBand) <= 2;
1374 1 : const double dfNoDataTrimming = data->m_dfTrimmingNodata;
1375 1 : const double dfNoDataRed = padfInNoData[nRedBand];
1376 1 : const double dfNoDataGreen = padfInNoData[nGreenBand];
1377 1 : const double dfNoDataBlue = padfInNoData[nBlueBand];
1378 7 : for (size_t i = 0; i < nElts; ++i)
1379 : {
1380 : // Extract local saturation value from trimming image
1381 6 : const double dfLocalMaxRGB = pafTrimming[i];
1382 : const double dfReducedRGB =
1383 6 : std::min((1.0 - dfTopMargin) * dfTopRGB / dfLocalMaxRGB, 1.0);
1384 :
1385 6 : const double dfRed = padfSrc[nRedBand];
1386 6 : const double dfGreen = padfSrc[nGreenBand];
1387 6 : const double dfBlue = padfSrc[nBlueBand];
1388 6 : bool bNoDataPixel = false;
1389 6 : if ((dfLocalMaxRGB != dfNoDataTrimming) && (dfRed != dfNoDataRed) &&
1390 6 : (dfGreen != dfNoDataGreen) && (dfBlue != dfNoDataBlue))
1391 : {
1392 : // RGB bands specific process
1393 6 : const double dfMaxRGB = std::max(std::max(dfRed, dfGreen), dfBlue);
1394 : #if !defined(trimming_non_optimized_version)
1395 6 : const double dfRedTimesToneRed = std::min(dfRed, dfToneCeil);
1396 6 : const double dfGreenTimesToneGreen = std::min(dfGreen, dfToneCeil);
1397 6 : const double dfBlueTimesToneBlue = std::min(dfBlue, dfToneCeil);
1398 : const double dfInvToneMaxRGB =
1399 6 : std::max(dfMaxRGB * dfInvToneCeil, 1.0);
1400 6 : const double dfReducedRGBTimesInvToneMaxRGB =
1401 : dfReducedRGB * dfInvToneMaxRGB;
1402 6 : padfDst[nRedBand] = std::min(
1403 6 : dfRedTimesToneRed * dfReducedRGBTimesInvToneMaxRGB, dfTopRGB);
1404 6 : padfDst[nGreenBand] =
1405 12 : std::min(dfGreenTimesToneGreen * dfReducedRGBTimesInvToneMaxRGB,
1406 6 : dfTopRGB);
1407 6 : padfDst[nBlueBand] = std::min(
1408 6 : dfBlueTimesToneBlue * dfReducedRGBTimesInvToneMaxRGB, dfTopRGB);
1409 : #else
1410 : // Original formulas. Slightly less optimized than the above ones.
1411 : const double dfToneMaxRGB = std::min(dfToneCeil / dfMaxRGB, 1.0);
1412 : const double dfToneRed = std::min(dfToneCeil / dfRed, 1.0);
1413 : const double dfToneGreen = std::min(dfToneCeil / dfGreen, 1.0);
1414 : const double dfToneBlue = std::min(dfToneCeil / dfBlue, 1.0);
1415 : padfDst[nRedBand] = std::min(
1416 : dfReducedRGB * dfRed * dfToneRed / dfToneMaxRGB, dfTopRGB);
1417 : padfDst[nGreenBand] = std::min(
1418 : dfReducedRGB * dfGreen * dfToneGreen / dfToneMaxRGB, dfTopRGB);
1419 : padfDst[nBlueBand] = std::min(
1420 : dfReducedRGB * dfBlue * dfToneBlue / dfToneMaxRGB, dfTopRGB);
1421 : #endif
1422 :
1423 : // Other bands processing (NIR, ...): only apply RGB reduction factor
1424 6 : if (bRGBBandsAreFirst)
1425 : {
1426 : // optimization
1427 12 : for (int iBand = 3; iBand < nInBands; ++iBand)
1428 : {
1429 6 : if (padfSrc[iBand] != padfInNoData[iBand])
1430 : {
1431 6 : padfDst[iBand] = dfReducedRGB * padfSrc[iBand];
1432 : }
1433 : else
1434 : {
1435 0 : bNoDataPixel = true;
1436 0 : break;
1437 : }
1438 : }
1439 : }
1440 : else
1441 : {
1442 0 : for (int iBand = 0; iBand < nInBands; ++iBand)
1443 : {
1444 0 : if (iBand != nRedBand && iBand != nGreenBand &&
1445 0 : iBand != nBlueBand)
1446 : {
1447 0 : if (padfSrc[iBand] != padfInNoData[iBand])
1448 : {
1449 0 : padfDst[iBand] = dfReducedRGB * padfSrc[iBand];
1450 : }
1451 : else
1452 : {
1453 0 : bNoDataPixel = true;
1454 0 : break;
1455 : }
1456 : }
1457 : }
1458 6 : }
1459 : }
1460 : else
1461 : {
1462 0 : bNoDataPixel = true;
1463 : }
1464 6 : if (bNoDataPixel)
1465 : {
1466 0 : for (int iBand = 0; iBand < nInBands; ++iBand)
1467 : {
1468 0 : padfDst[iBand] = padfOutNoData[iBand];
1469 : }
1470 : }
1471 :
1472 6 : padfSrc += nInBands;
1473 6 : padfDst += nInBands;
1474 : }
1475 :
1476 1 : return CE_None;
1477 : }
1478 :
1479 : /************************************************************************/
1480 : /* ExpressionInit() */
1481 : /************************************************************************/
1482 :
1483 : namespace
1484 : {
1485 :
1486 : class ExpressionData
1487 : {
1488 : public:
1489 7 : ExpressionData(int nInBands, int nBatchSize, std::string_view osExpression,
1490 : std::string_view osDialect)
1491 7 : : m_nInBands(nInBands), m_nNominalBatchSize(nBatchSize),
1492 7 : m_nBatchCount(DIV_ROUND_UP(nInBands, nBatchSize)), m_adfResults{},
1493 14 : m_osExpression(std::string(osExpression)),
1494 14 : m_osDialect(std::string(osDialect)), m_oNominalBatchEnv{},
1495 42 : m_oPartialBatchEnv{}
1496 : {
1497 7 : }
1498 :
1499 7 : CPLErr Compile()
1500 : {
1501 14 : auto eErr = m_oNominalBatchEnv.Initialize(m_osExpression, m_osDialect,
1502 7 : m_nNominalBatchSize);
1503 7 : if (eErr != CE_None)
1504 : {
1505 0 : return eErr;
1506 : }
1507 :
1508 7 : const auto nPartialBatchSize = m_nInBands % m_nNominalBatchSize;
1509 7 : if (nPartialBatchSize)
1510 : {
1511 1 : eErr = m_oPartialBatchEnv.Initialize(m_osExpression, m_osDialect,
1512 : nPartialBatchSize);
1513 : }
1514 :
1515 7 : return eErr;
1516 : }
1517 :
1518 14 : CPLErr Evaluate(const double *padfInputs, size_t nExpectedOutBands)
1519 : {
1520 14 : m_adfResults.clear();
1521 :
1522 60 : for (int iBatch = 0; iBatch < m_nBatchCount; iBatch++)
1523 : {
1524 46 : const auto nBandsRemaining =
1525 46 : static_cast<int>(m_nInBands - (m_nNominalBatchSize * iBatch));
1526 : const auto nBatchSize =
1527 46 : std::min(m_nNominalBatchSize, nBandsRemaining);
1528 :
1529 46 : auto &oEnv = GetEnv(nBatchSize);
1530 :
1531 46 : const double *pdfStart = padfInputs + iBatch * m_nNominalBatchSize;
1532 46 : const double *pdfEnd = pdfStart + nBatchSize;
1533 :
1534 46 : std::copy(pdfStart, pdfEnd, oEnv.m_adfValuesForPixel.begin());
1535 :
1536 46 : if (auto eErr = oEnv.m_poExpression->Evaluate(); eErr != CE_None)
1537 : {
1538 0 : return eErr;
1539 : }
1540 :
1541 46 : const auto &adfResults = oEnv.m_poExpression->Results();
1542 46 : if (m_nBatchCount > 1)
1543 : {
1544 : std::copy(adfResults.begin(), adfResults.end(),
1545 38 : std::back_inserter(m_adfResults));
1546 : }
1547 : }
1548 :
1549 14 : if (nExpectedOutBands > 0)
1550 : {
1551 10 : if (Results().size() != static_cast<std::size_t>(nExpectedOutBands))
1552 : {
1553 0 : CPLError(CE_Failure, CPLE_AppDefined,
1554 : "Expression returned %d values but "
1555 : "%d output bands were expected.",
1556 0 : static_cast<int>(Results().size()),
1557 : static_cast<int>(nExpectedOutBands));
1558 0 : return CE_Failure;
1559 : }
1560 : }
1561 :
1562 14 : return CE_None;
1563 : }
1564 :
1565 24 : const std::vector<double> &Results() const
1566 : {
1567 24 : if (m_nBatchCount == 1)
1568 : {
1569 15 : return m_oNominalBatchEnv.m_poExpression->Results();
1570 : }
1571 : else
1572 : {
1573 9 : return m_adfResults;
1574 : }
1575 : }
1576 :
1577 : private:
1578 : const int m_nInBands;
1579 : const int m_nNominalBatchSize;
1580 : const int m_nBatchCount;
1581 : std::vector<double> m_adfResults;
1582 :
1583 : const CPLString m_osExpression;
1584 : const CPLString m_osDialect;
1585 :
1586 : struct InvocationEnv
1587 : {
1588 : std::vector<double> m_adfValuesForPixel;
1589 : std::unique_ptr<gdal::MathExpression> m_poExpression;
1590 :
1591 8 : CPLErr Initialize(const CPLString &osExpression,
1592 : const CPLString &osDialect, int nBatchSize)
1593 : {
1594 : m_poExpression =
1595 8 : gdal::MathExpression::Create(osExpression, osDialect.c_str());
1596 : // cppcheck-suppress knownConditionTrueFalse
1597 8 : if (m_poExpression == nullptr)
1598 : {
1599 0 : return CE_Failure;
1600 : }
1601 :
1602 8 : m_adfValuesForPixel.resize(nBatchSize);
1603 :
1604 40 : for (int i = 0; i < nBatchSize; i++)
1605 : {
1606 64 : std::string osVar = "B" + std::to_string(i + 1);
1607 64 : m_poExpression->RegisterVariable(osVar,
1608 32 : &m_adfValuesForPixel[i]);
1609 : }
1610 :
1611 8 : if (osExpression.ifind("BANDS") != std::string::npos)
1612 : {
1613 6 : m_poExpression->RegisterVector("BANDS", &m_adfValuesForPixel);
1614 : }
1615 :
1616 8 : return m_poExpression->Compile();
1617 : }
1618 : };
1619 :
1620 46 : InvocationEnv &GetEnv(int nBatchSize)
1621 : {
1622 46 : if (nBatchSize == m_nNominalBatchSize)
1623 : {
1624 44 : return m_oNominalBatchEnv;
1625 : }
1626 : else
1627 : {
1628 2 : return m_oPartialBatchEnv;
1629 : }
1630 : }
1631 :
1632 : InvocationEnv m_oNominalBatchEnv;
1633 : InvocationEnv m_oPartialBatchEnv;
1634 : };
1635 :
1636 : } // namespace
1637 :
1638 7 : static CPLErr ExpressionInit(const char * /*pszFuncName*/, void * /*pUserData*/,
1639 : CSLConstList papszFunctionArgs, int nInBands,
1640 : GDALDataType eInDT, double * /* padfInNoData */,
1641 : int *pnOutBands, GDALDataType *peOutDT,
1642 : double ** /* ppadfOutNoData */,
1643 : const char * /* pszVRTPath */,
1644 : VRTPDWorkingDataPtr *ppWorkingData)
1645 : {
1646 7 : CPLAssert(eInDT == GDT_Float64);
1647 :
1648 7 : *peOutDT = eInDT;
1649 7 : *ppWorkingData = nullptr;
1650 :
1651 : const char *pszBatchSize =
1652 7 : CSLFetchNameValue(papszFunctionArgs, "batch_size");
1653 7 : auto nBatchSize = nInBands;
1654 :
1655 7 : if (pszBatchSize != nullptr)
1656 : {
1657 4 : nBatchSize = std::min(nInBands, std::atoi(pszBatchSize));
1658 : }
1659 :
1660 7 : if (nBatchSize < 1)
1661 : {
1662 0 : CPLError(CE_Failure, CPLE_IllegalArg, "batch_size must be at least 1");
1663 0 : return CE_Failure;
1664 : }
1665 :
1666 7 : const char *pszDialect = CSLFetchNameValue(papszFunctionArgs, "dialect");
1667 7 : if (pszDialect == nullptr)
1668 : {
1669 0 : pszDialect = "muparser";
1670 : }
1671 :
1672 : const char *pszExpression =
1673 7 : CSLFetchNameValue(papszFunctionArgs, "expression");
1674 :
1675 : auto data = std::make_unique<ExpressionData>(nInBands, nBatchSize,
1676 14 : pszExpression, pszDialect);
1677 :
1678 7 : if (auto eErr = data->Compile(); eErr != CE_None)
1679 : {
1680 0 : return eErr;
1681 : }
1682 :
1683 7 : if (*pnOutBands == 0)
1684 : {
1685 4 : std::vector<double> aDummyValues(nInBands);
1686 4 : if (auto eErr = data->Evaluate(aDummyValues.data(), 0); eErr != CE_None)
1687 : {
1688 0 : return eErr;
1689 : }
1690 :
1691 4 : *pnOutBands = static_cast<int>(data->Results().size());
1692 : }
1693 :
1694 7 : *ppWorkingData = data.release();
1695 :
1696 7 : return CE_None;
1697 : }
1698 :
1699 7 : static void ExpressionFree(const char * /* pszFuncName */,
1700 : void * /* pUserData */,
1701 : VRTPDWorkingDataPtr pWorkingData)
1702 : {
1703 7 : ExpressionData *data = static_cast<ExpressionData *>(pWorkingData);
1704 7 : delete data;
1705 7 : }
1706 :
1707 7 : static CPLErr ExpressionProcess(
1708 : const char * /* pszFuncName */, void * /* pUserData */,
1709 : VRTPDWorkingDataPtr pWorkingData, CSLConstList /* papszFunctionArgs */,
1710 : int nBufXSize, int nBufYSize, const void *pInBuffer,
1711 : size_t /* nInBufferSize */, GDALDataType eInDT, int nInBands,
1712 : const double *CPL_RESTRICT /* padfInNoData */, void *pOutBuffer,
1713 : size_t /* nOutBufferSize */, GDALDataType eOutDT, int nOutBands,
1714 : const double *CPL_RESTRICT /* padfOutNoData */, double /* dfSrcXOff */,
1715 : double /* dfSrcYOff */, double /* dfSrcXSize */, double /* dfSrcYSize */,
1716 : const double /* adfSrcGT */[], const char * /* pszVRTPath "*/,
1717 : CSLConstList /* papszExtra */)
1718 : {
1719 7 : ExpressionData *expr = static_cast<ExpressionData *>(pWorkingData);
1720 :
1721 7 : const size_t nElts = static_cast<size_t>(nBufXSize) * nBufYSize;
1722 :
1723 7 : CPL_IGNORE_RET_VAL(eInDT);
1724 7 : CPLAssert(eInDT == GDT_Float64);
1725 7 : const double *CPL_RESTRICT padfSrc = static_cast<const double *>(pInBuffer);
1726 :
1727 7 : CPLAssert(eOutDT == GDT_Float64);
1728 7 : CPL_IGNORE_RET_VAL(eOutDT);
1729 7 : double *CPL_RESTRICT padfDst = static_cast<double *>(pOutBuffer);
1730 :
1731 17 : for (size_t i = 0; i < nElts; i++)
1732 : {
1733 10 : if (auto eErr = expr->Evaluate(padfSrc, nOutBands); eErr != CE_None)
1734 : {
1735 0 : return eErr;
1736 : }
1737 :
1738 10 : const auto &adfResults = expr->Results();
1739 10 : std::copy(adfResults.begin(), adfResults.end(), padfDst);
1740 :
1741 10 : padfDst += nOutBands;
1742 10 : padfSrc += nInBands;
1743 : }
1744 :
1745 7 : return CE_None;
1746 : }
1747 :
1748 : /************************************************************************/
1749 : /* GDALVRTRegisterDefaultProcessedDatasetFuncs() */
1750 : /************************************************************************/
1751 :
1752 : /** Register builtin functions that can be used in a VRTProcessedDataset.
1753 : */
1754 1384 : void GDALVRTRegisterDefaultProcessedDatasetFuncs()
1755 : {
1756 1384 : GDALVRTRegisterProcessedDatasetFunc(
1757 : "BandAffineCombination", nullptr,
1758 : "<ProcessedDatasetFunctionArgumentsList>"
1759 : " <Argument name='src_nodata' type='double' "
1760 : "description='Override input nodata value'/>"
1761 : " <Argument name='dst_nodata' type='double' "
1762 : "description='Override output nodata value'/>"
1763 : " <Argument name='replacement_nodata' "
1764 : "description='value to substitute to a valid computed value that "
1765 : "would be nodata' type='double'/>"
1766 : " <Argument name='dst_intended_datatype' type='string' "
1767 : "description='Intented datatype of output (which might be "
1768 : "different than the working data type)'/>"
1769 : " <Argument name='coefficients_{band}' "
1770 : "description='Comma-separated coefficients for combining bands. "
1771 : "First one is constant term' "
1772 : "type='double_list' required='true'/>"
1773 : " <Argument name='min' description='clamp min value' type='double'/>"
1774 : " <Argument name='max' description='clamp max value' type='double'/>"
1775 : "</ProcessedDatasetFunctionArgumentsList>",
1776 : GDT_Float64, nullptr, 0, nullptr, 0, BandAffineCombinationInit,
1777 : BandAffineCombinationFree, BandAffineCombinationProcess, nullptr);
1778 :
1779 1384 : GDALVRTRegisterProcessedDatasetFunc(
1780 : "LUT", nullptr,
1781 : "<ProcessedDatasetFunctionArgumentsList>"
1782 : " <Argument name='src_nodata' type='double' "
1783 : "description='Override input nodata value'/>"
1784 : " <Argument name='dst_nodata' type='double' "
1785 : "description='Override output nodata value'/>"
1786 : " <Argument name='lut_{band}' "
1787 : "description='List of the form [src value 1]:[dest value 1],"
1788 : "[src value 2]:[dest value 2],...' "
1789 : "type='string' required='true'/>"
1790 : "</ProcessedDatasetFunctionArgumentsList>",
1791 : GDT_Float64, nullptr, 0, nullptr, 0, LUTInit, LUTFree, LUTProcess,
1792 : nullptr);
1793 :
1794 1384 : GDALVRTRegisterProcessedDatasetFunc(
1795 : "LocalScaleOffset", nullptr,
1796 : "<ProcessedDatasetFunctionArgumentsList>"
1797 : " <Argument name='relativeToVRT' "
1798 : "description='Whether gain and offset filenames are relative to "
1799 : "the VRT' type='boolean' default='false'/>"
1800 : " <Argument name='gain_dataset_filename_{band}' "
1801 : "description='Filename to the gain dataset' "
1802 : "type='string' required='true'/>"
1803 : " <Argument name='gain_dataset_band_{band}' "
1804 : "description='Band of the gain dataset' "
1805 : "type='integer' required='true'/>"
1806 : " <Argument name='offset_dataset_filename_{band}' "
1807 : "description='Filename to the offset dataset' "
1808 : "type='string' required='true'/>"
1809 : " <Argument name='offset_dataset_band_{band}' "
1810 : "description='Band of the offset dataset' "
1811 : "type='integer' required='true'/>"
1812 : " <Argument name='min' description='clamp min value' type='double'/>"
1813 : " <Argument name='max' description='clamp max value' type='double'/>"
1814 : " <Argument name='nodata' type='double' "
1815 : "description='Override dataset nodata value'/>"
1816 : " <Argument name='gain_nodata' type='double' "
1817 : "description='Override gain dataset nodata value'/>"
1818 : " <Argument name='offset_nodata' type='double' "
1819 : "description='Override offset dataset nodata value'/>"
1820 : "</ProcessedDatasetFunctionArgumentsList>",
1821 : GDT_Float64, nullptr, 0, nullptr, 0, LocalScaleOffsetInit,
1822 : LocalScaleOffsetFree, LocalScaleOffsetProcess, nullptr);
1823 :
1824 1384 : GDALVRTRegisterProcessedDatasetFunc(
1825 : "Trimming", nullptr,
1826 : "<ProcessedDatasetFunctionArgumentsList>"
1827 : " <Argument name='relativeToVRT' "
1828 : "description='Whether trimming_dataset_filename is relative to the VRT'"
1829 : " type='boolean' default='false'/>"
1830 : " <Argument name='trimming_dataset_filename' "
1831 : "description='Filename to the trimming dataset' "
1832 : "type='string' required='true'/>"
1833 : " <Argument name='red_band' type='integer' default='1'/>"
1834 : " <Argument name='green_band' type='integer' default='2'/>"
1835 : " <Argument name='blue_band' type='integer' default='3'/>"
1836 : " <Argument name='top_rgb' "
1837 : "description='Maximum saturating RGB output value' "
1838 : "type='double' required='true'/>"
1839 : " <Argument name='tone_ceil' "
1840 : "description='Maximum threshold beyond which we give up saturation' "
1841 : "type='double' required='true'/>"
1842 : " <Argument name='top_margin' "
1843 : "description='Margin to allow for dynamics in brighest areas "
1844 : "(between 0 and 1, should be close to 0)' "
1845 : "type='double' required='true'/>"
1846 : " <Argument name='nodata' type='double' "
1847 : "description='Override dataset nodata value'/>"
1848 : " <Argument name='trimming_nodata' type='double' "
1849 : "description='Override trimming dataset nodata value'/>"
1850 : "</ProcessedDatasetFunctionArgumentsList>",
1851 : GDT_Float64, nullptr, 0, nullptr, 0, TrimmingInit, TrimmingFree,
1852 : TrimmingProcess, nullptr);
1853 :
1854 1384 : GDALVRTRegisterProcessedDatasetFunc(
1855 : "Expression", nullptr,
1856 : "<ProcessedDatasetFunctionArgumentsList>"
1857 : " <Argument name='expression' description='the expression to "
1858 : "evaluate' type='string' required='true' />"
1859 : " <Argument name='dialect' description='expression dialect' "
1860 : "type='string' />"
1861 : " <Argument name='batch_size' description='batch size' "
1862 : "type='integer' />"
1863 : "</ProcessedDatasetFunctionArgumentsList>",
1864 : GDT_Float64, nullptr, 0, nullptr, 0, ExpressionInit, ExpressionFree,
1865 : ExpressionProcess, nullptr);
1866 1384 : }
|
