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  template <class Type> struct CArrayAccessor
50  {
51  Type *m_array;
52  size_t m_nXSize;
53 
54  CArrayAccessor(Type *array, size_t nXSize)
55  : m_array(array), m_nXSize(nXSize)
56  {
57  }
58 
59  inline Type Get(int nX, int nY, bool *pbSuccess = nullptr)
60  {
61  if (pbSuccess)
62  *pbSuccess = true;
63  return m_array[nY * m_nXSize + nX];
64  }
65 
66  inline bool Set(int nX, int nY, Type val)
67  {
68  m_array[nY * m_nXSize + nX] = val;
69  return true;
70  }
71  };
72 
73  CArrayAccessor<double> geolocXAccessor;
74  CArrayAccessor<double> geolocYAccessor;
75  CArrayAccessor<float> backMapXAccessor;
76  CArrayAccessor<float> backMapYAccessor;
77  CArrayAccessor<float> backMapWeightAccessor;
78 
79  explicit GDALGeoLocCArrayAccessors(GDALGeoLocTransformInfo *psTransform)
80  : m_psTransform(psTransform), geolocXAccessor(nullptr, 0),
81  geolocYAccessor(nullptr, 0), backMapXAccessor(nullptr, 0),
82  backMapYAccessor(nullptr, 0), backMapWeightAccessor(nullptr, 0)
83  {
84  }
85 
86  ~GDALGeoLocCArrayAccessors()
87  {
88  VSIFree(m_pafBackMapX);
89  VSIFree(m_pafBackMapY);
90  VSIFree(m_padfGeoLocX);
91  VSIFree(m_padfGeoLocY);
92  VSIFree(m_wgtsBackMap);
93  }
94 
95  GDALGeoLocCArrayAccessors(const GDALGeoLocCArrayAccessors &) = delete;
96  GDALGeoLocCArrayAccessors &
97  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  {
106  }
107  static void ReleaseBackmapDataset(GDALDataset *poDS)
108  {
109  delete poDS;
110  }
111 
112  void FreeWghtsBackMap();
113 };
114 
115 /************************************************************************/
116 /* AllocateBackMap() */
117 /************************************************************************/
118 
119 bool GDALGeoLocCArrayAccessors::AllocateBackMap()
120 {
121  m_pafBackMapX = static_cast<float *>(
122  VSI_MALLOC3_VERBOSE(m_psTransform->nBackMapWidth,
123  m_psTransform->nBackMapHeight, sizeof(float)));
124  m_pafBackMapY = static_cast<float *>(
125  VSI_MALLOC3_VERBOSE(m_psTransform->nBackMapWidth,
126  m_psTransform->nBackMapHeight, sizeof(float)));
127 
128  m_wgtsBackMap = static_cast<float *>(
129  VSI_MALLOC3_VERBOSE(m_psTransform->nBackMapWidth,
130  m_psTransform->nBackMapHeight, sizeof(float)));
131 
132  if (m_pafBackMapX == nullptr || m_pafBackMapY == nullptr ||
133  m_wgtsBackMap == nullptr)
134  {
135  return false;
136  }
137 
138  const size_t nBMXYCount =
139  static_cast<size_t>(m_psTransform->nBackMapWidth) *
140  m_psTransform->nBackMapHeight;
141  for (size_t i = 0; i < nBMXYCount; i++)
142  {
143  m_pafBackMapX[i] = 0;
144  m_pafBackMapY[i] = 0;
145  m_wgtsBackMap[i] = 0.0;
146  }
147 
148  backMapXAccessor.m_array = m_pafBackMapX;
149  backMapXAccessor.m_nXSize = m_psTransform->nBackMapWidth;
150 
151  backMapYAccessor.m_array = m_pafBackMapY;
152  backMapYAccessor.m_nXSize = m_psTransform->nBackMapWidth;
153 
154  backMapWeightAccessor.m_array = m_wgtsBackMap;
155  backMapWeightAccessor.m_nXSize = m_psTransform->nBackMapWidth;
156 
157  return true;
158 }
159 
160 /************************************************************************/
161 /* FreeWghtsBackMap() */
162 /************************************************************************/
163 
164 void GDALGeoLocCArrayAccessors::FreeWghtsBackMap()
165 {
166  VSIFree(m_wgtsBackMap);
167  m_wgtsBackMap = nullptr;
168  backMapWeightAccessor.m_array = nullptr;
169  backMapWeightAccessor.m_nXSize = 0;
170 }
171 
172 /************************************************************************/
173 /* GetBackmapDataset() */
174 /************************************************************************/
175 
176 GDALDataset *GDALGeoLocCArrayAccessors::GetBackmapDataset()
177 {
178  auto poMEMDS = MEMDataset::Create("", m_psTransform->nBackMapWidth,
179  m_psTransform->nBackMapHeight, 0,
180  GDT_Float32, nullptr);
181 
182  for (int i = 1; i <= 2; i++)
183  {
184  void *ptr = (i == 1) ? m_pafBackMapX : m_pafBackMapY;
185  GDALRasterBandH hMEMBand = MEMCreateRasterBandEx(
186  poMEMDS, i, static_cast<GByte *>(ptr), GDT_Float32, 0, 0, false);
187  poMEMDS->AddMEMBand(hMEMBand);
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 &&
202  GDALGeoLoc<AccessorType>::GenerateBackMap(m_psTransform)));
203 }
204 
205 /************************************************************************/
206 /* LoadGeoloc() */
207 /************************************************************************/
208 
209 bool GDALGeoLocCArrayAccessors::LoadGeoloc(bool bIsRegularGrid)
210 
211 {
212  const int nXSize = m_psTransform->nGeoLocXSize;
213  const int nYSize = m_psTransform->nGeoLocYSize;
214 
215  m_padfGeoLocY = static_cast<double *>(
216  VSI_MALLOC3_VERBOSE(sizeof(double), nXSize, nYSize));
217  m_padfGeoLocX = static_cast<double *>(
218  VSI_MALLOC3_VERBOSE(sizeof(double), nXSize, nYSize));
219 
220  if (m_padfGeoLocX == nullptr || 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 =
232  static_cast<double *>(VSI_MALLOC2_VERBOSE(nXSize, sizeof(double)));
233  double *padfTempY =
234  static_cast<double *>(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, 0, 0, nXSize, 1,
244  padfTempX, nXSize, 1, GDT_Float64, 0, 0);
245 
246  for (size_t j = 0; j < static_cast<size_t>(nYSize); j++)
247  {
248  memcpy(m_padfGeoLocX + j * nXSize, padfTempX,
249  nXSize * sizeof(double));
250  }
251 
252  if (eErr == CE_None)
253  {
254  eErr = GDALRasterIO(m_psTransform->hBand_Y, GF_Read, 0, 0, nYSize,
255  1, padfTempY, nYSize, 1, GDT_Float64, 0, 0);
256 
257  for (size_t j = 0; j < static_cast<size_t>(nYSize); j++)
258  {
259  for (size_t i = 0; i < static_cast<size_t>(nXSize); i++)
260  {
261  m_padfGeoLocY[j * nXSize + i] = padfTempY[j];
262  }
263  }
264  }
265 
266  CPLFree(padfTempX);
267  CPLFree(padfTempY);
268 
269  if (eErr != CE_None)
270  return false;
271  }
272  else
273  {
274  if (GDALRasterIO(m_psTransform->hBand_X, GF_Read, 0, 0, nXSize, nYSize,
275  m_padfGeoLocX, nXSize, nYSize, GDT_Float64, 0,
276  0) != CE_None ||
277  GDALRasterIO(m_psTransform->hBand_Y, GF_Read, 0, 0, nXSize, nYSize,
278  m_padfGeoLocY, nXSize, nYSize, GDT_Float64, 0,
279  0) != CE_None)
280  return false;
281  }
282 
283  geolocXAccessor.m_array = m_padfGeoLocX;
284  geolocXAccessor.m_nXSize = m_psTransform->nGeoLocXSize;
285 
286  geolocYAccessor.m_array = m_padfGeoLocY;
287  geolocYAccessor.m_nXSize = m_psTransform->nGeoLocXSize;
288 
289  GDALGeoLoc<GDALGeoLocCArrayAccessors>::LoadGeolocFinish(m_psTransform);
290  return true;
291 }
292 
GByte
unsigned char GByte
Unsigned byte type.
Definition: cpl_port.h:196
GDT_Float32
@ GDT_Float32
Definition: gdal.h:74
GF_Read
@ GF_Read
Definition: gdal.h:133
GDALDataset
A set of associated raster bands, usually from one file.
Definition: gdal_priv.h:348
VSI_MALLOC2_VERBOSE
#define VSI_MALLOC2_VERBOSE(nSize1, nSize2)
VSI_MALLOC2_VERBOSE.
Definition: cpl_vsi.h:337
GDT_Float64
@ GDT_Float64
Definition: gdal.h:75
VSIFree
void VSIFree(void *)
Analog of free() for data allocated with VSIMalloc(), VSICalloc(), VSIRealloc()
Definition: cpl_vsisimple.cpp:843
VSI_MALLOC3_VERBOSE
#define VSI_MALLOC3_VERBOSE(nSize1, nSize2, nSize3)
VSI_MALLOC3_VERBOSE.
Definition: cpl_vsi.h:345
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:435
GDALRasterBandH
void * GDALRasterBandH
Opaque type used for the C bindings of the C++ GDALRasterBand class.
Definition: gdal.h:294
CPLFree
#define CPLFree
Alias of VSIFree()
Definition: cpl_conv.h:98