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