Line data Source code
1 : /******************************************************************************
2 : * Name: gdalmultidim_array_gltorthorectification.cpp
3 : * Project: GDAL Core
4 : * Purpose: GLT orthorectification implementation
5 : * Author: Even Rouault <even.rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2023, Even Rouault <even.rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "gdal_priv.h"
14 : #include "gdal_pam_multidim.h"
15 :
16 : #include <algorithm>
17 : #include <limits>
18 :
19 : /************************************************************************/
20 : /* GLTOrthoRectifiedArray */
21 : /************************************************************************/
22 :
23 : class GLTOrthoRectifiedArray final : public GDALPamMDArray
24 : {
25 : private:
26 : std::shared_ptr<GDALMDArray> m_poParent{};
27 : std::vector<std::shared_ptr<GDALDimension>> m_apoDims;
28 : std::vector<GUInt64> m_anBlockSize;
29 : GDALExtendedDataType m_dt;
30 : std::shared_ptr<OGRSpatialReference> m_poSRS{};
31 : std::shared_ptr<GDALMDArray> m_poVarX{};
32 : std::shared_ptr<GDALMDArray> m_poVarY{};
33 : std::shared_ptr<GDALMDArray> m_poGLTX{};
34 : std::shared_ptr<GDALMDArray> m_poGLTY{};
35 : int m_nGLTIndexOffset = 0;
36 : std::vector<GByte> m_abyBandValidity{};
37 :
38 : protected:
39 9 : GLTOrthoRectifiedArray(
40 : const std::shared_ptr<GDALMDArray> &poParent,
41 : const std::vector<std::shared_ptr<GDALDimension>> &apoDims,
42 : const std::vector<GUInt64> &anBlockSize)
43 18 : : GDALAbstractMDArray(std::string(), "GLTOrthoRectifiedArray view of " +
44 9 : poParent->GetFullName()),
45 18 : GDALPamMDArray(std::string(),
46 18 : "GLTOrthoRectifiedArray view of " +
47 9 : poParent->GetFullName(),
48 18 : GDALPamMultiDim::GetPAM(poParent)),
49 9 : m_poParent(std::move(poParent)), m_apoDims(apoDims),
50 45 : m_anBlockSize(anBlockSize), m_dt(m_poParent->GetDataType())
51 : {
52 9 : CPLAssert(apoDims.size() == m_poParent->GetDimensionCount());
53 9 : CPLAssert(anBlockSize.size() == m_poParent->GetDimensionCount());
54 9 : }
55 :
56 : bool IRead(const GUInt64 *arrayStartIdx, const size_t *count,
57 : const GInt64 *arrayStep, const GPtrDiff_t *bufferStride,
58 : const GDALExtendedDataType &bufferDataType,
59 : void *pDstBuffer) const override;
60 :
61 : public:
62 : static std::shared_ptr<GDALMDArray>
63 9 : Create(const std::shared_ptr<GDALMDArray> &poParent,
64 : const std::shared_ptr<GDALGroup> &poRootGroup,
65 : const std::shared_ptr<GDALMDArray> &poGLTX,
66 : const std::shared_ptr<GDALMDArray> &poGLTY, int nGLTIndexOffset,
67 : const std::vector<double> &adfGeoTransform,
68 : CSLConstList papszOptions)
69 : {
70 18 : std::vector<std::shared_ptr<GDALDimension>> apoNewDims;
71 :
72 : auto poDimY = std::make_shared<GDALDimensionWeakIndexingVar>(
73 0 : std::string(), "lat", GDAL_DIM_TYPE_HORIZONTAL_Y, "NORTH",
74 18 : poGLTX->GetDimensions()[0]->GetSize());
75 : auto varY = GDALMDArrayRegularlySpaced::Create(
76 27 : std::string(), poDimY->GetName(), poDimY,
77 36 : adfGeoTransform[3] + adfGeoTransform[5] / 2, adfGeoTransform[5], 0);
78 9 : poDimY->SetIndexingVariable(varY);
79 9 : apoNewDims.emplace_back(poDimY);
80 :
81 : auto poDimX = std::make_shared<GDALDimensionWeakIndexingVar>(
82 0 : std::string(), "lon", GDAL_DIM_TYPE_HORIZONTAL_X, "EAST",
83 18 : poGLTX->GetDimensions()[1]->GetSize());
84 : auto varX = GDALMDArrayRegularlySpaced::Create(
85 27 : std::string(), poDimX->GetName(), poDimX,
86 36 : adfGeoTransform[0] + adfGeoTransform[1] / 2, adfGeoTransform[1], 0);
87 9 : poDimX->SetIndexingVariable(varX);
88 9 : apoNewDims.emplace_back(poDimX);
89 :
90 9 : if (poParent->GetDimensionCount() == 3)
91 8 : apoNewDims.emplace_back(poParent->GetDimensions()[2]);
92 :
93 : std::vector<GUInt64> anBlockSize = std::vector<GUInt64>{
94 9 : std::min<GUInt64>(apoNewDims[0]->GetSize(), 512),
95 18 : std::min<GUInt64>(apoNewDims[1]->GetSize(), 512)};
96 9 : if (poParent->GetDimensionCount() == 3)
97 : {
98 8 : anBlockSize.push_back(poParent->GetDimensions()[2]->GetSize());
99 : }
100 :
101 : auto newAr(std::shared_ptr<GLTOrthoRectifiedArray>(
102 18 : new GLTOrthoRectifiedArray(poParent, apoNewDims, anBlockSize)));
103 9 : newAr->SetSelf(newAr);
104 9 : newAr->m_poVarX = varX;
105 9 : newAr->m_poVarY = varY;
106 9 : newAr->m_poGLTX = poGLTX;
107 9 : newAr->m_poGLTY = poGLTY;
108 9 : newAr->m_nGLTIndexOffset = nGLTIndexOffset;
109 18 : OGRSpatialReference oSRS;
110 9 : oSRS.importFromEPSG(4326);
111 9 : newAr->m_poSRS.reset(oSRS.Clone());
112 9 : newAr->m_poSRS->SetDataAxisToSRSAxisMapping(
113 : {/*latIdx = */ 1, /* lonIdx = */ 2});
114 9 : if (CPLTestBool(CSLFetchNameValueDef(papszOptions,
115 17 : "USE_GOOD_WAVELENGTHS", "YES")) &&
116 8 : newAr->m_poParent->GetDimensionCount() == 3)
117 : {
118 : const auto poGoodWaveLengths = poRootGroup->OpenMDArrayFromFullname(
119 21 : "/sensor_band_parameters/good_wavelengths");
120 10 : if (poGoodWaveLengths &&
121 13 : poGoodWaveLengths->GetDimensionCount() == 1 &&
122 3 : poGoodWaveLengths->GetDimensions()[0]->GetSize() ==
123 6 : newAr->m_poParent->GetDimensions()[2]->GetSize() &&
124 3 : poGoodWaveLengths->GetDimensions()[0]->GetSize() <
125 10 : 1000 * 1000 &&
126 3 : poGoodWaveLengths->GetDataType().GetClass() == GEDTC_NUMERIC)
127 : {
128 3 : const GUInt64 arrayStartIdx[] = {0};
129 : const size_t count[] = {static_cast<size_t>(
130 3 : poGoodWaveLengths->GetDimensions()[0]->GetSize())};
131 3 : const GInt64 arrayStep[] = {1};
132 3 : const GPtrDiff_t bufferStride[] = {1};
133 3 : newAr->m_abyBandValidity.resize(count[0], 1);
134 6 : CPL_IGNORE_RET_VAL(poGoodWaveLengths->Read(
135 : arrayStartIdx, count, arrayStep, bufferStride,
136 6 : GDALExtendedDataType::Create(GDT_UInt8),
137 3 : newAr->m_abyBandValidity.data()));
138 : }
139 : }
140 18 : return newAr;
141 : }
142 :
143 2 : bool IsWritable() const override
144 : {
145 2 : return false;
146 : }
147 :
148 16 : const std::string &GetFilename() const override
149 : {
150 16 : return m_poParent->GetFilename();
151 : }
152 :
153 : const std::vector<std::shared_ptr<GDALDimension>> &
154 107 : GetDimensions() const override
155 : {
156 107 : return m_apoDims;
157 : }
158 :
159 77 : const GDALExtendedDataType &GetDataType() const override
160 : {
161 77 : return m_dt;
162 : }
163 :
164 4 : std::shared_ptr<OGRSpatialReference> GetSpatialRef() const override
165 : {
166 4 : return m_poSRS;
167 : }
168 :
169 5 : std::vector<GUInt64> GetBlockSize() const override
170 : {
171 5 : return m_anBlockSize;
172 : }
173 :
174 : std::shared_ptr<GDALAttribute>
175 3 : GetAttribute(const std::string &osName) const override
176 : {
177 3 : return m_poParent->GetAttribute(osName);
178 : }
179 :
180 : std::vector<std::shared_ptr<GDALAttribute>>
181 2 : GetAttributes(CSLConstList papszOptions = nullptr) const override
182 : {
183 2 : return m_poParent->GetAttributes(papszOptions);
184 : }
185 :
186 3 : const std::string &GetUnit() const override
187 : {
188 3 : return m_poParent->GetUnit();
189 : }
190 :
191 19 : const void *GetRawNoDataValue() const override
192 : {
193 19 : return m_poParent->GetRawNoDataValue();
194 : }
195 :
196 0 : double GetOffset(bool *pbHasOffset,
197 : GDALDataType *peStorageType) const override
198 : {
199 0 : return m_poParent->GetOffset(pbHasOffset, peStorageType);
200 : }
201 :
202 0 : double GetScale(bool *pbHasScale,
203 : GDALDataType *peStorageType) const override
204 : {
205 0 : return m_poParent->GetScale(pbHasScale, peStorageType);
206 : }
207 : };
208 :
209 : /************************************************************************/
210 : /* GDALMDArrayResampled::IRead() */
211 : /************************************************************************/
212 :
213 16 : bool GLTOrthoRectifiedArray::IRead(const GUInt64 *arrayStartIdx,
214 : const size_t *count, const GInt64 *arrayStep,
215 : const GPtrDiff_t *bufferStride,
216 : const GDALExtendedDataType &bufferDataType,
217 : void *pDstBuffer) const
218 : {
219 16 : if (bufferDataType.GetClass() != GEDTC_NUMERIC)
220 0 : return false;
221 :
222 16 : const auto eBufferDT = bufferDataType.GetNumericDataType();
223 16 : auto pRawNoDataValue = GetRawNoDataValue();
224 32 : std::vector<GByte> abyNoData(16);
225 16 : if (pRawNoDataValue)
226 16 : GDALCopyWords64(pRawNoDataValue, GetDataType().GetNumericDataType(), 0,
227 16 : abyNoData.data(), eBufferDT, 0, 1);
228 :
229 16 : const auto nBufferDTSize = bufferDataType.GetSize();
230 : const int nCopyWordsDstStride =
231 16 : m_apoDims.size() == 3
232 16 : ? static_cast<int>(bufferStride[2] * nBufferDTSize)
233 16 : : 0;
234 : const int nOutBandCount =
235 16 : m_apoDims.size() == 3 ? static_cast<int>(count[2]) : 1;
236 :
237 : const auto FillBufferWithNodata =
238 5 : [count, bufferStride, pDstBuffer, nBufferDTSize, &eBufferDT,
239 48 : nCopyWordsDstStride, nOutBandCount, &abyNoData]()
240 : {
241 12 : for (size_t iY = 0; iY < count[0]; ++iY)
242 : {
243 13 : if (nOutBandCount == 1 && bufferStride[1] > 0 &&
244 6 : static_cast<int>(nBufferDTSize) <
245 6 : std::numeric_limits<int>::max() / bufferStride[1])
246 : {
247 6 : GDALCopyWords64(
248 6 : abyNoData.data(), eBufferDT, 0,
249 : static_cast<GByte *>(pDstBuffer) +
250 6 : iY * bufferStride[0] * static_cast<int>(nBufferDTSize),
251 : eBufferDT,
252 6 : static_cast<int>(bufferStride[1] * nBufferDTSize),
253 6 : count[1]);
254 : }
255 : else
256 : {
257 3 : for (size_t iX = 0; iX < count[1]; ++iX)
258 : {
259 2 : GByte *pabyDstBuffer =
260 : static_cast<GByte *>(pDstBuffer) +
261 2 : (iY * bufferStride[0] + iX * bufferStride[1]) *
262 2 : static_cast<int>(nBufferDTSize);
263 2 : GDALCopyWords64(abyNoData.data(), eBufferDT, 0,
264 : pabyDstBuffer, eBufferDT,
265 : nCopyWordsDstStride, nOutBandCount);
266 : }
267 : }
268 : }
269 21 : };
270 :
271 16 : bool bSomeBandInvalids = false;
272 16 : if (!m_abyBandValidity.empty())
273 : {
274 3 : size_t nCountInvalid = 0;
275 7 : for (size_t i = 0; i < count[2]; ++i)
276 : {
277 4 : const size_t iBand =
278 4 : static_cast<size_t>(arrayStartIdx[2] + i * arrayStep[2]);
279 4 : CPLAssert(iBand < m_abyBandValidity.size());
280 4 : if (!m_abyBandValidity[iBand])
281 2 : nCountInvalid++;
282 : }
283 3 : if (nCountInvalid == count[2])
284 : {
285 1 : FillBufferWithNodata();
286 1 : return true;
287 : }
288 2 : bSomeBandInvalids = true;
289 : }
290 :
291 15 : const size_t nXYValsCount = count[0] * count[1];
292 30 : const auto eInt32DT = GDALExtendedDataType::Create(GDT_Int32);
293 30 : std::vector<int32_t> anGLTX;
294 30 : std::vector<int32_t> anGLTY;
295 : try
296 : {
297 15 : anGLTX.resize(nXYValsCount);
298 15 : anGLTY.resize(nXYValsCount);
299 : }
300 0 : catch (const std::bad_alloc &e)
301 : {
302 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
303 0 : "GLTOrthoRectifiedArray::IRead(): %s", e.what());
304 0 : return false;
305 : }
306 30 : if (!m_poGLTX->Read(arrayStartIdx, count, arrayStep, nullptr, eInt32DT,
307 45 : anGLTX.data()) ||
308 30 : !m_poGLTY->Read(arrayStartIdx, count, arrayStep, nullptr, eInt32DT,
309 15 : anGLTY.data()))
310 : {
311 0 : return false;
312 : }
313 :
314 15 : int32_t nMinX = std::numeric_limits<int32_t>::max();
315 15 : int32_t nMaxX = std::numeric_limits<int32_t>::min();
316 15 : const auto nXSize = m_poParent->GetDimensions()[0]->GetSize();
317 104 : for (size_t i = 0; i < nXYValsCount; ++i)
318 : {
319 : const int64_t nX64 =
320 89 : static_cast<int64_t>(anGLTX[i]) + m_nGLTIndexOffset;
321 89 : if (nX64 >= 0 && static_cast<uint64_t>(nX64) < nXSize)
322 : {
323 38 : const int32_t nX = static_cast<int32_t>(nX64);
324 38 : if (nX < nMinX)
325 11 : nMinX = nX;
326 38 : if (nX > nMaxX)
327 20 : nMaxX = nX;
328 : }
329 : }
330 :
331 15 : int32_t nMinY = std::numeric_limits<int32_t>::max();
332 15 : int32_t nMaxY = std::numeric_limits<int32_t>::min();
333 15 : const auto nYSize = m_poParent->GetDimensions()[0]->GetSize();
334 104 : for (size_t i = 0; i < nXYValsCount; ++i)
335 : {
336 : const int64_t nY64 =
337 89 : static_cast<int64_t>(anGLTY[i]) + m_nGLTIndexOffset;
338 89 : if (nY64 >= 0 && static_cast<uint64_t>(nY64) < nYSize)
339 : {
340 38 : const int32_t nY = static_cast<int32_t>(nY64);
341 38 : if (nY < nMinY)
342 20 : nMinY = nY;
343 38 : if (nY > nMaxY)
344 11 : nMaxY = nY;
345 : }
346 : }
347 :
348 15 : if (nMinX > nMaxX || nMinY > nMaxY)
349 : {
350 4 : FillBufferWithNodata();
351 4 : return true;
352 : }
353 :
354 : GUInt64 parentArrayIdxStart[3] = {
355 11 : static_cast<GUInt64>(nMinY), static_cast<GUInt64>(nMinX),
356 11 : m_apoDims.size() == 3 ? arrayStartIdx[2] : 0};
357 11 : size_t parentCount[3] = {static_cast<size_t>(nMaxY - nMinY + 1),
358 11 : static_cast<size_t>(nMaxX - nMinX + 1),
359 11 : m_apoDims.size() == 3 ? count[2] : 1};
360 11 : GInt64 parentArrayStep[3] = {1, 1,
361 11 : m_apoDims.size() == 3 ? arrayStep[2] : 0};
362 :
363 11 : size_t nParentValueSize = nBufferDTSize;
364 44 : for (int i = 0; i < 3; ++i)
365 : {
366 66 : if (parentCount[i] >
367 33 : std::numeric_limits<size_t>::max() / nParentValueSize)
368 : {
369 0 : CPLError(
370 : CE_Failure, CPLE_OutOfMemory,
371 : "GLTOrthoRectifiedArray::IRead(): too big temporary array");
372 0 : return false;
373 : }
374 33 : nParentValueSize *= parentCount[i];
375 : }
376 :
377 11 : GPtrDiff_t parentStride[3] = {
378 11 : static_cast<GPtrDiff_t>(parentCount[1] * parentCount[2]),
379 11 : static_cast<GPtrDiff_t>(parentCount[2]), 1};
380 22 : std::vector<GByte> parentValues;
381 : try
382 : {
383 11 : parentValues.resize(nParentValueSize);
384 : }
385 0 : catch (const std::bad_alloc &e)
386 : {
387 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
388 0 : "GLTOrthoRectifiedArray::IRead(): %s", e.what());
389 0 : return false;
390 : }
391 22 : if (!m_poParent->Read(parentArrayIdxStart, parentCount, parentArrayStep,
392 11 : parentStride, bufferDataType, parentValues.data()))
393 : {
394 0 : return false;
395 : }
396 :
397 11 : size_t iGLTIndex = 0;
398 11 : const size_t nXCount = parentCount[1];
399 11 : const size_t nBandCount = m_apoDims.size() == 3 ? parentCount[2] : 1;
400 40 : for (size_t iY = 0; iY < count[0]; ++iY)
401 : {
402 112 : for (size_t iX = 0; iX < count[1]; ++iX, ++iGLTIndex)
403 : {
404 : const int64_t nX64 =
405 83 : static_cast<int64_t>(anGLTX[iGLTIndex]) + m_nGLTIndexOffset;
406 : const int64_t nY64 =
407 83 : static_cast<int64_t>(anGLTY[iGLTIndex]) + m_nGLTIndexOffset;
408 83 : GByte *pabyDstBuffer =
409 : static_cast<GByte *>(pDstBuffer) +
410 83 : (iY * bufferStride[0] + iX * bufferStride[1]) *
411 83 : static_cast<int>(nBufferDTSize);
412 83 : if (nX64 >= nMinX && nX64 <= nMaxX && nY64 >= nMinY &&
413 38 : nY64 <= nMaxY)
414 : {
415 38 : const int32_t iSrcX = static_cast<int32_t>(nX64) - nMinX;
416 38 : const int32_t iSrcY = static_cast<int32_t>(nY64) - nMinY;
417 : const GByte *pabySrcBuffer =
418 38 : parentValues.data() +
419 38 : (iSrcY * nXCount + iSrcX) * nBandCount * nBufferDTSize;
420 38 : GDALCopyWords64(pabySrcBuffer, eBufferDT,
421 : static_cast<int>(nBufferDTSize), pabyDstBuffer,
422 : eBufferDT, nCopyWordsDstStride, nOutBandCount);
423 38 : if (bSomeBandInvalids)
424 : {
425 20 : for (size_t i = 0; i < count[2]; ++i)
426 : {
427 12 : const size_t iBand = static_cast<size_t>(
428 12 : arrayStartIdx[2] + i * arrayStep[2]);
429 12 : if (!m_abyBandValidity[iBand])
430 : {
431 4 : GDALCopyWords64(abyNoData.data(), eBufferDT, 0,
432 4 : pabyDstBuffer +
433 4 : i * nCopyWordsDstStride,
434 : eBufferDT, 0, 1);
435 : }
436 : }
437 38 : }
438 : }
439 : else
440 : {
441 45 : GDALCopyWords64(abyNoData.data(), eBufferDT, 0, pabyDstBuffer,
442 : eBufferDT, nCopyWordsDstStride, nOutBandCount);
443 : }
444 : }
445 : }
446 :
447 11 : return true;
448 : }
449 :
450 : /************************************************************************/
451 : /* CreateGLTOrthorectified() */
452 : /************************************************************************/
453 :
454 : //! @cond Doxygen_Suppress
455 :
456 9 : /* static */ std::shared_ptr<GDALMDArray> GDALMDArray::CreateGLTOrthorectified(
457 : const std::shared_ptr<GDALMDArray> &poParent,
458 : const std::shared_ptr<GDALGroup> &poRootGroup,
459 : const std::shared_ptr<GDALMDArray> &poGLTX,
460 : const std::shared_ptr<GDALMDArray> &poGLTY, int nGLTIndexOffset,
461 : const std::vector<double> &adfGeoTransform, CSLConstList papszOptions)
462 : {
463 : return GLTOrthoRectifiedArray::Create(poParent, poRootGroup, poGLTX, poGLTY,
464 : nGLTIndexOffset, adfGeoTransform,
465 9 : papszOptions);
466 : }
467 :
468 : //! @endcond
|
