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 : * 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 :
29 : /*! @cond Doxygen_Suppress */
30 :
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 230 : CArrayAccessor(Type *array, size_t nXSize)
55 230 : : m_array(array), m_nXSize(nXSize)
56 : {
57 230 : }
58 :
59 23309440 : inline Type Get(int nX, int nY, bool *pbSuccess = nullptr)
60 : {
61 23309440 : if (pbSuccess)
62 0 : *pbSuccess = true;
63 23309440 : return m_array[nY * m_nXSize + nX];
64 : }
65 :
66 1930990 : inline bool Set(int nX, int nY, Type val)
67 : {
68 1930990 : m_array[nY * m_nXSize + nX] = val;
69 1930990 : 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 46 : explicit GDALGeoLocCArrayAccessors(GDALGeoLocTransformInfo *psTransform)
80 46 : : m_psTransform(psTransform), geolocXAccessor(nullptr, 0),
81 : geolocYAccessor(nullptr, 0), backMapXAccessor(nullptr, 0),
82 46 : backMapYAccessor(nullptr, 0), backMapWeightAccessor(nullptr, 0)
83 : {
84 46 : }
85 :
86 46 : ~GDALGeoLocCArrayAccessors()
87 46 : {
88 46 : VSIFree(m_pafBackMapX);
89 46 : VSIFree(m_pafBackMapY);
90 46 : VSIFree(m_padfGeoLocX);
91 46 : VSIFree(m_padfGeoLocY);
92 46 : VSIFree(m_wgtsBackMap);
93 46 : }
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 :
105 42 : static void FlushBackmapCaches()
106 : {
107 42 : }
108 :
109 42 : static void ReleaseBackmapDataset(GDALDataset *poDS)
110 : {
111 42 : delete poDS;
112 42 : }
113 :
114 : void FreeWghtsBackMap();
115 : };
116 :
117 : /************************************************************************/
118 : /* AllocateBackMap() */
119 : /************************************************************************/
120 :
121 42 : bool GDALGeoLocCArrayAccessors::AllocateBackMap()
122 : {
123 42 : m_pafBackMapX = static_cast<float *>(
124 42 : VSI_MALLOC3_VERBOSE(m_psTransform->nBackMapWidth,
125 : m_psTransform->nBackMapHeight, sizeof(float)));
126 42 : m_pafBackMapY = static_cast<float *>(
127 42 : VSI_MALLOC3_VERBOSE(m_psTransform->nBackMapWidth,
128 : m_psTransform->nBackMapHeight, sizeof(float)));
129 :
130 42 : m_wgtsBackMap = static_cast<float *>(
131 42 : VSI_MALLOC3_VERBOSE(m_psTransform->nBackMapWidth,
132 : m_psTransform->nBackMapHeight, sizeof(float)));
133 :
134 42 : if (m_pafBackMapX == nullptr || m_pafBackMapY == nullptr ||
135 42 : m_wgtsBackMap == nullptr)
136 : {
137 0 : return false;
138 : }
139 :
140 42 : const size_t nBMXYCount =
141 42 : static_cast<size_t>(m_psTransform->nBackMapWidth) *
142 42 : m_psTransform->nBackMapHeight;
143 186014 : for (size_t i = 0; i < nBMXYCount; i++)
144 : {
145 185972 : m_pafBackMapX[i] = 0;
146 185972 : m_pafBackMapY[i] = 0;
147 185972 : m_wgtsBackMap[i] = 0.0;
148 : }
149 :
150 42 : backMapXAccessor.m_array = m_pafBackMapX;
151 42 : backMapXAccessor.m_nXSize = m_psTransform->nBackMapWidth;
152 :
153 42 : backMapYAccessor.m_array = m_pafBackMapY;
154 42 : backMapYAccessor.m_nXSize = m_psTransform->nBackMapWidth;
155 :
156 42 : backMapWeightAccessor.m_array = m_wgtsBackMap;
157 42 : backMapWeightAccessor.m_nXSize = m_psTransform->nBackMapWidth;
158 :
159 42 : return true;
160 : }
161 :
162 : /************************************************************************/
163 : /* FreeWghtsBackMap() */
164 : /************************************************************************/
165 :
166 42 : void GDALGeoLocCArrayAccessors::FreeWghtsBackMap()
167 : {
168 42 : VSIFree(m_wgtsBackMap);
169 42 : m_wgtsBackMap = nullptr;
170 42 : backMapWeightAccessor.m_array = nullptr;
171 42 : backMapWeightAccessor.m_nXSize = 0;
172 42 : }
173 :
174 : /************************************************************************/
175 : /* GetBackmapDataset() */
176 : /************************************************************************/
177 :
178 42 : GDALDataset *GDALGeoLocCArrayAccessors::GetBackmapDataset()
179 : {
180 84 : auto poMEMDS = MEMDataset::Create("", m_psTransform->nBackMapWidth,
181 42 : m_psTransform->nBackMapHeight, 0,
182 : GDT_Float32, nullptr);
183 :
184 126 : for (int i = 1; i <= 2; i++)
185 : {
186 84 : void *ptr = (i == 1) ? m_pafBackMapX : m_pafBackMapY;
187 84 : GDALRasterBandH hMEMBand = MEMCreateRasterBandEx(
188 : poMEMDS, i, static_cast<GByte *>(ptr), GDT_Float32, 0, 0, false);
189 84 : poMEMDS->AddMEMBand(hMEMBand);
190 84 : poMEMDS->GetRasterBand(i)->SetNoDataValue(INVALID_BMXY);
191 : }
192 42 : return poMEMDS;
193 : }
194 :
195 : /************************************************************************/
196 : /* Load() */
197 : /************************************************************************/
198 :
199 46 : bool GDALGeoLocCArrayAccessors::Load(bool bIsRegularGrid, bool bUseQuadtree)
200 : {
201 92 : return LoadGeoloc(bIsRegularGrid) &&
202 4 : ((bUseQuadtree && GDALGeoLocBuildQuadTree(m_psTransform)) ||
203 42 : (!bUseQuadtree &&
204 88 : GDALGeoLoc<AccessorType>::GenerateBackMap(m_psTransform)));
205 : }
206 :
207 : /************************************************************************/
208 : /* LoadGeoloc() */
209 : /************************************************************************/
210 :
211 46 : bool GDALGeoLocCArrayAccessors::LoadGeoloc(bool bIsRegularGrid)
212 :
213 : {
214 46 : const int nXSize = m_psTransform->nGeoLocXSize;
215 46 : const int nYSize = m_psTransform->nGeoLocYSize;
216 :
217 46 : m_padfGeoLocY = static_cast<double *>(
218 46 : VSI_MALLOC3_VERBOSE(sizeof(double), nXSize, nYSize));
219 46 : m_padfGeoLocX = static_cast<double *>(
220 46 : VSI_MALLOC3_VERBOSE(sizeof(double), nXSize, nYSize));
221 :
222 46 : if (m_padfGeoLocX == nullptr || m_padfGeoLocY == nullptr)
223 : {
224 0 : return false;
225 : }
226 :
227 46 : if (bIsRegularGrid)
228 : {
229 : // Case of regular grid.
230 : // The XBAND contains the x coordinates for all lines.
231 : // The YBAND contains the y coordinates for all columns.
232 :
233 : double *padfTempX =
234 14 : static_cast<double *>(VSI_MALLOC2_VERBOSE(nXSize, sizeof(double)));
235 : double *padfTempY =
236 14 : static_cast<double *>(VSI_MALLOC2_VERBOSE(nYSize, sizeof(double)));
237 14 : if (padfTempX == nullptr || padfTempY == nullptr)
238 : {
239 0 : CPLFree(padfTempX);
240 0 : CPLFree(padfTempY);
241 0 : return false;
242 : }
243 :
244 : CPLErr eErr =
245 14 : GDALRasterIO(m_psTransform->hBand_X, GF_Read, 0, 0, nXSize, 1,
246 : padfTempX, nXSize, 1, GDT_Float64, 0, 0);
247 :
248 262 : for (size_t j = 0; j < static_cast<size_t>(nYSize); j++)
249 : {
250 248 : memcpy(m_padfGeoLocX + j * nXSize, padfTempX,
251 248 : nXSize * sizeof(double));
252 : }
253 :
254 14 : if (eErr == CE_None)
255 : {
256 14 : eErr = GDALRasterIO(m_psTransform->hBand_Y, GF_Read, 0, 0, nYSize,
257 : 1, padfTempY, nYSize, 1, GDT_Float64, 0, 0);
258 :
259 262 : for (size_t j = 0; j < static_cast<size_t>(nYSize); j++)
260 : {
261 31480 : for (size_t i = 0; i < static_cast<size_t>(nXSize); i++)
262 : {
263 31232 : m_padfGeoLocY[j * nXSize + i] = padfTempY[j];
264 : }
265 : }
266 : }
267 :
268 14 : CPLFree(padfTempX);
269 14 : CPLFree(padfTempY);
270 :
271 14 : if (eErr != CE_None)
272 0 : return false;
273 : }
274 : else
275 : {
276 96 : if (GDALRasterIO(m_psTransform->hBand_X, GF_Read, 0, 0, nXSize, nYSize,
277 32 : m_padfGeoLocX, nXSize, nYSize, GDT_Float64, 0,
278 64 : 0) != CE_None ||
279 32 : GDALRasterIO(m_psTransform->hBand_Y, GF_Read, 0, 0, nXSize, nYSize,
280 32 : m_padfGeoLocY, nXSize, nYSize, GDT_Float64, 0,
281 : 0) != CE_None)
282 0 : return false;
283 : }
284 :
285 46 : geolocXAccessor.m_array = m_padfGeoLocX;
286 46 : geolocXAccessor.m_nXSize = m_psTransform->nGeoLocXSize;
287 :
288 46 : geolocYAccessor.m_array = m_padfGeoLocY;
289 46 : geolocYAccessor.m_nXSize = m_psTransform->nGeoLocXSize;
290 :
291 46 : GDALGeoLoc<GDALGeoLocCArrayAccessors>::LoadGeolocFinish(m_psTransform);
292 46 : return true;
293 : }
294 :
295 : /*! @endcond */
|
