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 void* ptr = (i == 1) ? m_pafBackMapX : m_pafBackMapY;
181 GDALRasterBandH hMEMBand = MEMCreateRasterBandEx( poMEMDS,
182 i, static_cast<GByte*>(ptr),
183 GDT_Float32,
184 0, 0,
185 false );
186 poMEMDS->AddMEMBand(hMEMBand);
187 poMEMDS->GetRasterBand(i)->SetNoDataValue(INVALID_BMXY);
188 }
189 return poMEMDS;
190 }
191
192 /************************************************************************/
193 /* Load() */
194 /************************************************************************/
195
196 bool GDALGeoLocCArrayAccessors::Load(bool bIsRegularGrid, bool bUseQuadtree)
197 {
198 return LoadGeoloc(bIsRegularGrid) &&
199 ((bUseQuadtree && GDALGeoLocBuildQuadTree(m_psTransform)) ||
200 (!bUseQuadtree && GDALGeoLoc<AccessorType>::GenerateBackMap(m_psTransform)));
201 }
202
203 /************************************************************************/
204 /* LoadGeoloc() */
205 /************************************************************************/
206
207 bool GDALGeoLocCArrayAccessors::LoadGeoloc(bool bIsRegularGrid)
208
209 {
210 const int nXSize = m_psTransform->nGeoLocXSize;
211 const int nYSize = m_psTransform->nGeoLocYSize;
212
213 m_padfGeoLocY = static_cast<double *>(
214 VSI_MALLOC3_VERBOSE(sizeof(double), nXSize, nYSize));
215 m_padfGeoLocX = static_cast<double *>(
216 VSI_MALLOC3_VERBOSE(sizeof(double), nXSize, nYSize));
217
218 if( m_padfGeoLocX == nullptr ||
219 m_padfGeoLocY == nullptr )
220 {
221 return false;
222 }
223
224 if( bIsRegularGrid )
225 {
226 // Case of regular grid.
227 // The XBAND contains the x coordinates for all lines.
228 // The YBAND contains the y coordinates for all columns.
229
230 double* padfTempX = static_cast<double *>(
231 VSI_MALLOC2_VERBOSE(nXSize, sizeof(double)));
232 double* padfTempY = static_cast<double *>(
233 VSI_MALLOC2_VERBOSE(nYSize, sizeof(double)));
234 if( padfTempX == nullptr || padfTempY == nullptr )
235 {
236 CPLFree(padfTempX);
237 CPLFree(padfTempY);
238 return false;
239 }
240
241 CPLErr eErr =
242 GDALRasterIO( m_psTransform->hBand_X, GF_Read,
243 0, 0, nXSize, 1,
244 padfTempX, nXSize, 1,
245 GDT_Float64, 0, 0 );
246
247 for( size_t j = 0; j < static_cast<size_t>(nYSize); j++ )
248 {
249 memcpy( m_padfGeoLocX + j * nXSize,
250 padfTempX,
251 nXSize * sizeof(double) );
252 }
253
254 if( eErr == CE_None )
255 {
256 eErr = GDALRasterIO( m_psTransform->hBand_Y, GF_Read,
257 0, 0, nYSize, 1,
258 padfTempY, nYSize, 1,
259 GDT_Float64, 0, 0 );
260
261 for( size_t j = 0; j < static_cast<size_t>(nYSize); j++ )
262 {
263 for( size_t i = 0; i < static_cast<size_t>(nXSize); i++ )
264 {
265 m_padfGeoLocY[j * nXSize + i] = padfTempY[j];
266 }
267 }
268 }
269
270 CPLFree(padfTempX);
271 CPLFree(padfTempY);
272
273 if( eErr != CE_None )
274 return false;
275 }
276 else
277 {
278 if( GDALRasterIO( m_psTransform->hBand_X, GF_Read,
279 0, 0, nXSize, nYSize,
280 m_padfGeoLocX, nXSize, nYSize,
281 GDT_Float64, 0, 0 ) != CE_None
282 || GDALRasterIO( m_psTransform->hBand_Y, GF_Read,
283 0, 0, nXSize, nYSize,
284 m_padfGeoLocY, nXSize, nYSize,
285 GDT_Float64, 0, 0 ) != CE_None )
286 return false;
287 }
288
289 geolocXAccessor.m_array = m_padfGeoLocX;
290 geolocXAccessor.m_nXSize = m_psTransform->nGeoLocXSize;
291
292 geolocYAccessor.m_array = m_padfGeoLocY;
293 geolocYAccessor.m_nXSize = m_psTransform->nGeoLocXSize;
294
295 return GDALGeoLoc<GDALGeoLocCArrayAccessors>::LoadGeolocFinish(m_psTransform);
296 }
297
GByte
unsigned char GByte
Unsigned byte type.
Definition: cpl_port.h:203
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:342
VSI_MALLOC2_VERBOSE
#define VSI_MALLOC2_VERBOSE(nSize1, nSize2)
VSI_MALLOC2_VERBOSE.
Definition: cpl_vsi.h:318
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:826
VSI_MALLOC3_VERBOSE
#define VSI_MALLOC3_VERBOSE(nSize1, nSize2, nSize3)
VSI_MALLOC3_VERBOSE.
Definition: cpl_vsi.h:323
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:400
GDALRasterBandH
void * GDALRasterBandH
Opaque type used for the C bindings of the C++ GDALRasterBand class.
Definition: gdal.h:272
CPLFree
#define CPLFree
Alias of VSIFree()
Definition: cpl_conv.h:83