Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL
       4             :  * Purpose:  C-Array storage of geolocation array and backmap
       5             :  * Author:   Even Rouault, <even.rouault at spatialys.com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2022, Planet Labs
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : /*! @cond Doxygen_Suppress */
      14             : 
      15             : /************************************************************************/
      16             : /*                        GDALGeoLocCArrayAccessors                     */
      17             : /************************************************************************/
      18             : 
      19             : class GDALGeoLocCArrayAccessors
      20             : {
      21             :     typedef class GDALGeoLocCArrayAccessors AccessorType;
      22             : 
      23             :     GDALGeoLocTransformInfo *m_psTransform;
      24             :     double *m_padfGeoLocX = nullptr;
      25             :     double *m_padfGeoLocY = nullptr;
      26             :     float *m_pafBackMapX = nullptr;
      27             :     float *m_pafBackMapY = nullptr;
      28             :     float *m_wgtsBackMap = nullptr;
      29             : 
      30             :     bool LoadGeoloc(bool bIsRegularGrid);
      31             : 
      32             :   public:
      33             :     template <class Type> struct CArrayAccessor
      34             :     {
      35             :         Type *m_array;
      36             :         size_t m_nXSize;
      37             : 
      38         230 :         CArrayAccessor(Type *array, size_t nXSize)
      39         230 :             : m_array(array), m_nXSize(nXSize)
      40             :         {
      41         230 :         }
      42             : 
      43    23334350 :         inline Type Get(int nX, int nY, bool *pbSuccess = nullptr)
      44             :         {
      45    23334350 :             if (pbSuccess)
      46           0 :                 *pbSuccess = true;
      47    23334350 :             return m_array[nY * m_nXSize + nX];
      48             :         }
      49             : 
      50     1932970 :         inline bool Set(int nX, int nY, Type val)
      51             :         {
      52     1932970 :             m_array[nY * m_nXSize + nX] = val;
      53     1932970 :             return true;
      54             :         }
      55             :     };
      56             : 
      57             :     CArrayAccessor<double> geolocXAccessor;
      58             :     CArrayAccessor<double> geolocYAccessor;
      59             :     CArrayAccessor<float> backMapXAccessor;
      60             :     CArrayAccessor<float> backMapYAccessor;
      61             :     CArrayAccessor<float> backMapWeightAccessor;
      62             : 
      63          46 :     explicit GDALGeoLocCArrayAccessors(GDALGeoLocTransformInfo *psTransform)
      64          46 :         : m_psTransform(psTransform), geolocXAccessor(nullptr, 0),
      65             :           geolocYAccessor(nullptr, 0), backMapXAccessor(nullptr, 0),
      66          46 :           backMapYAccessor(nullptr, 0), backMapWeightAccessor(nullptr, 0)
      67             :     {
      68          46 :     }
      69             : 
      70          46 :     ~GDALGeoLocCArrayAccessors()
      71          46 :     {
      72          46 :         VSIFree(m_pafBackMapX);
      73          46 :         VSIFree(m_pafBackMapY);
      74          46 :         VSIFree(m_padfGeoLocX);
      75          46 :         VSIFree(m_padfGeoLocY);
      76          46 :         VSIFree(m_wgtsBackMap);
      77          46 :     }
      78             : 
      79             :     GDALGeoLocCArrayAccessors(const GDALGeoLocCArrayAccessors &) = delete;
      80             :     GDALGeoLocCArrayAccessors &
      81             :     operator=(const GDALGeoLocCArrayAccessors &) = delete;
      82             : 
      83             :     bool Load(bool bIsRegularGrid, bool bUseQuadtree);
      84             : 
      85             :     bool AllocateBackMap();
      86             : 
      87             :     GDALDataset *GetBackmapDataset();
      88             : 
      89          42 :     static void FlushBackmapCaches()
      90             :     {
      91          42 :     }
      92             : 
      93          42 :     static void ReleaseBackmapDataset(GDALDataset *poDS)
      94             :     {
      95          42 :         delete poDS;
      96          42 :     }
      97             : 
      98             :     void FreeWghtsBackMap();
      99             : };
     100             : 
     101             : /************************************************************************/
     102             : /*                         AllocateBackMap()                            */
     103             : /************************************************************************/
     104             : 
     105          42 : bool GDALGeoLocCArrayAccessors::AllocateBackMap()
     106             : {
     107          42 :     m_pafBackMapX = static_cast<float *>(
     108          42 :         VSI_MALLOC3_VERBOSE(m_psTransform->nBackMapWidth,
     109             :                             m_psTransform->nBackMapHeight, sizeof(float)));
     110          42 :     m_pafBackMapY = static_cast<float *>(
     111          42 :         VSI_MALLOC3_VERBOSE(m_psTransform->nBackMapWidth,
     112             :                             m_psTransform->nBackMapHeight, sizeof(float)));
     113             : 
     114          42 :     m_wgtsBackMap = static_cast<float *>(
     115          42 :         VSI_MALLOC3_VERBOSE(m_psTransform->nBackMapWidth,
     116             :                             m_psTransform->nBackMapHeight, sizeof(float)));
     117             : 
     118          42 :     if (m_pafBackMapX == nullptr || m_pafBackMapY == nullptr ||
     119          42 :         m_wgtsBackMap == nullptr)
     120             :     {
     121           0 :         return false;
     122             :     }
     123             : 
     124          42 :     const size_t nBMXYCount =
     125          42 :         static_cast<size_t>(m_psTransform->nBackMapWidth) *
     126          42 :         m_psTransform->nBackMapHeight;
     127      186014 :     for (size_t i = 0; i < nBMXYCount; i++)
     128             :     {
     129      185972 :         m_pafBackMapX[i] = 0;
     130      185972 :         m_pafBackMapY[i] = 0;
     131      185972 :         m_wgtsBackMap[i] = 0.0;
     132             :     }
     133             : 
     134          42 :     backMapXAccessor.m_array = m_pafBackMapX;
     135          42 :     backMapXAccessor.m_nXSize = m_psTransform->nBackMapWidth;
     136             : 
     137          42 :     backMapYAccessor.m_array = m_pafBackMapY;
     138          42 :     backMapYAccessor.m_nXSize = m_psTransform->nBackMapWidth;
     139             : 
     140          42 :     backMapWeightAccessor.m_array = m_wgtsBackMap;
     141          42 :     backMapWeightAccessor.m_nXSize = m_psTransform->nBackMapWidth;
     142             : 
     143          42 :     return true;
     144             : }
     145             : 
     146             : /************************************************************************/
     147             : /*                         FreeWghtsBackMap()                           */
     148             : /************************************************************************/
     149             : 
     150          42 : void GDALGeoLocCArrayAccessors::FreeWghtsBackMap()
     151             : {
     152          42 :     VSIFree(m_wgtsBackMap);
     153          42 :     m_wgtsBackMap = nullptr;
     154          42 :     backMapWeightAccessor.m_array = nullptr;
     155          42 :     backMapWeightAccessor.m_nXSize = 0;
     156          42 : }
     157             : 
     158             : /************************************************************************/
     159             : /*                        GetBackmapDataset()                           */
     160             : /************************************************************************/
     161             : 
     162          42 : GDALDataset *GDALGeoLocCArrayAccessors::GetBackmapDataset()
     163             : {
     164          84 :     auto poMEMDS = MEMDataset::Create("", m_psTransform->nBackMapWidth,
     165          42 :                                       m_psTransform->nBackMapHeight, 0,
     166             :                                       GDT_Float32, nullptr);
     167             : 
     168         126 :     for (int i = 1; i <= 2; i++)
     169             :     {
     170          84 :         void *ptr = (i == 1) ? m_pafBackMapX : m_pafBackMapY;
     171          84 :         GDALRasterBandH hMEMBand = MEMCreateRasterBandEx(
     172             :             poMEMDS, i, static_cast<GByte *>(ptr), GDT_Float32, 0, 0, false);
     173          84 :         poMEMDS->AddMEMBand(hMEMBand);
     174          84 :         poMEMDS->GetRasterBand(i)->SetNoDataValue(INVALID_BMXY);
     175             :     }
     176          42 :     return poMEMDS;
     177             : }
     178             : 
     179             : /************************************************************************/
     180             : /*                             Load()                                   */
     181             : /************************************************************************/
     182             : 
     183          46 : bool GDALGeoLocCArrayAccessors::Load(bool bIsRegularGrid, bool bUseQuadtree)
     184             : {
     185          92 :     return LoadGeoloc(bIsRegularGrid) &&
     186           4 :            ((bUseQuadtree && GDALGeoLocBuildQuadTree(m_psTransform)) ||
     187          42 :             (!bUseQuadtree &&
     188          88 :              GDALGeoLoc<AccessorType>::GenerateBackMap(m_psTransform)));
     189             : }
     190             : 
     191             : /************************************************************************/
     192             : /*                          LoadGeoloc()                                */
     193             : /************************************************************************/
     194             : 
     195          46 : bool GDALGeoLocCArrayAccessors::LoadGeoloc(bool bIsRegularGrid)
     196             : 
     197             : {
     198          46 :     const int nXSize = m_psTransform->nGeoLocXSize;
     199          46 :     const int nYSize = m_psTransform->nGeoLocYSize;
     200             : 
     201          46 :     m_padfGeoLocY = static_cast<double *>(
     202          46 :         VSI_MALLOC3_VERBOSE(sizeof(double), nXSize, nYSize));
     203          46 :     m_padfGeoLocX = static_cast<double *>(
     204          46 :         VSI_MALLOC3_VERBOSE(sizeof(double), nXSize, nYSize));
     205             : 
     206          46 :     if (m_padfGeoLocX == nullptr || m_padfGeoLocY == nullptr)
     207             :     {
     208           0 :         return false;
     209             :     }
     210             : 
     211          46 :     if (bIsRegularGrid)
     212             :     {
     213             :         // Case of regular grid.
     214             :         // The XBAND contains the x coordinates for all lines.
     215             :         // The YBAND contains the y coordinates for all columns.
     216             : 
     217             :         double *padfTempX =
     218          14 :             static_cast<double *>(VSI_MALLOC2_VERBOSE(nXSize, sizeof(double)));
     219             :         double *padfTempY =
     220          14 :             static_cast<double *>(VSI_MALLOC2_VERBOSE(nYSize, sizeof(double)));
     221          14 :         if (padfTempX == nullptr || padfTempY == nullptr)
     222             :         {
     223           0 :             CPLFree(padfTempX);
     224           0 :             CPLFree(padfTempY);
     225           0 :             return false;
     226             :         }
     227             : 
     228             :         CPLErr eErr =
     229          14 :             GDALRasterIO(m_psTransform->hBand_X, GF_Read, 0, 0, nXSize, 1,
     230             :                          padfTempX, nXSize, 1, GDT_Float64, 0, 0);
     231             : 
     232         262 :         for (size_t j = 0; j < static_cast<size_t>(nYSize); j++)
     233             :         {
     234         248 :             memcpy(m_padfGeoLocX + j * nXSize, padfTempX,
     235         248 :                    nXSize * sizeof(double));
     236             :         }
     237             : 
     238          14 :         if (eErr == CE_None)
     239             :         {
     240          14 :             eErr = GDALRasterIO(m_psTransform->hBand_Y, GF_Read, 0, 0, nYSize,
     241             :                                 1, padfTempY, nYSize, 1, GDT_Float64, 0, 0);
     242             : 
     243         262 :             for (size_t j = 0; j < static_cast<size_t>(nYSize); j++)
     244             :             {
     245       31480 :                 for (size_t i = 0; i < static_cast<size_t>(nXSize); i++)
     246             :                 {
     247       31232 :                     m_padfGeoLocY[j * nXSize + i] = padfTempY[j];
     248             :                 }
     249             :             }
     250             :         }
     251             : 
     252          14 :         CPLFree(padfTempX);
     253          14 :         CPLFree(padfTempY);
     254             : 
     255          14 :         if (eErr != CE_None)
     256           0 :             return false;
     257             :     }
     258             :     else
     259             :     {
     260          96 :         if (GDALRasterIO(m_psTransform->hBand_X, GF_Read, 0, 0, nXSize, nYSize,
     261          32 :                          m_padfGeoLocX, nXSize, nYSize, GDT_Float64, 0,
     262          64 :                          0) != CE_None ||
     263          32 :             GDALRasterIO(m_psTransform->hBand_Y, GF_Read, 0, 0, nXSize, nYSize,
     264          32 :                          m_padfGeoLocY, nXSize, nYSize, GDT_Float64, 0,
     265             :                          0) != CE_None)
     266           0 :             return false;
     267             :     }
     268             : 
     269          46 :     geolocXAccessor.m_array = m_padfGeoLocX;
     270          46 :     geolocXAccessor.m_nXSize = m_psTransform->nGeoLocXSize;
     271             : 
     272          46 :     geolocYAccessor.m_array = m_padfGeoLocY;
     273          46 :     geolocYAccessor.m_nXSize = m_psTransform->nGeoLocXSize;
     274             : 
     275          46 :     GDALGeoLoc<GDALGeoLocCArrayAccessors>::LoadGeolocFinish(m_psTransform);
     276          46 :     return true;
     277             : }
     278             : 
     279             : /*! @endcond */