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 * SPDX-License-Identifier: MIT
11 ****************************************************************************/
12
15/************************************************************************/
16/* GDALGeoLocCArrayAccessors */
17/************************************************************************/
18
19class 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 CArrayAccessor(Type *array, size_t nXSize)
39 : m_array(array), m_nXSize(nXSize)
40 {
41 }
42
43 inline Type Get(int nX, int nY, bool *pbSuccess = nullptr)
44 {
45 if (pbSuccess)
46 *pbSuccess = true;
47 return m_array[nY * m_nXSize + nX];
48 }
49
50 inline bool Set(int nX, int nY, Type val)
51 {
52 m_array[nY * m_nXSize + nX] = val;
53 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 explicit GDALGeoLocCArrayAccessors(GDALGeoLocTransformInfo *psTransform)
64 : m_psTransform(psTransform), geolocXAccessor(nullptr, 0),
65 geolocYAccessor(nullptr, 0), backMapXAccessor(nullptr, 0),
66 backMapYAccessor(nullptr, 0), backMapWeightAccessor(nullptr, 0)
67 {
68 }
69
70 ~GDALGeoLocCArrayAccessors()
71 {
72 VSIFree(m_pafBackMapX);
73 VSIFree(m_pafBackMapY);
74 VSIFree(m_padfGeoLocX);
75 VSIFree(m_padfGeoLocY);
76 VSIFree(m_wgtsBackMap);
77 }
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 static void FlushBackmapCaches()
90 {
91 }
92
93 static void ReleaseBackmapDataset(GDALDataset *poDS)
94 {
95 delete poDS;
96 }
97
98 void FreeWghtsBackMap();
99};
100
101/************************************************************************/
102/* AllocateBackMap() */
103/************************************************************************/
104
105bool GDALGeoLocCArrayAccessors::AllocateBackMap()
106{
107 m_pafBackMapX = static_cast<float *>(
108 VSI_MALLOC3_VERBOSE(m_psTransform->nBackMapWidth,
109 m_psTransform->nBackMapHeight, sizeof(float)));
110 m_pafBackMapY = static_cast<float *>(
111 VSI_MALLOC3_VERBOSE(m_psTransform->nBackMapWidth,
112 m_psTransform->nBackMapHeight, sizeof(float)));
113
114 m_wgtsBackMap = static_cast<float *>(
115 VSI_MALLOC3_VERBOSE(m_psTransform->nBackMapWidth,
116 m_psTransform->nBackMapHeight, sizeof(float)));
117
118 if (m_pafBackMapX == nullptr || m_pafBackMapY == nullptr ||
119 m_wgtsBackMap == nullptr)
120 {
121 return false;
122 }
123
124 const size_t nBMXYCount =
125 static_cast<size_t>(m_psTransform->nBackMapWidth) *
126 m_psTransform->nBackMapHeight;
127 for (size_t i = 0; i < nBMXYCount; i++)
128 {
129 m_pafBackMapX[i] = 0;
130 m_pafBackMapY[i] = 0;
131 m_wgtsBackMap[i] = 0.0;
132 }
133
134 backMapXAccessor.m_array = m_pafBackMapX;
135 backMapXAccessor.m_nXSize = m_psTransform->nBackMapWidth;
136
137 backMapYAccessor.m_array = m_pafBackMapY;
138 backMapYAccessor.m_nXSize = m_psTransform->nBackMapWidth;
139
140 backMapWeightAccessor.m_array = m_wgtsBackMap;
141 backMapWeightAccessor.m_nXSize = m_psTransform->nBackMapWidth;
142
143 return true;
144}
145
146/************************************************************************/
147/* FreeWghtsBackMap() */
148/************************************************************************/
149
150void GDALGeoLocCArrayAccessors::FreeWghtsBackMap()
151{
152 VSIFree(m_wgtsBackMap);
153 m_wgtsBackMap = nullptr;
154 backMapWeightAccessor.m_array = nullptr;
155 backMapWeightAccessor.m_nXSize = 0;
156}
157
158/************************************************************************/
159/* GetBackmapDataset() */
160/************************************************************************/
161
162GDALDataset *GDALGeoLocCArrayAccessors::GetBackmapDataset()
163{
164 auto poMEMDS = MEMDataset::Create("", m_psTransform->nBackMapWidth,
165 m_psTransform->nBackMapHeight, 0,
166 GDT_Float32, nullptr);
167
168 for (int i = 1; i <= 2; i++)
169 {
170 void *ptr = (i == 1) ? m_pafBackMapX : m_pafBackMapY;
171 GDALRasterBandH hMEMBand = MEMCreateRasterBandEx(
172 poMEMDS, i, static_cast<GByte *>(ptr), GDT_Float32, 0, 0, false);
173 poMEMDS->AddMEMBand(hMEMBand);
174 poMEMDS->GetRasterBand(i)->SetNoDataValue(INVALID_BMXY);
175 }
176 return poMEMDS;
177}
178
179/************************************************************************/
180/* Load() */
181/************************************************************************/
182
183bool GDALGeoLocCArrayAccessors::Load(bool bIsRegularGrid, bool bUseQuadtree)
184{
185 return LoadGeoloc(bIsRegularGrid) &&
186 ((bUseQuadtree && GDALGeoLocBuildQuadTree(m_psTransform)) ||
187 (!bUseQuadtree &&
188 GDALGeoLoc<AccessorType>::GenerateBackMap(m_psTransform)));
189}
190
191/************************************************************************/
192/* LoadGeoloc() */
193/************************************************************************/
194
195bool GDALGeoLocCArrayAccessors::LoadGeoloc(bool bIsRegularGrid)
196
197{
198 const int nXSize = m_psTransform->nGeoLocXSize;
199 const int nYSize = m_psTransform->nGeoLocYSize;
200
201 m_padfGeoLocY = static_cast<double *>(
202 VSI_MALLOC3_VERBOSE(sizeof(double), nXSize, nYSize));
203 m_padfGeoLocX = static_cast<double *>(
204 VSI_MALLOC3_VERBOSE(sizeof(double), nXSize, nYSize));
205
206 if (m_padfGeoLocX == nullptr || m_padfGeoLocY == nullptr)
207 {
208 return false;
209 }
210
211 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 static_cast<double *>(VSI_MALLOC2_VERBOSE(nXSize, sizeof(double)));
219 double *padfTempY =
220 static_cast<double *>(VSI_MALLOC2_VERBOSE(nYSize, sizeof(double)));
221 if (padfTempX == nullptr || padfTempY == nullptr)
222 {
223 CPLFree(padfTempX);
224 CPLFree(padfTempY);
225 return false;
226 }
227
228 CPLErr eErr =
229 GDALRasterIO(m_psTransform->hBand_X, GF_Read, 0, 0, nXSize, 1,
230 padfTempX, nXSize, 1, GDT_Float64, 0, 0);
231
232 for (size_t j = 0; j < static_cast<size_t>(nYSize); j++)
233 {
234 memcpy(m_padfGeoLocX + j * nXSize, padfTempX,
235 nXSize * sizeof(double));
236 }
237
238 if (eErr == CE_None)
239 {
240 eErr = GDALRasterIO(m_psTransform->hBand_Y, GF_Read, 0, 0, nYSize,
241 1, padfTempY, nYSize, 1, GDT_Float64, 0, 0);
242
243 for (size_t j = 0; j < static_cast<size_t>(nYSize); j++)
244 {
245 for (size_t i = 0; i < static_cast<size_t>(nXSize); i++)
246 {
247 m_padfGeoLocY[j * nXSize + i] = padfTempY[j];
248 }
249 }
250 }
251
252 CPLFree(padfTempX);
253 CPLFree(padfTempY);
254
255 if (eErr != CE_None)
256 return false;
257 }
258 else
259 {
260 if (GDALRasterIO(m_psTransform->hBand_X, GF_Read, 0, 0, nXSize, nYSize,
261 m_padfGeoLocX, nXSize, nYSize, GDT_Float64, 0,
262 0) != CE_None ||
263 GDALRasterIO(m_psTransform->hBand_Y, GF_Read, 0, 0, nXSize, nYSize,
264 m_padfGeoLocY, nXSize, nYSize, GDT_Float64, 0,
265 0) != CE_None)
266 return false;
267 }
268
269 geolocXAccessor.m_array = m_padfGeoLocX;
270 geolocXAccessor.m_nXSize = m_psTransform->nGeoLocXSize;
271
272 geolocYAccessor.m_array = m_padfGeoLocY;
273 geolocYAccessor.m_nXSize = m_psTransform->nGeoLocXSize;
274
275 GDALGeoLoc<GDALGeoLocCArrayAccessors>::LoadGeolocFinish(m_psTransform);
276 return true;
277}
278
A set of associated raster bands, usually from one file.
Definition: gdal_priv.h:495
#define CPLFree
Alias of VSIFree()
Definition: cpl_conv.h:82
CPLErr
Error category.
Definition: cpl_error.h:37
unsigned char GByte
Unsigned byte type.
Definition: cpl_port.h:169
#define VSI_MALLOC3_VERBOSE(nSize1, nSize2, nSize3)
VSI_MALLOC3_VERBOSE.
Definition: cpl_vsi.h:347
#define VSI_MALLOC2_VERBOSE(nSize1, nSize2)
VSI_MALLOC2_VERBOSE.
Definition: cpl_vsi.h:339
void VSIFree(void *)
Analog of free() for data allocated with VSIMalloc(), VSICalloc(), VSIRealloc()
Definition: cpl_vsisimple.cpp:828
@ GDT_Float64
Definition: gdal.h:59
@ GDT_Float32
Definition: gdal.h:58
@ GF_Read
Definition: gdal.h:118
void * GDALRasterBandH
Opaque type used for the C bindings of the C++ GDALRasterBand class.
Definition: gdal.h:379
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:457