Line data    Source code

       1             : /******************************************************************************
       2             :  * Name:     gdalmultidim_gltorthorectification.cpp
       3             :  * Project:  GDAL Core
       4             :  * Purpose:  GLT orthorectification implementation
       5             :  * Author:   Even Rouault <even.rouault at spatialys.com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2023, Even Rouault <even.rouault at spatialys.com>
       9             :  *
      10             :  * Permission is hereby granted, free of charge, to any person obtaining a
      11             :  * copy of this software and associated documentation files (the "Software"),
      12             :  * to deal in the Software without restriction, including without limitation
      13             :  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
      14             :  * and/or sell copies of the Software, and to permit persons to whom the
      15             :  * Software is furnished to do so, subject to the following conditions:
      16             :  *
      17             :  * The above copyright notice and this permission notice shall be included
      18             :  * in all copies or substantial portions of the Software.
      19             :  *
      20             :  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
      21             :  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
      22             :  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
      23             :  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
      24             :  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
      25             :  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
      26             :  * DEALINGS IN THE SOFTWARE.
      27             :  ****************************************************************************/
      28             : 
      29             : #include "gdal_priv.h"
      30             : #include "gdal_pam.h"
      31             : 
      32             : #include <algorithm>
      33             : #include <limits>
      34             : 
      35             : /************************************************************************/
      36             : /*                        GLTOrthoRectifiedArray                        */
      37             : /************************************************************************/
      38             : 
      39             : class GLTOrthoRectifiedArray final : public GDALPamMDArray
      40             : {
      41             :   private:
      42             :     std::shared_ptr<GDALMDArray> m_poParent{};
      43             :     std::vector<std::shared_ptr<GDALDimension>> m_apoDims;
      44             :     std::vector<GUInt64> m_anBlockSize;
      45             :     GDALExtendedDataType m_dt;
      46             :     std::shared_ptr<OGRSpatialReference> m_poSRS{};
      47             :     std::shared_ptr<GDALMDArray> m_poVarX{};
      48             :     std::shared_ptr<GDALMDArray> m_poVarY{};
      49             :     std::shared_ptr<GDALMDArray> m_poGLTX{};
      50             :     std::shared_ptr<GDALMDArray> m_poGLTY{};
      51             :     int m_nGLTIndexOffset = 0;
      52             : 
      53             :   protected:
      54           5 :     GLTOrthoRectifiedArray(
      55             :         const std::shared_ptr<GDALMDArray> &poParent,
      56             :         const std::vector<std::shared_ptr<GDALDimension>> &apoDims,
      57             :         const std::vector<GUInt64> &anBlockSize)
      58          10 :         : GDALAbstractMDArray(std::string(), "GLTOrthoRectifiedArray view of " +
      59           5 :                                                  poParent->GetFullName()),
      60          10 :           GDALPamMDArray(std::string(),
      61          10 :                          "GLTOrthoRectifiedArray view of " +
      62           5 :                              poParent->GetFullName(),
      63          10 :                          GDALPamMultiDim::GetPAM(poParent)),
      64           5 :           m_poParent(std::move(poParent)), m_apoDims(apoDims),
      65          25 :           m_anBlockSize(anBlockSize), m_dt(m_poParent->GetDataType())
      66             :     {
      67           5 :         CPLAssert(apoDims.size() == m_poParent->GetDimensionCount());
      68           5 :         CPLAssert(anBlockSize.size() == m_poParent->GetDimensionCount());
      69           5 :     }
      70             : 
      71             :     bool IRead(const GUInt64 *arrayStartIdx, const size_t *count,
      72             :                const GInt64 *arrayStep, const GPtrDiff_t *bufferStride,
      73             :                const GDALExtendedDataType &bufferDataType,
      74             :                void *pDstBuffer) const override;
      75             : 
      76             :   public:
      77             :     static std::shared_ptr<GDALMDArray>
      78           5 :     Create(const std::shared_ptr<GDALMDArray> &poParent,
      79             :            const std::shared_ptr<GDALMDArray> &poGLTX,
      80             :            const std::shared_ptr<GDALMDArray> &poGLTY, int nGLTIndexOffset,
      81             :            const std::vector<double> &adfGeoTransform)
      82             :     {
      83          10 :         std::vector<std::shared_ptr<GDALDimension>> apoNewDims;
      84             : 
      85             :         auto poDimY = std::make_shared<GDALDimensionWeakIndexingVar>(
      86           0 :             std::string(), "lat", GDAL_DIM_TYPE_HORIZONTAL_Y, "NORTH",
      87          10 :             poGLTX->GetDimensions()[0]->GetSize());
      88             :         auto varY = GDALMDArrayRegularlySpaced::Create(
      89          15 :             std::string(), poDimY->GetName(), poDimY,
      90          20 :             adfGeoTransform[3] + adfGeoTransform[5] / 2, adfGeoTransform[5], 0);
      91           5 :         poDimY->SetIndexingVariable(varY);
      92           5 :         apoNewDims.emplace_back(poDimY);
      93             : 
      94             :         auto poDimX = std::make_shared<GDALDimensionWeakIndexingVar>(
      95           0 :             std::string(), "lon", GDAL_DIM_TYPE_HORIZONTAL_X, "EAST",
      96          10 :             poGLTX->GetDimensions()[1]->GetSize());
      97             :         auto varX = GDALMDArrayRegularlySpaced::Create(
      98          15 :             std::string(), poDimX->GetName(), poDimX,
      99          20 :             adfGeoTransform[0] + adfGeoTransform[1] / 2, adfGeoTransform[1], 0);
     100           5 :         poDimX->SetIndexingVariable(varX);
     101           5 :         apoNewDims.emplace_back(poDimX);
     102             : 
     103           5 :         if (poParent->GetDimensionCount() == 3)
     104           4 :             apoNewDims.emplace_back(poParent->GetDimensions()[2]);
     105             : 
     106             :         std::vector<GUInt64> anBlockSize = std::vector<GUInt64>{
     107           5 :             std::min<GUInt64>(apoNewDims[0]->GetSize(), 512),
     108          10 :             std::min<GUInt64>(apoNewDims[1]->GetSize(), 512)};
     109           5 :         if (poParent->GetDimensionCount() == 3)
     110             :         {
     111           4 :             anBlockSize.push_back(poParent->GetDimensions()[2]->GetSize());
     112             :         }
     113             : 
     114             :         auto newAr(std::shared_ptr<GLTOrthoRectifiedArray>(
     115          10 :             new GLTOrthoRectifiedArray(poParent, apoNewDims, anBlockSize)));
     116           5 :         newAr->SetSelf(newAr);
     117           5 :         newAr->m_poVarX = varX;
     118           5 :         newAr->m_poVarY = varY;
     119           5 :         newAr->m_poGLTX = poGLTX;
     120           5 :         newAr->m_poGLTY = poGLTY;
     121           5 :         newAr->m_nGLTIndexOffset = nGLTIndexOffset;
     122          10 :         OGRSpatialReference oSRS;
     123           5 :         oSRS.importFromEPSG(4326);
     124           5 :         newAr->m_poSRS.reset(oSRS.Clone());
     125           5 :         newAr->m_poSRS->SetDataAxisToSRSAxisMapping(
     126             :             {/*latIdx = */ 1, /* lonIdx = */ 2});
     127          10 :         return newAr;
     128             :     }
     129             : 
     130           2 :     bool IsWritable() const override
     131             :     {
     132           2 :         return false;
     133             :     }
     134             : 
     135          14 :     const std::string &GetFilename() const override
     136             :     {
     137          14 :         return m_poParent->GetFilename();
     138             :     }
     139             : 
     140             :     const std::vector<std::shared_ptr<GDALDimension>> &
     141          83 :     GetDimensions() const override
     142             :     {
     143          83 :         return m_apoDims;
     144             :     }
     145             : 
     146          61 :     const GDALExtendedDataType &GetDataType() const override
     147             :     {
     148          61 :         return m_dt;
     149             :     }
     150             : 
     151           4 :     std::shared_ptr<OGRSpatialReference> GetSpatialRef() const override
     152             :     {
     153           4 :         return m_poSRS;
     154             :     }
     155             : 
     156           5 :     std::vector<GUInt64> GetBlockSize() const override
     157             :     {
     158           5 :         return m_anBlockSize;
     159             :     }
     160             : 
     161             :     std::shared_ptr<GDALAttribute>
     162           3 :     GetAttribute(const std::string &osName) const override
     163             :     {
     164           3 :         return m_poParent->GetAttribute(osName);
     165             :     }
     166             : 
     167             :     std::vector<std::shared_ptr<GDALAttribute>>
     168           2 :     GetAttributes(CSLConstList papszOptions = nullptr) const override
     169             :     {
     170           2 :         return m_poParent->GetAttributes(papszOptions);
     171             :     }
     172             : 
     173           3 :     const std::string &GetUnit() const override
     174             :     {
     175           3 :         return m_poParent->GetUnit();
     176             :     }
     177             : 
     178          15 :     const void *GetRawNoDataValue() const override
     179             :     {
     180          15 :         return m_poParent->GetRawNoDataValue();
     181             :     }
     182             : 
     183           0 :     double GetOffset(bool *pbHasOffset,
     184             :                      GDALDataType *peStorageType) const override
     185             :     {
     186           0 :         return m_poParent->GetOffset(pbHasOffset, peStorageType);
     187             :     }
     188             : 
     189           0 :     double GetScale(bool *pbHasScale,
     190             :                     GDALDataType *peStorageType) const override
     191             :     {
     192           0 :         return m_poParent->GetScale(pbHasScale, peStorageType);
     193             :     }
     194             : };
     195             : 
     196             : /************************************************************************/
     197             : /*                   GDALMDArrayResampled::IRead()                      */
     198             : /************************************************************************/
     199             : 
     200          12 : bool GLTOrthoRectifiedArray::IRead(const GUInt64 *arrayStartIdx,
     201             :                                    const size_t *count, const GInt64 *arrayStep,
     202             :                                    const GPtrDiff_t *bufferStride,
     203             :                                    const GDALExtendedDataType &bufferDataType,
     204             :                                    void *pDstBuffer) const
     205             : {
     206          12 :     if (bufferDataType.GetClass() != GEDTC_NUMERIC)
     207           0 :         return false;
     208             : 
     209          12 :     const size_t nXYValsCount = count[0] * count[1];
     210          24 :     const auto eInt32DT = GDALExtendedDataType::Create(GDT_Int32);
     211          24 :     std::vector<int32_t> anGLTX;
     212          24 :     std::vector<int32_t> anGLTY;
     213             :     try
     214             :     {
     215          12 :         anGLTX.resize(nXYValsCount);
     216          12 :         anGLTY.resize(nXYValsCount);
     217             :     }
     218           0 :     catch (const std::bad_alloc &e)
     219             :     {
     220           0 :         CPLError(CE_Failure, CPLE_OutOfMemory,
     221           0 :                  "GLTOrthoRectifiedArray::IRead(): %s", e.what());
     222           0 :         return false;
     223             :     }
     224          24 :     if (!m_poGLTX->Read(arrayStartIdx, count, arrayStep, nullptr, eInt32DT,
     225          36 :                         anGLTX.data()) ||
     226          24 :         !m_poGLTY->Read(arrayStartIdx, count, arrayStep, nullptr, eInt32DT,
     227          12 :                         anGLTY.data()))
     228             :     {
     229           0 :         return false;
     230             :     }
     231             : 
     232          12 :     int32_t nMinX = std::numeric_limits<int32_t>::max();
     233          12 :     int32_t nMaxX = std::numeric_limits<int32_t>::min();
     234          12 :     const auto nXSize = m_poParent->GetDimensions()[0]->GetSize();
     235          74 :     for (size_t i = 0; i < nXYValsCount; ++i)
     236             :     {
     237             :         const int64_t nX64 =
     238          62 :             static_cast<int64_t>(anGLTX[i]) + m_nGLTIndexOffset;
     239          62 :         if (nX64 >= 0 && static_cast<uint64_t>(nX64) < nXSize)
     240             :         {
     241          26 :             const int32_t nX = static_cast<int32_t>(nX64);
     242          26 :             if (nX < nMinX)
     243           8 :                 nMinX = nX;
     244          26 :             if (nX > nMaxX)
     245          14 :                 nMaxX = nX;
     246             :         }
     247             :     }
     248             : 
     249          12 :     int32_t nMinY = std::numeric_limits<int32_t>::max();
     250          12 :     int32_t nMaxY = std::numeric_limits<int32_t>::min();
     251          12 :     const auto nYSize = m_poParent->GetDimensions()[0]->GetSize();
     252          74 :     for (size_t i = 0; i < nXYValsCount; ++i)
     253             :     {
     254             :         const int64_t nY64 =
     255          62 :             static_cast<int64_t>(anGLTY[i]) + m_nGLTIndexOffset;
     256          62 :         if (nY64 >= 0 && static_cast<uint64_t>(nY64) < nYSize)
     257             :         {
     258          26 :             const int32_t nY = static_cast<int32_t>(nY64);
     259          26 :             if (nY < nMinY)
     260          14 :                 nMinY = nY;
     261          26 :             if (nY > nMaxY)
     262           8 :                 nMaxY = nY;
     263             :         }
     264             :     }
     265             : 
     266          12 :     const auto eBufferDT = bufferDataType.GetNumericDataType();
     267          12 :     auto pRawNoDataValue = GetRawNoDataValue();
     268          24 :     std::vector<GByte> abyNoData(16);
     269          12 :     if (pRawNoDataValue)
     270          12 :         GDALCopyWords(pRawNoDataValue, GetDataType().GetNumericDataType(), 0,
     271          12 :                       abyNoData.data(), eBufferDT, 0, 1);
     272             : 
     273          12 :     const auto nBufferDTSize = bufferDataType.GetSize();
     274             :     const int nCopyWordsDstStride =
     275          12 :         m_apoDims.size() == 3
     276          12 :             ? static_cast<int>(bufferStride[2] * nBufferDTSize)
     277          12 :             : 0;
     278             :     const int nCopyWordsCount =
     279          12 :         m_apoDims.size() == 3 ? static_cast<int>(count[2]) : 1;
     280          12 :     if (nMinX > nMaxX || nMinY > nMaxY)
     281             :     {
     282           8 :         for (size_t iY = 0; iY < count[0]; ++iY)
     283             :         {
     284          10 :             for (size_t iX = 0; iX < count[1]; ++iX)
     285             :             {
     286           6 :                 GByte *pabyDstBuffer =
     287             :                     static_cast<GByte *>(pDstBuffer) +
     288           6 :                     (iY * bufferStride[0] + iX * bufferStride[1]) *
     289           6 :                         static_cast<int>(nBufferDTSize);
     290           6 :                 GDALCopyWords(abyNoData.data(), eBufferDT, 0, pabyDstBuffer,
     291             :                               eBufferDT, nCopyWordsDstStride, nCopyWordsCount);
     292             :             }
     293             :         }
     294           4 :         return true;
     295             :     }
     296             : 
     297             :     GUInt64 parentArrayIdxStart[3] = {
     298           8 :         static_cast<GUInt64>(nMinY), static_cast<GUInt64>(nMinX),
     299           8 :         m_apoDims.size() == 3 ? arrayStartIdx[2] : 0};
     300           8 :     size_t parentCount[3] = {static_cast<size_t>(nMaxY - nMinY + 1),
     301           8 :                              static_cast<size_t>(nMaxX - nMinX + 1),
     302           8 :                              m_apoDims.size() == 3 ? count[2] : 1};
     303           8 :     GInt64 parentArrayStep[3] = {1, 1,
     304           8 :                                  m_apoDims.size() == 3 ? arrayStep[2] : 0};
     305             : 
     306           8 :     size_t nParentValueSize = nBufferDTSize;
     307          32 :     for (int i = 0; i < 3; ++i)
     308             :     {
     309          48 :         if (parentCount[i] >
     310          24 :             std::numeric_limits<size_t>::max() / nParentValueSize)
     311             :         {
     312           0 :             CPLError(
     313             :                 CE_Failure, CPLE_OutOfMemory,
     314             :                 "GLTOrthoRectifiedArray::IRead(): too big temporary array");
     315           0 :             return false;
     316             :         }
     317          24 :         nParentValueSize *= parentCount[i];
     318             :     }
     319             : 
     320           8 :     GPtrDiff_t parentStride[3] = {
     321           8 :         static_cast<GPtrDiff_t>(parentCount[1] * parentCount[2]),
     322           8 :         static_cast<GPtrDiff_t>(parentCount[2]), 1};
     323          16 :     std::vector<GByte> parentValues;
     324             :     try
     325             :     {
     326           8 :         parentValues.resize(nParentValueSize);
     327             :     }
     328           0 :     catch (const std::bad_alloc &e)
     329             :     {
     330           0 :         CPLError(CE_Failure, CPLE_OutOfMemory,
     331           0 :                  "GLTOrthoRectifiedArray::IRead(): %s", e.what());
     332           0 :         return false;
     333             :     }
     334          16 :     if (!m_poParent->Read(parentArrayIdxStart, parentCount, parentArrayStep,
     335           8 :                           parentStride, bufferDataType, parentValues.data()))
     336             :     {
     337           0 :         return false;
     338             :     }
     339             : 
     340           8 :     size_t iGLTIndex = 0;
     341           8 :     const size_t nXCount = parentCount[1];
     342           8 :     const size_t nBandCount = m_apoDims.size() == 3 ? parentCount[2] : 1;
     343          28 :     for (size_t iY = 0; iY < count[0]; ++iY)
     344             :     {
     345          76 :         for (size_t iX = 0; iX < count[1]; ++iX, ++iGLTIndex)
     346             :         {
     347             :             const int64_t nX64 =
     348          56 :                 static_cast<int64_t>(anGLTX[iGLTIndex]) + m_nGLTIndexOffset;
     349             :             const int64_t nY64 =
     350          56 :                 static_cast<int64_t>(anGLTY[iGLTIndex]) + m_nGLTIndexOffset;
     351          56 :             GByte *pabyDstBuffer =
     352             :                 static_cast<GByte *>(pDstBuffer) +
     353          56 :                 (iY * bufferStride[0] + iX * bufferStride[1]) *
     354          56 :                     static_cast<int>(nBufferDTSize);
     355          56 :             if (nX64 >= nMinX && nX64 <= nMaxX && nY64 >= nMinY &&
     356          26 :                 nY64 <= nMaxY)
     357             :             {
     358          26 :                 const int32_t iSrcX = static_cast<int32_t>(nX64) - nMinX;
     359          26 :                 const int32_t iSrcY = static_cast<int32_t>(nY64) - nMinY;
     360             :                 const GByte *pabySrcBuffer =
     361          26 :                     parentValues.data() +
     362          26 :                     (iSrcY * nXCount + iSrcX) * nBandCount * nBufferDTSize;
     363          26 :                 GDALCopyWords(pabySrcBuffer, eBufferDT,
     364             :                               static_cast<int>(nBufferDTSize), pabyDstBuffer,
     365          26 :                               eBufferDT, nCopyWordsDstStride, nCopyWordsCount);
     366             :             }
     367             :             else
     368             :             {
     369          30 :                 GDALCopyWords(abyNoData.data(), eBufferDT, 0, pabyDstBuffer,
     370             :                               eBufferDT, nCopyWordsDstStride, nCopyWordsCount);
     371             :             }
     372             :         }
     373             :     }
     374             : 
     375           8 :     return true;
     376             : }
     377             : 
     378             : /************************************************************************/
     379             : /*                     CreateGLTOrthorectified()                        */
     380             : /************************************************************************/
     381             : 
     382             : //! @cond Doxygen_Suppress
     383             : 
     384           5 : /* static */ std::shared_ptr<GDALMDArray> GDALMDArray::CreateGLTOrthorectified(
     385             :     const std::shared_ptr<GDALMDArray> &poParent,
     386             :     const std::shared_ptr<GDALMDArray> &poGLTX,
     387             :     const std::shared_ptr<GDALMDArray> &poGLTY, int nGLTIndexOffset,
     388             :     const std::vector<double> &adfGeoTransform)
     389             : {
     390             :     return GLTOrthoRectifiedArray::Create(poParent, poGLTX, poGLTY,
     391           5 :                                           nGLTIndexOffset, adfGeoTransform);
     392             : }
     393             : 
     394             : //! @endcond