Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Virtual GDAL Datasets
4 : * Purpose: Implementation of Reclassifier
5 : * Author: Daniel Baston
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2025, ISciences LLC
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "cpl_conv.h"
14 : #include "vrtreclassifier.h"
15 :
16 : #include <algorithm>
17 : #include <cmath>
18 : #include <limits>
19 :
20 : namespace gdal
21 : {
22 :
23 2071 : bool Reclassifier::Interval::Overlaps(const Interval &other) const
24 : {
25 2071 : if (dfMin > other.dfMax || dfMax < other.dfMin)
26 : {
27 2069 : return false;
28 : }
29 :
30 2 : return true;
31 : }
32 :
33 2101 : CPLErr Reclassifier::Interval::Parse(const char *s, char **rest)
34 : {
35 2101 : const char *start = s;
36 : bool bMinIncluded;
37 : bool bMaxIncluded;
38 :
39 2101 : while (isspace(*start))
40 : {
41 0 : start++;
42 : }
43 :
44 : char *end;
45 :
46 2101 : if (*start == '(')
47 : {
48 46 : bMinIncluded = false;
49 : }
50 2055 : else if (*start == '[')
51 : {
52 2023 : bMinIncluded = true;
53 : }
54 : else
55 : {
56 32 : double dfVal = CPLStrtod(start, &end);
57 :
58 32 : if (end == start)
59 : {
60 2 : CPLError(CE_Failure, CPLE_AppDefined,
61 : "Interval must start with '(' or ']'");
62 2 : return CE_Failure;
63 : }
64 :
65 30 : SetToConstant(dfVal);
66 :
67 30 : if (rest != nullptr)
68 : {
69 30 : *rest = end;
70 : }
71 :
72 30 : return CE_None;
73 : }
74 2069 : start++;
75 :
76 2069 : while (isspace(*start))
77 : {
78 0 : start++;
79 : }
80 :
81 2069 : if (STARTS_WITH_CI(start, "-inf"))
82 : {
83 22 : dfMin = -std::numeric_limits<double>::infinity();
84 22 : end = const_cast<char *>(start + 4);
85 : }
86 : else
87 : {
88 2047 : dfMin = CPLStrtod(start, &end);
89 : }
90 :
91 2069 : if (end == start || *end != ',')
92 : {
93 0 : CPLError(CE_Failure, CPLE_AppDefined, "Expected a number");
94 0 : return CE_Failure;
95 : }
96 2069 : start = end + 1;
97 :
98 4129 : while (isspace(*start))
99 : {
100 2060 : start++;
101 : }
102 :
103 2069 : if (STARTS_WITH_CI(start, "inf"))
104 : {
105 14 : dfMax = std::numeric_limits<double>::infinity();
106 14 : end = const_cast<char *>(start + 3);
107 : }
108 : else
109 : {
110 2055 : dfMax = CPLStrtod(start, &end);
111 : }
112 :
113 2069 : if (end == start || (*end != ')' && *end != ']'))
114 : {
115 2 : CPLError(CE_Failure, CPLE_AppDefined,
116 : "Interval must end with ')' or ']");
117 2 : return CE_Failure;
118 : }
119 2067 : if (*end == ')')
120 : {
121 2026 : bMaxIncluded = false;
122 : }
123 : else
124 : {
125 41 : bMaxIncluded = true;
126 : }
127 :
128 2067 : if (rest != nullptr)
129 : {
130 2067 : *rest = end + 1;
131 : }
132 :
133 2067 : if (std::isnan(dfMin) || std::isnan(dfMax))
134 : {
135 2 : CPLError(CE_Failure, CPLE_AppDefined,
136 : "NaN is not a valid value for bounds of interval");
137 2 : return CE_Failure;
138 : }
139 :
140 2065 : if (dfMin > dfMax)
141 : {
142 1 : CPLError(
143 : CE_Failure, CPLE_AppDefined,
144 : "Lower bound of interval must be lower or equal to upper bound");
145 1 : return CE_Failure;
146 : }
147 :
148 2064 : if (!bMinIncluded)
149 : {
150 44 : dfMin = std::nextafter(dfMin, std::numeric_limits<double>::infinity());
151 : }
152 2064 : if (!bMaxIncluded)
153 : {
154 2026 : dfMax = std::nextafter(dfMax, -std::numeric_limits<double>::infinity());
155 : }
156 :
157 2064 : return CE_None;
158 : }
159 :
160 44 : void Reclassifier::Interval::SetToConstant(double dfVal)
161 : {
162 44 : dfMin = dfVal;
163 44 : dfMax = dfVal;
164 44 : }
165 :
166 30 : CPLErr Reclassifier::Finalize()
167 : {
168 30 : std::sort(m_aoIntervalMappings.begin(), m_aoIntervalMappings.end(),
169 21410 : [](const auto &a, const auto &b)
170 21410 : { return a.first.dfMin < b.first.dfMin; });
171 :
172 2099 : for (std::size_t i = 1; i < m_aoIntervalMappings.size(); i++)
173 : {
174 4142 : if (m_aoIntervalMappings[i - 1].first.Overlaps(
175 2071 : m_aoIntervalMappings[i].first))
176 : {
177 : // Don't use [, ) notation because we will have modified those values for an open interval
178 2 : CPLError(CE_Failure, CPLE_AppDefined,
179 : "Interval from %g to %g (mapped to %g) overlaps with "
180 : "interval from %g to %g (mapped to %g)",
181 2 : m_aoIntervalMappings[i - 1].first.dfMin,
182 2 : m_aoIntervalMappings[i - 1].first.dfMax,
183 2 : m_aoIntervalMappings[i - 1].second.value_or(
184 2 : std::numeric_limits<double>::quiet_NaN()),
185 2 : m_aoIntervalMappings[i].first.dfMin,
186 2 : m_aoIntervalMappings[i].first.dfMax,
187 2 : m_aoIntervalMappings[i].second.value_or(
188 2 : std::numeric_limits<double>::quiet_NaN()));
189 2 : return CE_Failure;
190 : }
191 : }
192 :
193 28 : return CE_None;
194 : }
195 :
196 2103 : void Reclassifier::AddMapping(const Interval &interval,
197 : std::optional<double> dfDstVal)
198 : {
199 2103 : m_aoIntervalMappings.emplace_back(interval, dfDstVal);
200 2103 : }
201 :
202 43 : CPLErr Reclassifier::Init(const char *pszText,
203 : std::optional<double> noDataValue,
204 : GDALDataType eBufType)
205 : {
206 43 : const char *start = pszText;
207 43 : char *end = const_cast<char *>(start);
208 :
209 2161 : while (*end != '\0')
210 : {
211 2197 : while (isspace(*start))
212 : {
213 66 : start++;
214 : }
215 :
216 2131 : Interval sInt{};
217 2131 : bool bFromIsDefault = false;
218 2131 : bool bPassThrough = false;
219 2131 : bool bFromNaN = false;
220 :
221 2131 : if (STARTS_WITH_CI(start, "DEFAULT"))
222 : {
223 14 : bFromIsDefault = true;
224 14 : end = const_cast<char *>(start + 7);
225 : }
226 2117 : else if (STARTS_WITH_CI(start, "NO_DATA"))
227 : {
228 15 : if (!noDataValue.has_value())
229 : {
230 1 : CPLError(
231 : CE_Failure, CPLE_AppDefined,
232 : "Value mapped from NO_DATA, but NoData value is not set");
233 13 : return CE_Failure;
234 : }
235 :
236 14 : sInt.SetToConstant(noDataValue.value());
237 14 : end = const_cast<char *>(start + 7);
238 : }
239 2102 : else if (STARTS_WITH_CI(start, "NAN"))
240 : {
241 1 : bFromNaN = true;
242 1 : end = const_cast<char *>(start + 3);
243 : }
244 : else
245 : {
246 2101 : if (auto eErr = sInt.Parse(start, &end); eErr != CE_None)
247 : {
248 7 : return eErr;
249 : }
250 : }
251 :
252 4295 : while (isspace(*end))
253 : {
254 2172 : end++;
255 : }
256 :
257 2123 : if (*end != MAPPING_FROMTO_SEP_CHAR)
258 : {
259 1 : CPLError(CE_Failure, CPLE_AppDefined,
260 : "Failed to parse mapping (expected '%c', got '%c')",
261 1 : MAPPING_FROMTO_SEP_CHAR, *end);
262 1 : return CE_Failure;
263 : }
264 :
265 2122 : start = end + 1;
266 :
267 4185 : while (isspace(*start))
268 : {
269 2063 : start++;
270 : }
271 :
272 2122 : std::optional<double> dfDstVal{};
273 2122 : if (STARTS_WITH(start, "NO_DATA"))
274 : {
275 16 : if (!noDataValue.has_value())
276 : {
277 1 : CPLError(
278 : CE_Failure, CPLE_AppDefined,
279 : "Value mapped to NO_DATA, but NoData value is not set");
280 1 : return CE_Failure;
281 : }
282 15 : dfDstVal = noDataValue.value();
283 15 : end = const_cast<char *>(start) + 7;
284 : }
285 2106 : else if (STARTS_WITH(start, "PASS_THROUGH"))
286 : {
287 15 : bPassThrough = true;
288 15 : end = const_cast<char *>(start + 12);
289 : }
290 : else
291 : {
292 2091 : dfDstVal = CPLStrtod(start, &end);
293 2091 : if (start == end)
294 : {
295 1 : CPLError(CE_Failure, CPLE_AppDefined,
296 : "Failed to parse output value (expected number or "
297 : "NO_DATA)");
298 1 : return CE_Failure;
299 : }
300 : }
301 :
302 2135 : while (isspace(*end))
303 : {
304 15 : end++;
305 : }
306 :
307 2120 : if (*end != '\0' && *end != MAPPING_INTERVAL_SEP_CHAR)
308 : {
309 1 : CPLError(CE_Failure, CPLE_AppDefined,
310 : "Failed to parse mapping (expected '%c' or end of string, "
311 : "got '%c')",
312 1 : MAPPING_INTERVAL_SEP_CHAR, *end);
313 1 : return CE_Failure;
314 : }
315 :
316 4223 : if (dfDstVal.has_value() &&
317 2104 : !GDALIsValueExactAs(dfDstVal.value(), eBufType))
318 : {
319 1 : CPLError(CE_Failure, CPLE_AppDefined,
320 : "Value %g cannot be represented as data type %s",
321 1 : dfDstVal.value(), GDALGetDataTypeName(eBufType));
322 1 : return CE_Failure;
323 : }
324 :
325 2118 : if (bFromNaN)
326 : {
327 2 : SetNaNValue(bPassThrough ? std::numeric_limits<double>::quiet_NaN()
328 1 : : dfDstVal.value());
329 : }
330 2117 : else if (bFromIsDefault)
331 : {
332 14 : if (bPassThrough)
333 : {
334 1 : SetDefaultPassThrough(true);
335 : }
336 : else
337 : {
338 13 : SetDefaultValue(dfDstVal.value());
339 : }
340 : }
341 : else
342 : {
343 2103 : AddMapping(sInt, dfDstVal);
344 : }
345 :
346 2118 : start = end + 1;
347 : }
348 :
349 30 : return Finalize();
350 : }
351 :
352 20388 : static std::optional<size_t> FindInterval(
353 : const std::vector<std::pair<Reclassifier::Interval, std::optional<double>>>
354 : &arr,
355 : double srcVal)
356 : {
357 20388 : if (arr.empty())
358 : {
359 1 : return std::nullopt;
360 : }
361 :
362 20387 : size_t low = 0;
363 20387 : size_t high = arr.size() - 1;
364 :
365 105629 : while (low <= high)
366 : {
367 104705 : auto mid = low + (high - low) / 2;
368 :
369 104705 : const auto &mid_interval = arr[mid].first;
370 104705 : if (mid_interval.Contains(srcVal))
371 : {
372 19462 : return mid;
373 : }
374 :
375 : // Could an interval exist to the left?
376 85243 : if (srcVal < mid_interval.dfMin)
377 : {
378 40366 : if (mid == 0)
379 : {
380 1 : return std::nullopt;
381 : }
382 40365 : high = mid - 1;
383 : }
384 : // Could an interval exist to the right?
385 44877 : else if (srcVal > mid_interval.dfMax)
386 : {
387 44877 : low = mid + 1;
388 : }
389 : else
390 : {
391 0 : return std::nullopt;
392 : }
393 : }
394 :
395 924 : return std::nullopt;
396 : }
397 :
398 20389 : double Reclassifier::Reclassify(double srcVal, bool &bFoundInterval) const
399 : {
400 20389 : bFoundInterval = false;
401 :
402 20389 : if (std::isnan(srcVal))
403 : {
404 1 : if (m_NaNValue.has_value())
405 : {
406 1 : bFoundInterval = true;
407 1 : return m_NaNValue.value();
408 : }
409 : }
410 : else
411 : {
412 20388 : auto nInterval = FindInterval(m_aoIntervalMappings, srcVal);
413 20388 : if (nInterval.has_value())
414 : {
415 19462 : bFoundInterval = true;
416 19462 : return m_aoIntervalMappings[nInterval.value()].second.value_or(
417 19462 : srcVal);
418 : }
419 : }
420 :
421 926 : if (m_defaultValue.has_value())
422 : {
423 921 : bFoundInterval = true;
424 921 : return m_defaultValue.value();
425 : }
426 :
427 5 : if (m_defaultPassThrough)
428 : {
429 3 : bFoundInterval = true;
430 3 : return srcVal;
431 : }
432 :
433 2 : return 0;
434 : }
435 :
436 : } // namespace gdal
|
