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