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