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