Line data Source code
1 : /******************************************************************************
2 : *
3 : * Name: gdalmultidim_array_mask.cpp
4 : * Project: GDAL Core
5 : * Purpose: GDALMDArray::GetMask() implementation
6 : * Author: Even Rouault <even.rouault at spatialys.com>
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2019, Even Rouault <even.rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_float.h"
15 : #include "gdal_multidim.h"
16 : #include "gdal_pam_multidim.h"
17 :
18 : #include <algorithm>
19 : #include <cmath>
20 :
21 : //! @cond Doxygen_Suppress
22 :
23 : /************************************************************************/
24 : /* GDALMDArrayMask */
25 : /************************************************************************/
26 :
27 : class GDALMDArrayMask final : public GDALPamMDArray
28 : {
29 : private:
30 : std::shared_ptr<GDALMDArray> m_poParent{};
31 : GDALExtendedDataType m_dt{GDALExtendedDataType::Create(GDT_UInt8)};
32 : double m_dfMissingValue = 0.0;
33 : bool m_bHasMissingValue = false;
34 : double m_dfFillValue = 0.0;
35 : bool m_bHasFillValue = false;
36 : double m_dfValidMin = 0.0;
37 : bool m_bHasValidMin = false;
38 : double m_dfValidMax = 0.0;
39 : bool m_bHasValidMax = false;
40 : std::vector<uint32_t> m_anValidFlagMasks{};
41 : std::vector<uint32_t> m_anValidFlagValues{};
42 :
43 : bool Init(CSLConstList papszOptions);
44 :
45 : template <typename Type>
46 : void
47 : ReadInternal(const size_t *count, const GPtrDiff_t *bufferStride,
48 : const GDALExtendedDataType &bufferDataType, void *pDstBuffer,
49 : const void *pTempBuffer,
50 : const GDALExtendedDataType &oTmpBufferDT,
51 : const std::vector<GPtrDiff_t> &tmpBufferStrideVector) const;
52 :
53 : protected:
54 48 : explicit GDALMDArrayMask(const std::shared_ptr<GDALMDArray> &poParent)
55 96 : : GDALAbstractMDArray(std::string(),
56 96 : "Mask of " + poParent->GetFullName()),
57 96 : GDALPamMDArray(std::string(), "Mask of " + poParent->GetFullName(),
58 96 : GDALPamMultiDim::GetPAM(poParent),
59 : poParent->GetContext()),
60 240 : m_poParent(std::move(poParent))
61 : {
62 48 : }
63 :
64 : bool IRead(const GUInt64 *arrayStartIdx, const size_t *count,
65 : const GInt64 *arrayStep, const GPtrDiff_t *bufferStride,
66 : const GDALExtendedDataType &bufferDataType,
67 : void *pDstBuffer) const override;
68 :
69 0 : bool IAdviseRead(const GUInt64 *arrayStartIdx, const size_t *count,
70 : CSLConstList papszOptions) const override
71 : {
72 0 : return m_poParent->AdviseRead(arrayStartIdx, count, papszOptions);
73 : }
74 :
75 : public:
76 : static std::shared_ptr<GDALMDArrayMask>
77 : Create(const std::shared_ptr<GDALMDArray> &poParent,
78 : CSLConstList papszOptions);
79 :
80 1 : bool IsWritable() const override
81 : {
82 1 : return false;
83 : }
84 :
85 53 : const std::string &GetFilename() const override
86 : {
87 53 : return m_poParent->GetFilename();
88 : }
89 :
90 : const std::vector<std::shared_ptr<GDALDimension>> &
91 383 : GetDimensions() const override
92 : {
93 383 : return m_poParent->GetDimensions();
94 : }
95 :
96 138 : const GDALExtendedDataType &GetDataType() const override
97 : {
98 138 : return m_dt;
99 : }
100 :
101 1 : std::shared_ptr<OGRSpatialReference> GetSpatialRef() const override
102 : {
103 1 : return m_poParent->GetSpatialRef();
104 : }
105 :
106 2 : std::vector<GUInt64> GetBlockSize() const override
107 : {
108 2 : return m_poParent->GetBlockSize();
109 : }
110 : };
111 :
112 : /************************************************************************/
113 : /* GDALMDArrayMask::Create() */
114 : /************************************************************************/
115 :
116 : /* static */ std::shared_ptr<GDALMDArrayMask>
117 48 : GDALMDArrayMask::Create(const std::shared_ptr<GDALMDArray> &poParent,
118 : CSLConstList papszOptions)
119 : {
120 96 : auto newAr(std::shared_ptr<GDALMDArrayMask>(new GDALMDArrayMask(poParent)));
121 48 : newAr->SetSelf(newAr);
122 48 : if (!newAr->Init(papszOptions))
123 6 : return nullptr;
124 42 : return newAr;
125 : }
126 :
127 : /************************************************************************/
128 : /* GDALMDArrayMask::Init() */
129 : /************************************************************************/
130 :
131 48 : bool GDALMDArrayMask::Init(CSLConstList papszOptions)
132 : {
133 : const auto GetSingleValNumericAttr =
134 192 : [this](const char *pszAttrName, bool &bHasVal, double &dfVal)
135 : {
136 576 : auto poAttr = m_poParent->GetAttribute(pszAttrName);
137 192 : if (poAttr && poAttr->GetDataType().GetClass() == GEDTC_NUMERIC)
138 : {
139 22 : const auto anDimSizes = poAttr->GetDimensionsSize();
140 21 : if (anDimSizes.empty() ||
141 10 : (anDimSizes.size() == 1 && anDimSizes[0] == 1))
142 : {
143 11 : bHasVal = true;
144 11 : dfVal = poAttr->ReadAsDouble();
145 : }
146 : }
147 192 : };
148 :
149 48 : GetSingleValNumericAttr("missing_value", m_bHasMissingValue,
150 48 : m_dfMissingValue);
151 48 : GetSingleValNumericAttr("_FillValue", m_bHasFillValue, m_dfFillValue);
152 48 : GetSingleValNumericAttr("valid_min", m_bHasValidMin, m_dfValidMin);
153 48 : GetSingleValNumericAttr("valid_max", m_bHasValidMax, m_dfValidMax);
154 :
155 : {
156 144 : auto poValidRange = m_poParent->GetAttribute("valid_range");
157 54 : if (poValidRange && poValidRange->GetDimensionsSize().size() == 1 &&
158 60 : poValidRange->GetDimensionsSize()[0] == 2 &&
159 6 : poValidRange->GetDataType().GetClass() == GEDTC_NUMERIC)
160 : {
161 6 : m_bHasValidMin = true;
162 6 : m_bHasValidMax = true;
163 6 : auto vals = poValidRange->ReadAsDoubleArray();
164 6 : CPLAssert(vals.size() == 2);
165 6 : m_dfValidMin = vals[0];
166 6 : m_dfValidMax = vals[1];
167 : }
168 : }
169 :
170 : // Take into account
171 : // https://cfconventions.org/cf-conventions/cf-conventions.html#flags
172 : // Cf GDALMDArray::GetMask() for semantics of UNMASK_FLAGS
173 : const char *pszUnmaskFlags =
174 48 : CSLFetchNameValue(papszOptions, "UNMASK_FLAGS");
175 48 : if (pszUnmaskFlags)
176 : {
177 : const auto IsScalarStringAttr =
178 13 : [](const std::shared_ptr<GDALAttribute> &poAttr)
179 : {
180 26 : return poAttr->GetDataType().GetClass() == GEDTC_STRING &&
181 26 : (poAttr->GetDimensionsSize().empty() ||
182 13 : (poAttr->GetDimensionsSize().size() == 1 &&
183 26 : poAttr->GetDimensionsSize()[0] == 1));
184 : };
185 :
186 28 : auto poFlagMeanings = m_poParent->GetAttribute("flag_meanings");
187 14 : if (!(poFlagMeanings && IsScalarStringAttr(poFlagMeanings)))
188 : {
189 1 : CPLError(CE_Failure, CPLE_AppDefined,
190 : "UNMASK_FLAGS option specified but array has no "
191 : "flag_meanings attribute");
192 1 : return false;
193 : }
194 13 : const char *pszFlagMeanings = poFlagMeanings->ReadAsString();
195 13 : if (!pszFlagMeanings)
196 : {
197 1 : CPLError(CE_Failure, CPLE_AppDefined,
198 : "Cannot read flag_meanings attribute");
199 1 : return false;
200 : }
201 :
202 : const auto IsSingleDimNumericAttr =
203 13 : [](const std::shared_ptr<GDALAttribute> &poAttr)
204 : {
205 26 : return poAttr->GetDataType().GetClass() == GEDTC_NUMERIC &&
206 26 : poAttr->GetDimensionsSize().size() == 1;
207 : };
208 :
209 24 : auto poFlagValues = m_poParent->GetAttribute("flag_values");
210 : const bool bHasFlagValues =
211 12 : poFlagValues && IsSingleDimNumericAttr(poFlagValues);
212 :
213 24 : auto poFlagMasks = m_poParent->GetAttribute("flag_masks");
214 : const bool bHasFlagMasks =
215 12 : poFlagMasks && IsSingleDimNumericAttr(poFlagMasks);
216 :
217 12 : if (!bHasFlagValues && !bHasFlagMasks)
218 : {
219 1 : CPLError(CE_Failure, CPLE_AppDefined,
220 : "Cannot find flag_values and/or flag_masks attribute");
221 1 : return false;
222 : }
223 :
224 : const CPLStringList aosUnmaskFlags(
225 11 : CSLTokenizeString2(pszUnmaskFlags, ",", 0));
226 : const CPLStringList aosFlagMeanings(
227 11 : CSLTokenizeString2(pszFlagMeanings, " ", 0));
228 :
229 11 : if (bHasFlagValues)
230 : {
231 7 : const auto eType = poFlagValues->GetDataType().GetNumericDataType();
232 : // We could support Int64 or UInt64, but more work...
233 7 : if (eType != GDT_UInt8 && eType != GDT_Int8 &&
234 7 : eType != GDT_UInt16 && eType != GDT_Int16 &&
235 0 : eType != GDT_UInt32 && eType != GDT_Int32)
236 : {
237 0 : CPLError(CE_Failure, CPLE_NotSupported,
238 : "Unsupported data type for flag_values attribute: %s",
239 : GDALGetDataTypeName(eType));
240 0 : return false;
241 : }
242 : }
243 :
244 11 : if (bHasFlagMasks)
245 : {
246 6 : const auto eType = poFlagMasks->GetDataType().GetNumericDataType();
247 : // We could support Int64 or UInt64, but more work...
248 6 : if (eType != GDT_UInt8 && eType != GDT_Int8 &&
249 6 : eType != GDT_UInt16 && eType != GDT_Int16 &&
250 0 : eType != GDT_UInt32 && eType != GDT_Int32)
251 : {
252 0 : CPLError(CE_Failure, CPLE_NotSupported,
253 : "Unsupported data type for flag_masks attribute: %s",
254 : GDALGetDataTypeName(eType));
255 0 : return false;
256 : }
257 : }
258 :
259 : const std::vector<double> adfValues(
260 : bHasFlagValues ? poFlagValues->ReadAsDoubleArray()
261 11 : : std::vector<double>());
262 : const std::vector<double> adfMasks(
263 : bHasFlagMasks ? poFlagMasks->ReadAsDoubleArray()
264 11 : : std::vector<double>());
265 :
266 18 : if (bHasFlagValues &&
267 7 : adfValues.size() != static_cast<size_t>(aosFlagMeanings.size()))
268 : {
269 1 : CPLError(CE_Failure, CPLE_AppDefined,
270 : "Number of values in flag_values attribute is different "
271 : "from the one in flag_meanings");
272 1 : return false;
273 : }
274 :
275 16 : if (bHasFlagMasks &&
276 6 : adfMasks.size() != static_cast<size_t>(aosFlagMeanings.size()))
277 : {
278 1 : CPLError(CE_Failure, CPLE_AppDefined,
279 : "Number of values in flag_masks attribute is different "
280 : "from the one in flag_meanings");
281 1 : return false;
282 : }
283 :
284 19 : for (int i = 0; i < aosUnmaskFlags.size(); ++i)
285 : {
286 11 : const int nIdxFlag = aosFlagMeanings.FindString(aosUnmaskFlags[i]);
287 11 : if (nIdxFlag < 0)
288 : {
289 1 : CPLError(
290 : CE_Failure, CPLE_AppDefined,
291 : "Cannot fing flag %s in flag_meanings = '%s' attribute",
292 : aosUnmaskFlags[i], pszFlagMeanings);
293 1 : return false;
294 : }
295 :
296 10 : if (bHasFlagValues && adfValues[nIdxFlag] < 0)
297 : {
298 0 : CPLError(CE_Failure, CPLE_AppDefined,
299 : "Invalid value in flag_values[%d] = %f", nIdxFlag,
300 0 : adfValues[nIdxFlag]);
301 0 : return false;
302 : }
303 :
304 10 : if (bHasFlagMasks && adfMasks[nIdxFlag] < 0)
305 : {
306 0 : CPLError(CE_Failure, CPLE_AppDefined,
307 : "Invalid value in flag_masks[%d] = %f", nIdxFlag,
308 0 : adfMasks[nIdxFlag]);
309 0 : return false;
310 : }
311 :
312 10 : if (bHasFlagValues)
313 : {
314 12 : m_anValidFlagValues.push_back(
315 6 : static_cast<uint32_t>(adfValues[nIdxFlag]));
316 : }
317 :
318 10 : if (bHasFlagMasks)
319 : {
320 12 : m_anValidFlagMasks.push_back(
321 6 : static_cast<uint32_t>(adfMasks[nIdxFlag]));
322 : }
323 : }
324 : }
325 :
326 42 : return true;
327 : }
328 :
329 : /************************************************************************/
330 : /* IRead() */
331 : /************************************************************************/
332 :
333 51 : bool GDALMDArrayMask::IRead(const GUInt64 *arrayStartIdx, const size_t *count,
334 : const GInt64 *arrayStep,
335 : const GPtrDiff_t *bufferStride,
336 : const GDALExtendedDataType &bufferDataType,
337 : void *pDstBuffer) const
338 : {
339 51 : if (bufferDataType.GetClass() != GEDTC_NUMERIC)
340 : {
341 0 : CPLError(CE_Failure, CPLE_AppDefined,
342 : "%s: only reading to a numeric data type is supported",
343 : __func__);
344 0 : return false;
345 : }
346 51 : size_t nElts = 1;
347 51 : const size_t nDims = GetDimensionCount();
348 102 : std::vector<GPtrDiff_t> tmpBufferStrideVector(nDims);
349 139 : for (size_t i = 0; i < nDims; i++)
350 88 : nElts *= count[i];
351 51 : if (nDims > 0)
352 : {
353 46 : tmpBufferStrideVector.back() = 1;
354 88 : for (size_t i = nDims - 1; i > 0;)
355 : {
356 42 : --i;
357 42 : tmpBufferStrideVector[i] =
358 42 : tmpBufferStrideVector[i + 1] * count[i + 1];
359 : }
360 : }
361 :
362 : /* Optimized case: if we are an integer data type and that there is no */
363 : /* attribute that can be used to set mask = 0, then fill the mask buffer */
364 : /* directly */
365 49 : if (!m_bHasMissingValue && !m_bHasFillValue && !m_bHasValidMin &&
366 74 : !m_bHasValidMax && m_anValidFlagValues.empty() &&
367 34 : m_anValidFlagMasks.empty() &&
368 111 : m_poParent->GetRawNoDataValue() == nullptr &&
369 11 : GDALDataTypeIsInteger(m_poParent->GetDataType().GetNumericDataType()))
370 : {
371 7 : const bool bBufferDataTypeIsByte = bufferDataType == m_dt;
372 7 : if (bBufferDataTypeIsByte) // Byte case
373 : {
374 4 : bool bContiguous = true;
375 10 : for (size_t i = 0; i < nDims; i++)
376 : {
377 7 : if (bufferStride[i] != tmpBufferStrideVector[i])
378 : {
379 1 : bContiguous = false;
380 1 : break;
381 : }
382 : }
383 4 : if (bContiguous)
384 : {
385 : // CPLDebug("GDAL", "GetMask(): contiguous case");
386 3 : memset(pDstBuffer, 1, nElts);
387 3 : return true;
388 : }
389 : }
390 :
391 : struct Stack
392 : {
393 : size_t nIters = 0;
394 : GByte *dst_ptr = nullptr;
395 : GPtrDiff_t dst_inc_offset = 0;
396 : };
397 :
398 4 : std::vector<Stack> stack(std::max(static_cast<size_t>(1), nDims));
399 4 : const size_t nBufferDTSize = bufferDataType.GetSize();
400 13 : for (size_t i = 0; i < nDims; i++)
401 : {
402 9 : stack[i].dst_inc_offset =
403 9 : static_cast<GPtrDiff_t>(bufferStride[i] * nBufferDTSize);
404 : }
405 4 : stack[0].dst_ptr = static_cast<GByte *>(pDstBuffer);
406 :
407 4 : size_t dimIdx = 0;
408 4 : const size_t nDimsMinus1 = nDims > 0 ? nDims - 1 : 0;
409 : GByte abyOne[16]; // 16 is sizeof GDT_CFloat64
410 4 : CPLAssert(nBufferDTSize <= 16);
411 4 : const GByte flag = 1;
412 4 : GDALCopyWords64(&flag, GDT_UInt8, 0, abyOne,
413 : bufferDataType.GetNumericDataType(), 0, 1);
414 :
415 28 : lbl_next_depth:
416 28 : if (dimIdx == nDimsMinus1)
417 : {
418 19 : auto nIters = nDims > 0 ? count[dimIdx] : 1;
419 19 : GByte *dst_ptr = stack[dimIdx].dst_ptr;
420 :
421 : while (true)
422 : {
423 : // cppcheck-suppress knownConditionTrueFalse
424 73 : if (bBufferDataTypeIsByte)
425 : {
426 24 : *dst_ptr = flag;
427 : }
428 : else
429 : {
430 49 : memcpy(dst_ptr, abyOne, nBufferDTSize);
431 : }
432 :
433 73 : if ((--nIters) == 0)
434 19 : break;
435 54 : dst_ptr += stack[dimIdx].dst_inc_offset;
436 : }
437 : }
438 : else
439 : {
440 9 : stack[dimIdx].nIters = count[dimIdx];
441 : while (true)
442 : {
443 24 : dimIdx++;
444 24 : stack[dimIdx].dst_ptr = stack[dimIdx - 1].dst_ptr;
445 24 : goto lbl_next_depth;
446 24 : lbl_return_to_caller:
447 24 : dimIdx--;
448 24 : if ((--stack[dimIdx].nIters) == 0)
449 9 : break;
450 15 : stack[dimIdx].dst_ptr += stack[dimIdx].dst_inc_offset;
451 : }
452 : }
453 28 : if (dimIdx > 0)
454 24 : goto lbl_return_to_caller;
455 :
456 4 : return true;
457 : }
458 :
459 : const auto oTmpBufferDT =
460 44 : GDALDataTypeIsComplex(m_poParent->GetDataType().GetNumericDataType())
461 : ? GDALExtendedDataType::Create(GDT_Float64)
462 88 : : m_poParent->GetDataType();
463 44 : const size_t nTmpBufferDTSize = oTmpBufferDT.GetSize();
464 44 : void *pTempBuffer = VSI_MALLOC2_VERBOSE(nTmpBufferDTSize, nElts);
465 44 : if (!pTempBuffer)
466 0 : return false;
467 88 : if (!m_poParent->Read(arrayStartIdx, count, arrayStep,
468 44 : tmpBufferStrideVector.data(), oTmpBufferDT,
469 : pTempBuffer))
470 : {
471 0 : VSIFree(pTempBuffer);
472 0 : return false;
473 : }
474 :
475 44 : switch (oTmpBufferDT.GetNumericDataType())
476 : {
477 7 : case GDT_UInt8:
478 7 : ReadInternal<GByte>(count, bufferStride, bufferDataType, pDstBuffer,
479 : pTempBuffer, oTmpBufferDT,
480 : tmpBufferStrideVector);
481 7 : break;
482 :
483 0 : case GDT_Int8:
484 0 : ReadInternal<GInt8>(count, bufferStride, bufferDataType, pDstBuffer,
485 : pTempBuffer, oTmpBufferDT,
486 : tmpBufferStrideVector);
487 0 : break;
488 :
489 1 : case GDT_UInt16:
490 1 : ReadInternal<GUInt16>(count, bufferStride, bufferDataType,
491 : pDstBuffer, pTempBuffer, oTmpBufferDT,
492 : tmpBufferStrideVector);
493 1 : break;
494 :
495 14 : case GDT_Int16:
496 14 : ReadInternal<GInt16>(count, bufferStride, bufferDataType,
497 : pDstBuffer, pTempBuffer, oTmpBufferDT,
498 : tmpBufferStrideVector);
499 14 : break;
500 :
501 1 : case GDT_UInt32:
502 1 : ReadInternal<GUInt32>(count, bufferStride, bufferDataType,
503 : pDstBuffer, pTempBuffer, oTmpBufferDT,
504 : tmpBufferStrideVector);
505 1 : break;
506 :
507 5 : case GDT_Int32:
508 5 : ReadInternal<GInt32>(count, bufferStride, bufferDataType,
509 : pDstBuffer, pTempBuffer, oTmpBufferDT,
510 : tmpBufferStrideVector);
511 5 : break;
512 :
513 0 : case GDT_UInt64:
514 0 : ReadInternal<std::uint64_t>(count, bufferStride, bufferDataType,
515 : pDstBuffer, pTempBuffer, oTmpBufferDT,
516 : tmpBufferStrideVector);
517 0 : break;
518 :
519 0 : case GDT_Int64:
520 0 : ReadInternal<std::int64_t>(count, bufferStride, bufferDataType,
521 : pDstBuffer, pTempBuffer, oTmpBufferDT,
522 : tmpBufferStrideVector);
523 0 : break;
524 :
525 0 : case GDT_Float16:
526 0 : ReadInternal<GFloat16>(count, bufferStride, bufferDataType,
527 : pDstBuffer, pTempBuffer, oTmpBufferDT,
528 : tmpBufferStrideVector);
529 0 : break;
530 :
531 7 : case GDT_Float32:
532 7 : ReadInternal<float>(count, bufferStride, bufferDataType, pDstBuffer,
533 : pTempBuffer, oTmpBufferDT,
534 : tmpBufferStrideVector);
535 7 : break;
536 :
537 9 : case GDT_Float64:
538 9 : ReadInternal<double>(count, bufferStride, bufferDataType,
539 : pDstBuffer, pTempBuffer, oTmpBufferDT,
540 : tmpBufferStrideVector);
541 9 : break;
542 0 : case GDT_Unknown:
543 : case GDT_CInt16:
544 : case GDT_CInt32:
545 : case GDT_CFloat16:
546 : case GDT_CFloat32:
547 : case GDT_CFloat64:
548 : case GDT_TypeCount:
549 0 : CPLAssert(false);
550 : break;
551 : }
552 :
553 44 : VSIFree(pTempBuffer);
554 :
555 44 : return true;
556 : }
557 :
558 : /************************************************************************/
559 : /* IsValidForDT() */
560 : /************************************************************************/
561 :
562 40 : template <typename Type> static bool IsValidForDT(double dfVal)
563 : {
564 40 : if (std::isnan(dfVal))
565 0 : return false;
566 40 : if (dfVal < static_cast<double>(cpl::NumericLimits<Type>::lowest()))
567 0 : return false;
568 40 : if (dfVal > static_cast<double>(cpl::NumericLimits<Type>::max()))
569 0 : return false;
570 40 : return static_cast<double>(static_cast<Type>(dfVal)) == dfVal;
571 : }
572 :
573 9 : template <> bool IsValidForDT<double>(double)
574 : {
575 9 : return true;
576 : }
577 :
578 : /************************************************************************/
579 : /* IsNan() */
580 : /************************************************************************/
581 :
582 1438 : template <typename Type> inline bool IsNan(Type)
583 : {
584 1438 : return false;
585 : }
586 :
587 65 : template <> bool IsNan<double>(double val)
588 : {
589 65 : return std::isnan(val);
590 : }
591 :
592 26 : template <> bool IsNan<float>(float val)
593 : {
594 26 : return std::isnan(val);
595 : }
596 :
597 : /************************************************************************/
598 : /* ReadInternal() */
599 : /************************************************************************/
600 :
601 : template <typename Type>
602 44 : void GDALMDArrayMask::ReadInternal(
603 : const size_t *count, const GPtrDiff_t *bufferStride,
604 : const GDALExtendedDataType &bufferDataType, void *pDstBuffer,
605 : const void *pTempBuffer, const GDALExtendedDataType &oTmpBufferDT,
606 : const std::vector<GPtrDiff_t> &tmpBufferStrideVector) const
607 : {
608 44 : const size_t nDims = GetDimensionCount();
609 :
610 220 : const auto castValue = [](bool &bHasVal, double dfVal) -> Type
611 : {
612 220 : if (bHasVal)
613 : {
614 49 : if (IsValidForDT<Type>(dfVal))
615 : {
616 49 : return static_cast<Type>(dfVal);
617 : }
618 : else
619 : {
620 0 : bHasVal = false;
621 : }
622 : }
623 171 : return 0;
624 : };
625 :
626 44 : const void *pSrcRawNoDataValue = m_poParent->GetRawNoDataValue();
627 44 : bool bHasNodataValue = pSrcRawNoDataValue != nullptr;
628 : const Type nNoDataValue =
629 44 : castValue(bHasNodataValue, m_poParent->GetNoDataValueAsDouble());
630 44 : bool bHasMissingValue = m_bHasMissingValue;
631 44 : const Type nMissingValue = castValue(bHasMissingValue, m_dfMissingValue);
632 44 : bool bHasFillValue = m_bHasFillValue;
633 44 : const Type nFillValue = castValue(bHasFillValue, m_dfFillValue);
634 44 : bool bHasValidMin = m_bHasValidMin;
635 44 : const Type nValidMin = castValue(bHasValidMin, m_dfValidMin);
636 44 : bool bHasValidMax = m_bHasValidMax;
637 44 : const Type nValidMax = castValue(bHasValidMax, m_dfValidMax);
638 44 : const bool bHasValidFlags =
639 44 : !m_anValidFlagValues.empty() || !m_anValidFlagMasks.empty();
640 :
641 351 : const auto IsValidFlag = [this](Type v)
642 : {
643 54 : if (!m_anValidFlagValues.empty() && !m_anValidFlagMasks.empty())
644 : {
645 20 : for (size_t i = 0; i < m_anValidFlagValues.size(); ++i)
646 : {
647 12 : if ((static_cast<uint32_t>(v) & m_anValidFlagMasks[i]) ==
648 : m_anValidFlagValues[i])
649 : {
650 4 : return true;
651 : }
652 : }
653 : }
654 42 : else if (!m_anValidFlagValues.empty())
655 : {
656 49 : for (size_t i = 0; i < m_anValidFlagValues.size(); ++i)
657 : {
658 29 : if (static_cast<uint32_t>(v) == m_anValidFlagValues[i])
659 : {
660 4 : return true;
661 : }
662 : }
663 : }
664 : else /* if( !m_anValidFlagMasks.empty() ) */
665 : {
666 31 : for (size_t i = 0; i < m_anValidFlagMasks.size(); ++i)
667 : {
668 22 : if ((static_cast<uint32_t>(v) & m_anValidFlagMasks[i]) != 0)
669 : {
670 9 : return true;
671 : }
672 : }
673 : }
674 37 : return false;
675 : };
676 :
677 : #define GET_MASK_FOR_SAMPLE(v) \
678 : static_cast<GByte>(!IsNan(v) && !(bHasNodataValue && v == nNoDataValue) && \
679 : !(bHasMissingValue && v == nMissingValue) && \
680 : !(bHasFillValue && v == nFillValue) && \
681 : !(bHasValidMin && v < nValidMin) && \
682 : !(bHasValidMax && v > nValidMax) && \
683 : (!bHasValidFlags || IsValidFlag(v)));
684 :
685 44 : const bool bBufferDataTypeIsByte = bufferDataType == m_dt;
686 : /* Optimized case: Byte output and output buffer is contiguous */
687 44 : if (bBufferDataTypeIsByte)
688 : {
689 40 : bool bContiguous = true;
690 103 : for (size_t i = 0; i < nDims; i++)
691 : {
692 64 : if (bufferStride[i] != tmpBufferStrideVector[i])
693 : {
694 1 : bContiguous = false;
695 1 : break;
696 : }
697 : }
698 40 : if (bContiguous)
699 : {
700 39 : size_t nElts = 1;
701 102 : for (size_t i = 0; i < nDims; i++)
702 63 : nElts *= count[i];
703 :
704 1113 : for (size_t i = 0; i < nElts; i++)
705 : {
706 1074 : const Type *pSrc = static_cast<const Type *>(pTempBuffer) + i;
707 1074 : static_cast<GByte *>(pDstBuffer)[i] =
708 1074 : GET_MASK_FOR_SAMPLE(*pSrc);
709 : }
710 39 : return;
711 : }
712 : }
713 :
714 5 : const size_t nTmpBufferDTSize = oTmpBufferDT.GetSize();
715 :
716 : struct Stack
717 : {
718 : size_t nIters = 0;
719 : const GByte *src_ptr = nullptr;
720 : GByte *dst_ptr = nullptr;
721 : GPtrDiff_t src_inc_offset = 0;
722 : GPtrDiff_t dst_inc_offset = 0;
723 : };
724 :
725 10 : std::vector<Stack> stack(std::max(static_cast<size_t>(1), nDims));
726 5 : const size_t nBufferDTSize = bufferDataType.GetSize();
727 15 : for (size_t i = 0; i < nDims; i++)
728 : {
729 20 : stack[i].src_inc_offset = static_cast<GPtrDiff_t>(
730 10 : tmpBufferStrideVector[i] * nTmpBufferDTSize);
731 10 : stack[i].dst_inc_offset =
732 10 : static_cast<GPtrDiff_t>(bufferStride[i] * nBufferDTSize);
733 : }
734 5 : stack[0].src_ptr = static_cast<const GByte *>(pTempBuffer);
735 5 : stack[0].dst_ptr = static_cast<GByte *>(pDstBuffer);
736 :
737 5 : size_t dimIdx = 0;
738 5 : const size_t nDimsMinus1 = nDims > 0 ? nDims - 1 : 0;
739 : GByte abyZeroOrOne[2][16]; // 16 is sizeof GDT_CFloat64
740 5 : CPLAssert(nBufferDTSize <= 16);
741 15 : for (GByte flag = 0; flag <= 1; flag++)
742 : {
743 10 : GDALCopyWords64(&flag, m_dt.GetNumericDataType(), 0, abyZeroOrOne[flag],
744 : bufferDataType.GetNumericDataType(), 0, 1);
745 : }
746 :
747 43 : lbl_next_depth:
748 43 : if (dimIdx == nDimsMinus1)
749 : {
750 35 : auto nIters = nDims > 0 ? count[dimIdx] : 1;
751 35 : const GByte *src_ptr = stack[dimIdx].src_ptr;
752 35 : GByte *dst_ptr = stack[dimIdx].dst_ptr;
753 :
754 420 : while (true)
755 : {
756 455 : const Type *pSrc = reinterpret_cast<const Type *>(src_ptr);
757 455 : const GByte flag = GET_MASK_FOR_SAMPLE(*pSrc);
758 :
759 455 : if (bBufferDataTypeIsByte)
760 : {
761 24 : *dst_ptr = flag;
762 : }
763 : else
764 : {
765 431 : memcpy(dst_ptr, abyZeroOrOne[flag], nBufferDTSize);
766 : }
767 :
768 455 : if ((--nIters) == 0)
769 35 : break;
770 420 : src_ptr += stack[dimIdx].src_inc_offset;
771 420 : dst_ptr += stack[dimIdx].dst_inc_offset;
772 : }
773 : }
774 : else
775 : {
776 8 : stack[dimIdx].nIters = count[dimIdx];
777 : while (true)
778 : {
779 38 : dimIdx++;
780 38 : stack[dimIdx].src_ptr = stack[dimIdx - 1].src_ptr;
781 38 : stack[dimIdx].dst_ptr = stack[dimIdx - 1].dst_ptr;
782 38 : goto lbl_next_depth;
783 38 : lbl_return_to_caller:
784 38 : dimIdx--;
785 38 : if ((--stack[dimIdx].nIters) == 0)
786 8 : break;
787 30 : stack[dimIdx].src_ptr += stack[dimIdx].src_inc_offset;
788 30 : stack[dimIdx].dst_ptr += stack[dimIdx].dst_inc_offset;
789 : }
790 : }
791 43 : if (dimIdx > 0)
792 38 : goto lbl_return_to_caller;
793 : }
794 :
795 : //! @endcond
796 :
797 : /************************************************************************/
798 : /* GetMask() */
799 : /************************************************************************/
800 :
801 : /** Return an array that is a mask for the current array
802 :
803 : This array will be of type Byte, with values set to 0 to indicate invalid
804 : pixels of the current array, and values set to 1 to indicate valid pixels.
805 :
806 : The generic implementation honours the NoDataValue, as well as various
807 : netCDF CF attributes: missing_value, _FillValue, valid_min, valid_max
808 : and valid_range.
809 :
810 : Starting with GDAL 3.8, option UNMASK_FLAGS=flag_meaning_1[,flag_meaning_2,...]
811 : can be used to specify strings of the "flag_meanings" attribute
812 : (cf https://cfconventions.org/cf-conventions/cf-conventions.html#flags)
813 : for which pixels matching any of those flags will be set at 1 in the mask array,
814 : and pixels matching none of those flags will be set at 0.
815 : For example, let's consider the following netCDF variable defined with:
816 : \verbatim
817 : l2p_flags:valid_min = 0s ;
818 : l2p_flags:valid_max = 256s ;
819 : l2p_flags:flag_meanings = "microwave land ice lake river reserved_for_future_use unused_currently unused_currently unused_currently" ;
820 : l2p_flags:flag_masks = 1s, 2s, 4s, 8s, 16s, 32s, 64s, 128s, 256s ;
821 : \endverbatim
822 :
823 : GetMask(["UNMASK_FLAGS=microwave,land"]) will return an array such that:
824 : - for pixel values *outside* valid_range [0,256], the mask value will be 0.
825 : - for a pixel value with bit 0 or bit 1 at 1 within [0,256], the mask value
826 : will be 1.
827 : - for a pixel value with bit 0 and bit 1 at 0 within [0,256], the mask value
828 : will be 0.
829 :
830 : This is the same as the C function GDALMDArrayGetMask().
831 :
832 : @param papszOptions NULL-terminated list of options, or NULL.
833 :
834 : @return a new array, that holds a reference to the original one, and thus is
835 : a view of it (not a copy), or nullptr in case of error.
836 : */
837 : std::shared_ptr<GDALMDArray>
838 49 : GDALMDArray::GetMask(CSLConstList papszOptions) const
839 : {
840 98 : auto self = std::dynamic_pointer_cast<GDALMDArray>(m_pSelf.lock());
841 49 : if (!self)
842 : {
843 0 : CPLError(CE_Failure, CPLE_AppDefined,
844 : "Driver implementation issue: m_pSelf not set !");
845 0 : return nullptr;
846 : }
847 49 : if (GetDataType().GetClass() != GEDTC_NUMERIC)
848 : {
849 1 : CPLError(CE_Failure, CPLE_AppDefined,
850 : "GetMask() only supports numeric data type");
851 1 : return nullptr;
852 : }
853 48 : return GDALMDArrayMask::Create(self, papszOptions);
854 : }
|
