Line data Source code
1 : /******************************************************************************
2 : *
3 : * Name: gdalmultidim_array_unscaled.cpp
4 : * Project: GDAL Core
5 : * Purpose: GDALMDArray::GetUnscaled() 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 "gdal_multidim.h"
15 : #include "gdalmultidim_array_unscaled.h"
16 :
17 : #include <cmath>
18 :
19 : //! @cond Doxygen_Suppress
20 :
21 : /************************************************************************/
22 : /* IRead() */
23 : /************************************************************************/
24 :
25 16 : bool GDALMDArrayUnscaled::IRead(const GUInt64 *arrayStartIdx,
26 : const size_t *count, const GInt64 *arrayStep,
27 : const GPtrDiff_t *bufferStride,
28 : const GDALExtendedDataType &bufferDataType,
29 : void *pDstBuffer) const
30 : {
31 16 : const double dfScale = m_dfScale;
32 16 : const double dfOffset = m_dfOffset;
33 : const bool bDTIsComplex =
34 16 : CPL_TO_BOOL(GDALDataTypeIsComplex(m_dt.GetNumericDataType()));
35 : const auto dtDouble =
36 32 : GDALExtendedDataType::Create(bDTIsComplex ? GDT_CFloat64 : GDT_Float64);
37 16 : const size_t nDTSize = dtDouble.GetSize();
38 16 : const bool bTempBufferNeeded = (dtDouble != bufferDataType);
39 :
40 16 : double adfSrcNoData[2] = {0, 0};
41 16 : if (m_bHasNoData)
42 : {
43 9 : GDALExtendedDataType::CopyValue(m_poParent->GetRawNoDataValue(),
44 9 : m_poParent->GetDataType(),
45 : &adfSrcNoData[0], dtDouble);
46 : }
47 :
48 16 : const auto nDims = GetDimensions().size();
49 16 : if (nDims == 0)
50 : {
51 : double adfVal[2];
52 9 : if (!m_poParent->Read(arrayStartIdx, count, arrayStep, bufferStride,
53 : dtDouble, &adfVal[0]))
54 : {
55 0 : return false;
56 : }
57 9 : if (!m_bHasNoData || adfVal[0] != adfSrcNoData[0])
58 : {
59 6 : adfVal[0] = adfVal[0] * dfScale + dfOffset;
60 6 : if (bDTIsComplex)
61 : {
62 2 : adfVal[1] = adfVal[1] * dfScale + dfOffset;
63 : }
64 6 : GDALExtendedDataType::CopyValue(&adfVal[0], dtDouble, pDstBuffer,
65 : bufferDataType);
66 : }
67 : else
68 : {
69 3 : GDALExtendedDataType::CopyValue(m_abyRawNoData.data(), m_dt,
70 : pDstBuffer, bufferDataType);
71 : }
72 9 : return true;
73 : }
74 :
75 14 : std::vector<GPtrDiff_t> actualBufferStrideVector;
76 7 : const GPtrDiff_t *actualBufferStridePtr = bufferStride;
77 7 : void *pTempBuffer = pDstBuffer;
78 7 : if (bTempBufferNeeded)
79 : {
80 2 : size_t nElts = 1;
81 2 : actualBufferStrideVector.resize(nDims);
82 7 : for (size_t i = 0; i < nDims; i++)
83 5 : nElts *= count[i];
84 2 : actualBufferStrideVector.back() = 1;
85 5 : for (size_t i = nDims - 1; i > 0;)
86 : {
87 3 : --i;
88 3 : actualBufferStrideVector[i] =
89 3 : actualBufferStrideVector[i + 1] * count[i + 1];
90 : }
91 2 : actualBufferStridePtr = actualBufferStrideVector.data();
92 2 : pTempBuffer = VSI_MALLOC2_VERBOSE(nDTSize, nElts);
93 2 : if (!pTempBuffer)
94 0 : return false;
95 : }
96 7 : if (!m_poParent->Read(arrayStartIdx, count, arrayStep,
97 : actualBufferStridePtr, dtDouble, pTempBuffer))
98 : {
99 0 : if (bTempBufferNeeded)
100 0 : VSIFree(pTempBuffer);
101 0 : return false;
102 : }
103 :
104 : struct Stack
105 : {
106 : size_t nIters = 0;
107 : double *src_ptr = nullptr;
108 : GByte *dst_ptr = nullptr;
109 : GPtrDiff_t src_inc_offset = 0;
110 : GPtrDiff_t dst_inc_offset = 0;
111 : };
112 :
113 7 : std::vector<Stack> stack(nDims);
114 7 : const size_t nBufferDTSize = bufferDataType.GetSize();
115 23 : for (size_t i = 0; i < nDims; i++)
116 : {
117 32 : stack[i].src_inc_offset =
118 16 : actualBufferStridePtr[i] * (bDTIsComplex ? 2 : 1);
119 16 : stack[i].dst_inc_offset =
120 16 : static_cast<GPtrDiff_t>(bufferStride[i] * nBufferDTSize);
121 : }
122 7 : stack[0].src_ptr = static_cast<double *>(pTempBuffer);
123 7 : stack[0].dst_ptr = static_cast<GByte *>(pDstBuffer);
124 :
125 7 : size_t dimIdx = 0;
126 7 : const size_t nDimsMinus1 = nDims - 1;
127 : GByte abyDstNoData[16];
128 7 : CPLAssert(nBufferDTSize <= sizeof(abyDstNoData));
129 7 : GDALExtendedDataType::CopyValue(m_abyRawNoData.data(), m_dt, abyDstNoData,
130 : bufferDataType);
131 :
132 37 : lbl_next_depth:
133 37 : if (dimIdx == nDimsMinus1)
134 : {
135 25 : auto nIters = count[dimIdx];
136 25 : double *padfVal = stack[dimIdx].src_ptr;
137 25 : GByte *dst_ptr = stack[dimIdx].dst_ptr;
138 : while (true)
139 : {
140 92 : if (!m_bHasNoData || padfVal[0] != adfSrcNoData[0])
141 : {
142 88 : padfVal[0] = padfVal[0] * dfScale + dfOffset;
143 88 : if (bDTIsComplex)
144 : {
145 1 : padfVal[1] = padfVal[1] * dfScale + dfOffset;
146 : }
147 88 : if (bTempBufferNeeded)
148 : {
149 29 : GDALExtendedDataType::CopyValue(&padfVal[0], dtDouble,
150 : dst_ptr, bufferDataType);
151 : }
152 : }
153 : else
154 : {
155 4 : memcpy(dst_ptr, abyDstNoData, nBufferDTSize);
156 : }
157 :
158 92 : if ((--nIters) == 0)
159 25 : break;
160 67 : padfVal += stack[dimIdx].src_inc_offset;
161 67 : dst_ptr += stack[dimIdx].dst_inc_offset;
162 : }
163 : }
164 : else
165 : {
166 12 : stack[dimIdx].nIters = count[dimIdx];
167 : while (true)
168 : {
169 30 : dimIdx++;
170 30 : stack[dimIdx].src_ptr = stack[dimIdx - 1].src_ptr;
171 30 : stack[dimIdx].dst_ptr = stack[dimIdx - 1].dst_ptr;
172 30 : goto lbl_next_depth;
173 30 : lbl_return_to_caller:
174 30 : dimIdx--;
175 30 : if ((--stack[dimIdx].nIters) == 0)
176 12 : break;
177 18 : stack[dimIdx].src_ptr += stack[dimIdx].src_inc_offset;
178 18 : stack[dimIdx].dst_ptr += stack[dimIdx].dst_inc_offset;
179 : }
180 : }
181 37 : if (dimIdx > 0)
182 30 : goto lbl_return_to_caller;
183 :
184 7 : if (bTempBufferNeeded)
185 2 : VSIFree(pTempBuffer);
186 7 : return true;
187 : }
188 :
189 : /************************************************************************/
190 : /* IWrite() */
191 : /************************************************************************/
192 :
193 16 : bool GDALMDArrayUnscaled::IWrite(const GUInt64 *arrayStartIdx,
194 : const size_t *count, const GInt64 *arrayStep,
195 : const GPtrDiff_t *bufferStride,
196 : const GDALExtendedDataType &bufferDataType,
197 : const void *pSrcBuffer)
198 : {
199 16 : const double dfScale = m_dfScale;
200 16 : const double dfOffset = m_dfOffset;
201 : const bool bDTIsComplex =
202 16 : CPL_TO_BOOL(GDALDataTypeIsComplex(m_dt.GetNumericDataType()));
203 : const auto dtDouble =
204 32 : GDALExtendedDataType::Create(bDTIsComplex ? GDT_CFloat64 : GDT_Float64);
205 16 : const size_t nDTSize = dtDouble.GetSize();
206 16 : const bool bIsBufferDataTypeNativeDataType = (dtDouble == bufferDataType);
207 : const bool bSelfAndParentHaveNoData =
208 16 : m_bHasNoData && m_poParent->GetRawNoDataValue() != nullptr;
209 16 : double dfNoData = 0;
210 16 : if (m_bHasNoData)
211 : {
212 7 : GDALCopyWords64(m_abyRawNoData.data(), m_dt.GetNumericDataType(), 0,
213 : &dfNoData, GDT_Float64, 0, 1);
214 : }
215 :
216 16 : double adfSrcNoData[2] = {0, 0};
217 16 : if (bSelfAndParentHaveNoData)
218 : {
219 7 : GDALExtendedDataType::CopyValue(m_poParent->GetRawNoDataValue(),
220 7 : m_poParent->GetDataType(),
221 : &adfSrcNoData[0], dtDouble);
222 : }
223 :
224 16 : const auto nDims = GetDimensions().size();
225 16 : if (nDims == 0)
226 : {
227 : double adfVal[2];
228 10 : GDALExtendedDataType::CopyValue(pSrcBuffer, bufferDataType, &adfVal[0],
229 : dtDouble);
230 16 : if (bSelfAndParentHaveNoData &&
231 6 : (std::isnan(adfVal[0]) || adfVal[0] == dfNoData))
232 : {
233 4 : return m_poParent->Write(arrayStartIdx, count, arrayStep,
234 2 : bufferStride, m_poParent->GetDataType(),
235 4 : m_poParent->GetRawNoDataValue());
236 : }
237 : else
238 : {
239 8 : adfVal[0] = (adfVal[0] - dfOffset) / dfScale;
240 8 : if (bDTIsComplex)
241 : {
242 2 : adfVal[1] = (adfVal[1] - dfOffset) / dfScale;
243 : }
244 8 : return m_poParent->Write(arrayStartIdx, count, arrayStep,
245 8 : bufferStride, dtDouble, &adfVal[0]);
246 : }
247 : }
248 :
249 12 : std::vector<GPtrDiff_t> tmpBufferStrideVector;
250 6 : size_t nElts = 1;
251 6 : tmpBufferStrideVector.resize(nDims);
252 20 : for (size_t i = 0; i < nDims; i++)
253 14 : nElts *= count[i];
254 6 : tmpBufferStrideVector.back() = 1;
255 14 : for (size_t i = nDims - 1; i > 0;)
256 : {
257 8 : --i;
258 8 : tmpBufferStrideVector[i] = tmpBufferStrideVector[i + 1] * count[i + 1];
259 : }
260 6 : const GPtrDiff_t *tmpBufferStridePtr = tmpBufferStrideVector.data();
261 6 : void *pTempBuffer = VSI_MALLOC2_VERBOSE(nDTSize, nElts);
262 6 : if (!pTempBuffer)
263 0 : return false;
264 :
265 : struct Stack
266 : {
267 : size_t nIters = 0;
268 : double *dst_ptr = nullptr;
269 : const GByte *src_ptr = nullptr;
270 : GPtrDiff_t src_inc_offset = 0;
271 : GPtrDiff_t dst_inc_offset = 0;
272 : };
273 :
274 6 : std::vector<Stack> stack(nDims);
275 6 : const size_t nBufferDTSize = bufferDataType.GetSize();
276 20 : for (size_t i = 0; i < nDims; i++)
277 : {
278 28 : stack[i].dst_inc_offset =
279 14 : tmpBufferStridePtr[i] * (bDTIsComplex ? 2 : 1);
280 14 : stack[i].src_inc_offset =
281 14 : static_cast<GPtrDiff_t>(bufferStride[i] * nBufferDTSize);
282 : }
283 6 : stack[0].dst_ptr = static_cast<double *>(pTempBuffer);
284 6 : stack[0].src_ptr = static_cast<const GByte *>(pSrcBuffer);
285 :
286 6 : size_t dimIdx = 0;
287 6 : const size_t nDimsMinus1 = nDims - 1;
288 :
289 34 : lbl_next_depth:
290 34 : if (dimIdx == nDimsMinus1)
291 : {
292 23 : auto nIters = count[dimIdx];
293 23 : double *dst_ptr = stack[dimIdx].dst_ptr;
294 23 : const GByte *src_ptr = stack[dimIdx].src_ptr;
295 : while (true)
296 : {
297 : double adfVal[2];
298 : const double *padfSrcVal;
299 86 : if (bIsBufferDataTypeNativeDataType)
300 : {
301 50 : padfSrcVal = reinterpret_cast<const double *>(src_ptr);
302 : }
303 : else
304 : {
305 36 : GDALExtendedDataType::CopyValue(src_ptr, bufferDataType,
306 : &adfVal[0], dtDouble);
307 36 : padfSrcVal = adfVal;
308 : }
309 :
310 148 : if (bSelfAndParentHaveNoData &&
311 62 : (std::isnan(padfSrcVal[0]) || padfSrcVal[0] == dfNoData))
312 : {
313 3 : dst_ptr[0] = adfSrcNoData[0];
314 3 : if (bDTIsComplex)
315 : {
316 1 : dst_ptr[1] = adfSrcNoData[1];
317 : }
318 : }
319 : else
320 : {
321 83 : dst_ptr[0] = (padfSrcVal[0] - dfOffset) / dfScale;
322 83 : if (bDTIsComplex)
323 : {
324 1 : dst_ptr[1] = (padfSrcVal[1] - dfOffset) / dfScale;
325 : }
326 : }
327 :
328 86 : if ((--nIters) == 0)
329 23 : break;
330 63 : dst_ptr += stack[dimIdx].dst_inc_offset;
331 63 : src_ptr += stack[dimIdx].src_inc_offset;
332 63 : }
333 : }
334 : else
335 : {
336 11 : stack[dimIdx].nIters = count[dimIdx];
337 : while (true)
338 : {
339 28 : dimIdx++;
340 28 : stack[dimIdx].src_ptr = stack[dimIdx - 1].src_ptr;
341 28 : stack[dimIdx].dst_ptr = stack[dimIdx - 1].dst_ptr;
342 28 : goto lbl_next_depth;
343 28 : lbl_return_to_caller:
344 28 : dimIdx--;
345 28 : if ((--stack[dimIdx].nIters) == 0)
346 11 : break;
347 17 : stack[dimIdx].src_ptr += stack[dimIdx].src_inc_offset;
348 17 : stack[dimIdx].dst_ptr += stack[dimIdx].dst_inc_offset;
349 : }
350 : }
351 34 : if (dimIdx > 0)
352 28 : goto lbl_return_to_caller;
353 :
354 : // If the parent array is not double/complex-double, then convert the
355 : // values to it, before calling Write(), as some implementations can be
356 : // very slow when doing the type conversion.
357 6 : const auto &eParentDT = m_poParent->GetDataType();
358 6 : const size_t nParentDTSize = eParentDT.GetSize();
359 6 : if (nParentDTSize <= nDTSize / 2)
360 : {
361 : // Copy in-place by making sure that source and target do not overlap
362 6 : const auto eNumericDT = dtDouble.GetNumericDataType();
363 6 : const auto eParentNumericDT = eParentDT.GetNumericDataType();
364 :
365 : // Copy first element
366 : {
367 6 : std::vector<GByte> abyTemp(nParentDTSize);
368 6 : GDALCopyWords64(static_cast<GByte *>(pTempBuffer), eNumericDT,
369 6 : static_cast<int>(nDTSize), &abyTemp[0],
370 : eParentNumericDT, static_cast<int>(nParentDTSize),
371 : 1);
372 6 : memcpy(pTempBuffer, abyTemp.data(), abyTemp.size());
373 : }
374 : // Remaining elements
375 86 : for (size_t i = 1; i < nElts; ++i)
376 : {
377 80 : GDALCopyWords64(
378 80 : static_cast<GByte *>(pTempBuffer) + i * nDTSize, eNumericDT, 0,
379 80 : static_cast<GByte *>(pTempBuffer) + i * nParentDTSize,
380 : eParentNumericDT, 0, 1);
381 : }
382 : }
383 :
384 : const bool ret =
385 6 : m_poParent->Write(arrayStartIdx, count, arrayStep, tmpBufferStridePtr,
386 : eParentDT, pTempBuffer);
387 :
388 6 : VSIFree(pTempBuffer);
389 6 : return ret;
390 : }
391 :
392 : //! @endcond
393 :
394 : /************************************************************************/
395 : /* GetUnscaled() */
396 : /************************************************************************/
397 :
398 : /** Return an array that is the unscaled version of the current one.
399 : *
400 : * That is each value of the unscaled array will be
401 : * unscaled_value = raw_value * GetScale() + GetOffset()
402 : *
403 : * Starting with GDAL 3.3, the Write() method is implemented and will convert
404 : * from unscaled values to raw values.
405 : *
406 : * This is the same as the C function GDALMDArrayGetUnscaled().
407 : *
408 : * @param dfOverriddenScale Custom scale value instead of GetScale()
409 : * @param dfOverriddenOffset Custom offset value instead of GetOffset()
410 : * @param dfOverriddenDstNodata Custom target nodata value instead of NaN
411 : * @return a new array, that holds a reference to the original one, and thus is
412 : * a view of it (not a copy), or nullptr in case of error.
413 : */
414 : std::shared_ptr<GDALMDArray>
415 17 : GDALMDArray::GetUnscaled(double dfOverriddenScale, double dfOverriddenOffset,
416 : double dfOverriddenDstNodata) const
417 : {
418 34 : auto self = std::dynamic_pointer_cast<GDALMDArray>(m_pSelf.lock());
419 17 : if (!self)
420 : {
421 0 : CPLError(CE_Failure, CPLE_AppDefined,
422 : "Driver implementation issue: m_pSelf not set !");
423 0 : return nullptr;
424 : }
425 17 : if (GetDataType().GetClass() != GEDTC_NUMERIC)
426 : {
427 0 : CPLError(CE_Failure, CPLE_AppDefined,
428 : "GetUnscaled() only supports numeric data type");
429 0 : return nullptr;
430 : }
431 : const double dfScale =
432 17 : std::isnan(dfOverriddenScale) ? GetScale() : dfOverriddenScale;
433 : const double dfOffset =
434 17 : std::isnan(dfOverriddenOffset) ? GetOffset() : dfOverriddenOffset;
435 17 : if (dfScale == 1.0 && dfOffset == 0.0)
436 4 : return self;
437 :
438 13 : GDALDataType eDT = GDALDataTypeIsComplex(GetDataType().GetNumericDataType())
439 13 : ? GDT_CFloat64
440 13 : : GDT_Float64;
441 13 : if (dfOverriddenScale == -1 && dfOverriddenOffset == 0)
442 : {
443 1 : if (GetDataType().GetNumericDataType() == GDT_Float16)
444 0 : eDT = GDT_Float16;
445 1 : if (GetDataType().GetNumericDataType() == GDT_Float32)
446 1 : eDT = GDT_Float32;
447 : }
448 :
449 26 : return GDALMDArrayUnscaled::Create(self, dfScale, dfOffset,
450 13 : dfOverriddenDstNodata, eDT);
451 : }
|
