GDAL
gdalgeoloc_carray_accessor.h
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  * 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 
31 /************************************************************************/
32 /* GDALGeoLocCArrayAccessors */
33 /************************************************************************/
34 
35 class GDALGeoLocCArrayAccessors
36 {
37  typedef class GDALGeoLocCArrayAccessors AccessorType;
38 
39  GDALGeoLocTransformInfo* m_psTransform;
40  double *m_padfGeoLocX = nullptr;
41  double *m_padfGeoLocY = nullptr;
42  float *m_pafBackMapX = nullptr;
43  float *m_pafBackMapY = nullptr;
44  float *m_wgtsBackMap = nullptr;
45 
46  bool LoadGeoloc(bool bIsRegularGrid);
47 
48 public:
49 
50  template<class Type> struct CArrayAccessor
51  {
52  Type* m_array;
53  size_t m_nXSize;
54 
55  CArrayAccessor(Type* array, size_t nXSize): m_array(array), m_nXSize(nXSize) {}
56 
57  inline Type Get(int nX, int nY, bool* pbSuccess = nullptr)
58  {
59  if( pbSuccess )
60  *pbSuccess = true;
61  return m_array[nY * m_nXSize + nX];
62  }
63 
64  inline bool Set(int nX, int nY, Type val)
65  {
66  m_array[nY * m_nXSize + nX] = val;
67  return true;
68  }
69  };
70 
71  CArrayAccessor<double> geolocXAccessor;
72  CArrayAccessor<double> geolocYAccessor;
73  CArrayAccessor<float> backMapXAccessor;
74  CArrayAccessor<float> backMapYAccessor;
75  CArrayAccessor<float> backMapWeightAccessor;
76 
77  explicit GDALGeoLocCArrayAccessors(GDALGeoLocTransformInfo* psTransform):
78  m_psTransform(psTransform),
79  geolocXAccessor(nullptr, 0),
80  geolocYAccessor(nullptr, 0),
81  backMapXAccessor(nullptr, 0),
82  backMapYAccessor(nullptr, 0),
83  backMapWeightAccessor(nullptr, 0)
84  {
85  }
86 
87  ~GDALGeoLocCArrayAccessors()
88  {
89  VSIFree(m_pafBackMapX);
90  VSIFree(m_pafBackMapY);
91  VSIFree(m_padfGeoLocX);
92  VSIFree(m_padfGeoLocY);
93  VSIFree(m_wgtsBackMap);
94  }
95 
96  GDALGeoLocCArrayAccessors(const GDALGeoLocCArrayAccessors&) = delete;
97  GDALGeoLocCArrayAccessors& operator= (const GDALGeoLocCArrayAccessors&) = delete;
98 
99  bool Load(bool bIsRegularGrid, bool bUseQuadtree);
100 
101  bool AllocateBackMap();
102 
103  GDALDataset* GetBackmapDataset();
104  static void FlushBackmapCaches() {}
105  static void ReleaseBackmapDataset(GDALDataset* poDS) { delete poDS; }
106 
107  void FreeWghtsBackMap();
108 };
109 
110 /************************************************************************/
111 /* AllocateBackMap() */
112 /************************************************************************/
113 
114 bool GDALGeoLocCArrayAccessors::AllocateBackMap()
115 {
116  m_pafBackMapX = static_cast<float *>(
117  VSI_MALLOC3_VERBOSE(m_psTransform->nBackMapWidth,
118  m_psTransform->nBackMapHeight, sizeof(float)));
119  m_pafBackMapY = static_cast<float *>(
120  VSI_MALLOC3_VERBOSE(m_psTransform->nBackMapWidth,
121  m_psTransform->nBackMapHeight, sizeof(float)));
122 
123  m_wgtsBackMap = static_cast<float *>(
124  VSI_MALLOC3_VERBOSE(m_psTransform->nBackMapWidth,
125  m_psTransform->nBackMapHeight, sizeof(float)));
126 
127  if( m_pafBackMapX == nullptr ||
128  m_pafBackMapY == nullptr ||
129  m_wgtsBackMap == nullptr)
130  {
131  return false;
132  }
133 
134  const size_t nBMXYCount = static_cast<size_t>(m_psTransform->nBackMapWidth) *
135  m_psTransform->nBackMapHeight;
136  for( size_t i = 0; i < nBMXYCount; i++ )
137  {
138  m_pafBackMapX[i] = 0;
139  m_pafBackMapY[i] = 0;
140  m_wgtsBackMap[i] = 0.0;
141  }
142 
143  backMapXAccessor.m_array = m_pafBackMapX;
144  backMapXAccessor.m_nXSize = m_psTransform->nBackMapWidth;
145 
146  backMapYAccessor.m_array = m_pafBackMapY;
147  backMapYAccessor.m_nXSize = m_psTransform->nBackMapWidth;
148 
149  backMapWeightAccessor.m_array = m_wgtsBackMap;
150  backMapWeightAccessor.m_nXSize = m_psTransform->nBackMapWidth;
151 
152  return true;
153 }
154 
155 /************************************************************************/
156 /* FreeWghtsBackMap() */
157 /************************************************************************/
158 
159 void GDALGeoLocCArrayAccessors::FreeWghtsBackMap()
160 {
161  VSIFree(m_wgtsBackMap);
162  m_wgtsBackMap = nullptr;
163  backMapWeightAccessor.m_array = nullptr;
164  backMapWeightAccessor.m_nXSize = 0;
165 }
166 
167 /************************************************************************/
168 /* GetBackmapDataset() */
169 /************************************************************************/
170 
171 GDALDataset* GDALGeoLocCArrayAccessors::GetBackmapDataset()
172 {
173  auto poMEMDS = MEMDataset::Create( "",
174  m_psTransform->nBackMapWidth,
175  m_psTransform->nBackMapHeight,
176  0, GDT_Float32, nullptr );
177 
178  for( int i = 1; i <= 2; i++ )
179  {
180  char szBuffer[32] = { '\0' };
181  char szBuffer0[64] = { '\0' };
182  char* apszOptions[] = { szBuffer0, nullptr };
183 
184  void* ptr = (i == 1) ? m_pafBackMapX : m_pafBackMapY;
185  szBuffer[CPLPrintPointer(szBuffer, ptr, sizeof(szBuffer))] = '\0';
186  snprintf(szBuffer0, sizeof(szBuffer0), "DATAPOINTER=%s", szBuffer);
187  poMEMDS->AddBand(GDT_Float32, apszOptions);
188  poMEMDS->GetRasterBand(i)->SetNoDataValue(INVALID_BMXY);
189  }
190  return poMEMDS;
191 }
192 
193 /************************************************************************/
194 /* Load() */
195 /************************************************************************/
196 
197 bool GDALGeoLocCArrayAccessors::Load(bool bIsRegularGrid, bool bUseQuadtree)
198 {
199  return LoadGeoloc(bIsRegularGrid) &&
200  ((bUseQuadtree && GDALGeoLocBuildQuadTree(m_psTransform)) ||
201  (!bUseQuadtree && GDALGeoLoc<AccessorType>::GenerateBackMap(m_psTransform)));
202 }
203 
204 /************************************************************************/
205 /* LoadGeoloc() */
206 /************************************************************************/
207 
208 bool GDALGeoLocCArrayAccessors::LoadGeoloc(bool bIsRegularGrid)
209 
210 {
211  const int nXSize = m_psTransform->nGeoLocXSize;
212  const int nYSize = m_psTransform->nGeoLocYSize;
213 
214  m_padfGeoLocY = static_cast<double *>(
215  VSI_MALLOC3_VERBOSE(sizeof(double), nXSize, nYSize));
216  m_padfGeoLocX = static_cast<double *>(
217  VSI_MALLOC3_VERBOSE(sizeof(double), nXSize, nYSize));
218 
219  if( m_padfGeoLocX == nullptr ||
220  m_padfGeoLocY == nullptr )
221  {
222  return false;
223  }
224 
225  if( bIsRegularGrid )
226  {
227  // Case of regular grid.
228  // The XBAND contains the x coordinates for all lines.
229  // The YBAND contains the y coordinates for all columns.
230 
231  double* padfTempX = static_cast<double *>(
232  VSI_MALLOC2_VERBOSE(nXSize, sizeof(double)));
233  double* padfTempY = static_cast<double *>(
234  VSI_MALLOC2_VERBOSE(nYSize, sizeof(double)));
235  if( padfTempX == nullptr || padfTempY == nullptr )
236  {
237  CPLFree(padfTempX);
238  CPLFree(padfTempY);
239  return false;
240  }
241 
242  CPLErr eErr =
243  GDALRasterIO( m_psTransform->hBand_X, GF_Read,
244  0, 0, nXSize, 1,
245  padfTempX, nXSize, 1,
246  GDT_Float64, 0, 0 );
247 
248  for( size_t j = 0; j < static_cast<size_t>(nYSize); j++ )
249  {
250  memcpy( m_padfGeoLocX + j * nXSize,
251  padfTempX,
252  nXSize * sizeof(double) );
253  }
254 
255  if( eErr == CE_None )
256  {
257  eErr = GDALRasterIO( m_psTransform->hBand_Y, GF_Read,
258  0, 0, nYSize, 1,
259  padfTempY, nYSize, 1,
260  GDT_Float64, 0, 0 );
261 
262  for( size_t j = 0; j < static_cast<size_t>(nYSize); j++ )
263  {
264  for( size_t i = 0; i < static_cast<size_t>(nXSize); i++ )
265  {
266  m_padfGeoLocY[j * nXSize + i] = padfTempY[j];
267  }
268  }
269  }
270 
271  CPLFree(padfTempX);
272  CPLFree(padfTempY);
273 
274  if( eErr != CE_None )
275  return false;
276  }
277  else
278  {
279  if( GDALRasterIO( m_psTransform->hBand_X, GF_Read,
280  0, 0, nXSize, nYSize,
281  m_padfGeoLocX, nXSize, nYSize,
282  GDT_Float64, 0, 0 ) != CE_None
283  || GDALRasterIO( m_psTransform->hBand_Y, GF_Read,
284  0, 0, nXSize, nYSize,
285  m_padfGeoLocY, nXSize, nYSize,
286  GDT_Float64, 0, 0 ) != CE_None )
287  return false;
288  }
289 
290  geolocXAccessor.m_array = m_padfGeoLocX;
291  geolocXAccessor.m_nXSize = m_psTransform->nGeoLocXSize;
292 
293  geolocYAccessor.m_array = m_padfGeoLocY;
294  geolocYAccessor.m_nXSize = m_psTransform->nGeoLocXSize;
295 
296  return GDALGeoLoc<GDALGeoLocCArrayAccessors>::LoadGeolocFinish(m_psTransform);
297 }
298 
GDT_Float32
@ GDT_Float32
Definition: gdal.h:72
GF_Read
@ GF_Read
Definition: gdal.h:123
GDALDataset
A set of associated raster bands, usually from one file.
Definition: gdal_priv.h:339
VSI_MALLOC2_VERBOSE
#define VSI_MALLOC2_VERBOSE(nSize1, nSize2)
VSI_MALLOC2_VERBOSE.
Definition: cpl_vsi.h:292
GDT_Float64
@ GDT_Float64
Definition: gdal.h:73
VSIFree
void VSIFree(void *)
Analog of free() for data allocated with VSIMalloc(), VSICalloc(), VSIRealloc()
Definition: cpl_vsisimple.cpp:827
VSI_MALLOC3_VERBOSE
#define VSI_MALLOC3_VERBOSE(nSize1, nSize2, nSize3)
VSI_MALLOC3_VERBOSE.
Definition: cpl_vsi.h:297
CPLErr
CPLErr
Error category.
Definition: cpl_error.h:52
GDALRasterIO
CPLErr GDALRasterIO(GDALRasterBandH hRBand, GDALRWFlag eRWFlag, int nDSXOff, int nDSYOff, int nDSXSize, int nDSYSize, void *pBuffer, int nBXSize, int nBYSize, GDALDataType eBDataType, int nPixelSpace, int nLineSpace)
Read/write a region of image data for this band.
Definition: gdalrasterband.cpp:401
CPLPrintPointer
int CPLPrintPointer(char *, void *, int)
Print pointer value into specified string buffer.
Definition: cpl_conv.cpp:1386
CPLFree
#define CPLFree
Alias of VSIFree()
Definition: cpl_conv.h:83