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 return GDALGeoLoc<GDALGeoLocCArrayAccessors>::LoadGeolocFinish(
290 m_psTransform);
291 }
292
GByte
unsigned char GByte
Unsigned byte type.
Definition: cpl_port.h:205
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:347
VSI_MALLOC2_VERBOSE
#define VSI_MALLOC2_VERBOSE(nSize1, nSize2)
VSI_MALLOC2_VERBOSE.
Definition: cpl_vsi.h:335
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:832
VSI_MALLOC3_VERBOSE
#define VSI_MALLOC3_VERBOSE(nSize1, nSize2, nSize3)
VSI_MALLOC3_VERBOSE.
Definition: cpl_vsi.h:343
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:403
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:96