Line data Source code
1 : /******************************************************************************
2 : * $Id$
3 : *
4 : * Project: GDAL DEM Interpolation
5 : * Purpose: Interpolation algorithms with cache
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2001, Frank Warmerdam
10 : * Copyright (c) 2008-2012, Even Rouault <even dot rouault at spatialys.com>
11 : * Copyright (c) 2024, Javier Jimenez Shaw
12 : *
13 : * SPDX-License-Identifier: MIT
14 : ****************************************************************************/
15 :
16 : #include "gdal_interpolateatpoint.h"
17 :
18 : #include "gdalresamplingkernels.h"
19 :
20 : #include "gdal_vectorx.h"
21 :
22 : #include <algorithm>
23 : #include <complex>
24 :
25 : template <typename T> bool areEqualReal(double dfNoDataValue, T dfOut);
26 :
27 16750 : template <> bool areEqualReal(double dfNoDataValue, double dfOut)
28 : {
29 16750 : return ARE_REAL_EQUAL(dfNoDataValue, dfOut);
30 : }
31 :
32 76 : template <> bool areEqualReal(double dfNoDataValue, std::complex<double> dfOut)
33 : {
34 76 : return ARE_REAL_EQUAL(dfNoDataValue, dfOut.real());
35 : }
36 :
37 : // Only valid for T = double or std::complex<double>
38 : template <typename T>
39 254600 : bool GDALInterpExtractValuesWindow(GDALRasterBand *pBand,
40 : std::unique_ptr<DoublePointsCache> &cache,
41 : gdal::Vector2i point,
42 : gdal::Vector2i dimensions, T *padfOut)
43 : {
44 254600 : constexpr int BLOCK_SIZE = 64;
45 :
46 254600 : const int nX = point.x();
47 254600 : const int nY = point.y();
48 254600 : const int nWidth = dimensions.x();
49 254600 : const int nHeight = dimensions.y();
50 :
51 : // Request the DEM by blocks of BLOCK_SIZE * BLOCK_SIZE and put them
52 : // in cache
53 254600 : if (!cache)
54 100 : cache.reset(new DoublePointsCache{});
55 :
56 254600 : const int nXIters = (nX + nWidth - 1) / BLOCK_SIZE - nX / BLOCK_SIZE + 1;
57 254600 : const int nYIters = (nY + nHeight - 1) / BLOCK_SIZE - nY / BLOCK_SIZE + 1;
58 254600 : const int nRasterXSize = pBand->GetXSize();
59 254600 : const int nRasterYSize = pBand->GetYSize();
60 : const bool bIsComplex =
61 254600 : CPL_TO_BOOL(GDALDataTypeIsComplex(pBand->GetRasterDataType()));
62 :
63 509636 : for (int iY = 0; iY < nYIters; iY++)
64 : {
65 255036 : const int nBlockY = nY / BLOCK_SIZE + iY;
66 255036 : const int nReqYSize =
67 255036 : std::min(nRasterYSize - nBlockY * BLOCK_SIZE, BLOCK_SIZE);
68 255036 : const int nFirstLineInCachedBlock = (iY == 0) ? nY % BLOCK_SIZE : 0;
69 255036 : const int nFirstLineInOutput =
70 : (iY == 0) ? 0
71 436 : : BLOCK_SIZE - (nY % BLOCK_SIZE) + (iY - 1) * BLOCK_SIZE;
72 255908 : const int nLinesToCopy = (nYIters == 1) ? nHeight
73 436 : : (iY == 0) ? BLOCK_SIZE - (nY % BLOCK_SIZE)
74 436 : : (iY == nYIters - 1)
75 436 : ? 1 + (nY + nHeight - 1) % BLOCK_SIZE
76 : : BLOCK_SIZE;
77 514121 : for (int iX = 0; iX < nXIters; iX++)
78 : {
79 259085 : const int nBlockX = nX / BLOCK_SIZE + iX;
80 259085 : const int nReqXSize =
81 259085 : std::min(nRasterXSize - nBlockX * BLOCK_SIZE, BLOCK_SIZE);
82 259085 : const uint64_t nKey =
83 259085 : (static_cast<uint64_t>(nBlockY) << 32) | nBlockX;
84 259085 : const int nFirstColInCachedBlock = (iX == 0) ? nX % BLOCK_SIZE : 0;
85 259085 : const int nFirstColInOutput =
86 : (iX == 0)
87 : ? 0
88 4049 : : BLOCK_SIZE - (nX % BLOCK_SIZE) + (iX - 1) * BLOCK_SIZE;
89 267183 : const int nColsToCopy = (nXIters == 1) ? nWidth
90 4049 : : (iX == 0) ? BLOCK_SIZE - (nX % BLOCK_SIZE)
91 4049 : : (iX == nXIters - 1)
92 4049 : ? 1 + (nX + nWidth - 1) % BLOCK_SIZE
93 : : BLOCK_SIZE;
94 :
95 : #if 0
96 : CPLDebug("RPC", "nY=%d nX=%d nBlockY=%d nBlockX=%d "
97 : "nFirstLineInCachedBlock=%d nFirstLineInOutput=%d nLinesToCopy=%d "
98 : "nFirstColInCachedBlock=%d nFirstColInOutput=%d nColsToCopy=%d",
99 : nY, nX, nBlockY, nBlockX, nFirstLineInCachedBlock, nFirstLineInOutput, nLinesToCopy,
100 : nFirstColInCachedBlock, nFirstColInOutput, nColsToCopy);
101 : #endif
102 :
103 259085 : constexpr int nTypeFactor = sizeof(T) / sizeof(double);
104 0 : std::shared_ptr<std::vector<double>> poValue;
105 259085 : if (!cache->tryGet(nKey, poValue))
106 : {
107 739 : const GDALDataType eDataType =
108 : bIsComplex ? GDT_CFloat64 : GDT_Float64;
109 739 : const size_t nVectorSize =
110 739 : size_t(nReqXSize) * nReqYSize * nTypeFactor;
111 739 : poValue = std::make_shared<std::vector<double>>(nVectorSize);
112 739 : CPLErr eErr = pBand->RasterIO(
113 : GF_Read, nBlockX * BLOCK_SIZE, nBlockY * BLOCK_SIZE,
114 739 : nReqXSize, nReqYSize, poValue->data(), nReqXSize, nReqYSize,
115 : eDataType, 0, 0, nullptr);
116 739 : if (eErr != CE_None)
117 : {
118 0 : return false;
119 : }
120 739 : cache->insert(nKey, poValue);
121 : }
122 :
123 259085 : double *padfAsDouble = reinterpret_cast<double *>(padfOut);
124 : // Compose the cached block to the final buffer
125 776350 : for (int j = 0; j < nLinesToCopy; j++)
126 : {
127 517265 : const size_t dstOffset =
128 517265 : (static_cast<size_t>(nFirstLineInOutput + j) * nWidth +
129 517265 : nFirstColInOutput) *
130 : nTypeFactor;
131 517265 : const size_t srcOffset =
132 517265 : (static_cast<size_t>(nFirstLineInCachedBlock + j) *
133 517265 : nReqXSize +
134 517265 : nFirstColInCachedBlock) *
135 : nTypeFactor;
136 517265 : if (srcOffset + nColsToCopy * nTypeFactor > poValue->size())
137 : {
138 0 : return false;
139 : }
140 517265 : memcpy(padfAsDouble + dstOffset, poValue->data() + srcOffset,
141 517265 : nColsToCopy * sizeof(T));
142 : }
143 : }
144 : }
145 :
146 : #if 0
147 : CPLDebug("RPC_DEM", "DEM for %d,%d,%d,%d", nX, nY, nWidth, nHeight);
148 : for(int j = 0; j < nHeight; j++)
149 : {
150 : std::string osLine;
151 : for(int i = 0; i < nWidth; ++i )
152 : {
153 : if( !osLine.empty() )
154 : osLine += ", ";
155 : osLine += std::to_string(padfOut[j * nWidth + i]);
156 : }
157 : CPLDebug("RPC_DEM", "%s", osLine.c_str());
158 : }
159 : #endif
160 :
161 254600 : return true;
162 : }
163 :
164 : template <typename T>
165 307910 : bool GDALInterpolateAtPointImpl(GDALRasterBand *pBand,
166 : GDALRIOResampleAlg eResampleAlg,
167 : std::unique_ptr<DoublePointsCache> &cache,
168 : double dfXIn, double dfYIn, T &out)
169 : {
170 307910 : const gdal::Vector2i rasterSize{pBand->GetXSize(), pBand->GetYSize()};
171 :
172 307910 : if (eResampleAlg == GRIORA_NearestNeighbour)
173 : {
174 : // Allow input coordinates right at the bottom or right edge
175 : // with GRIORA_NearestNeighbour.
176 : // "introduce" them in the pixel of the image.
177 100 : if (dfXIn >= rasterSize.x() && dfXIn <= rasterSize.x() + 1e-5)
178 5 : dfXIn -= 0.25;
179 100 : if (dfYIn >= rasterSize.y() && dfYIn <= rasterSize.y() + 1e-5)
180 5 : dfYIn -= 0.25;
181 : }
182 307910 : const gdal::Vector2d inLoc{dfXIn, dfYIn};
183 :
184 307910 : int bGotNoDataValue = FALSE;
185 307910 : const double dfNoDataValue = pBand->GetNoDataValue(&bGotNoDataValue);
186 :
187 562530 : if (inLoc.x() < 0 || inLoc.x() > rasterSize.x() || inLoc.y() < 0 ||
188 254620 : inLoc.y() > rasterSize.y())
189 : {
190 53310 : return FALSE;
191 : }
192 :
193 : // Downgrade the interpolation algorithm if the image is too small
194 254613 : if ((rasterSize.x() < 4 || rasterSize.y() < 4) &&
195 13 : (eResampleAlg == GRIORA_CubicSpline || eResampleAlg == GRIORA_Cubic))
196 : {
197 0 : eResampleAlg = GRIORA_Bilinear;
198 : }
199 254605 : if ((rasterSize.x() < 2 || rasterSize.y() < 2) &&
200 5 : eResampleAlg == GRIORA_Bilinear)
201 : {
202 0 : eResampleAlg = GRIORA_NearestNeighbour;
203 : }
204 :
205 : auto outOfBorderCorrectionSimple =
206 509018 : [](int dNew, int nRasterSize, int nKernelsize)
207 : {
208 509018 : int dOutOfBorder = 0;
209 509018 : if (dNew < 0)
210 : {
211 34586 : dOutOfBorder = dNew;
212 : }
213 509018 : if (dNew + nKernelsize >= nRasterSize)
214 : {
215 45848 : dOutOfBorder = dNew + nKernelsize - nRasterSize;
216 : }
217 509018 : return dOutOfBorder;
218 : };
219 :
220 254600 : auto outOfBorderCorrection =
221 509018 : [&outOfBorderCorrectionSimple,
222 : &rasterSize](gdal::Vector2i input, int nKernelsize) -> gdal::Vector2i
223 : {
224 : return {
225 254509 : outOfBorderCorrectionSimple(input.x(), rasterSize.x(), nKernelsize),
226 254509 : outOfBorderCorrectionSimple(input.y(), rasterSize.y(),
227 763527 : nKernelsize)};
228 : };
229 :
230 : auto dragReadDataInBorderSimple =
231 1052638 : [](T *adfElevData, int dOutOfBorder, int nKernelSize, bool bIsX)
232 : {
233 543616 : while (dOutOfBorder < 0)
234 : {
235 103848 : for (int j = 0; j < nKernelSize; j++)
236 138716 : for (int ii = 0; ii < nKernelSize - 1; ii++)
237 : {
238 69466 : const int i = nKernelSize - ii - 2; // iterate in reverse
239 69466 : const int row_src = bIsX ? j : i;
240 69466 : const int row_dst = bIsX ? j : i + 1;
241 69466 : const int col_src = bIsX ? i : j;
242 69466 : const int col_dst = bIsX ? i + 1 : j;
243 69466 : adfElevData[nKernelSize * row_dst + col_dst] =
244 69466 : adfElevData[nKernelSize * row_src + col_src];
245 : }
246 34598 : dOutOfBorder++;
247 : }
248 543482 : while (dOutOfBorder > 0)
249 : {
250 103410 : for (int j = 0; j < nKernelSize; j++)
251 137964 : for (int i = 0; i < nKernelSize - 1; i++)
252 : {
253 69018 : const int row_src = bIsX ? j : i + 1;
254 69018 : const int row_dst = bIsX ? j : i;
255 69018 : const int col_src = bIsX ? i + 1 : j;
256 69018 : const int col_dst = bIsX ? i : j;
257 69018 : adfElevData[nKernelSize * row_dst + col_dst] =
258 69018 : adfElevData[nKernelSize * row_src + col_src];
259 : }
260 34464 : dOutOfBorder--;
261 : }
262 : };
263 1018127 : auto dragReadDataInBorder = [&dragReadDataInBorderSimple](
264 : T *adfElevData, gdal::Vector2i dOutOfBorder,
265 : int nKernelSize) -> void
266 : {
267 254509 : dragReadDataInBorderSimple(adfElevData, dOutOfBorder.x(), nKernelSize,
268 : true);
269 254509 : dragReadDataInBorderSimple(adfElevData, dOutOfBorder.y(), nKernelSize,
270 : false);
271 : };
272 :
273 525756 : auto applyBilinearKernel = [&](gdal::Vector2d dfDelta, T *adfValues,
274 : T &pdfRes) -> bool
275 : {
276 254468 : if (bGotNoDataValue)
277 : {
278 : // TODO: We could perhaps use a valid sample if there's one.
279 4172 : bool bFoundNoDataElev = false;
280 20860 : for (int k_i = 0; k_i < 4; k_i++)
281 : {
282 16688 : if (areEqualReal(dfNoDataValue, adfValues[k_i]))
283 16 : bFoundNoDataElev = true;
284 : }
285 4172 : if (bFoundNoDataElev)
286 : {
287 4 : return FALSE;
288 : }
289 : }
290 254464 : const gdal::Vector2d dfDelta1 = 1.0 - dfDelta;
291 :
292 254459 : const T dfXZ1 =
293 254464 : adfValues[0] * dfDelta1.x() + adfValues[1] * dfDelta.x();
294 254459 : const T dfXZ2 =
295 254464 : adfValues[2] * dfDelta1.x() + adfValues[3] * dfDelta.x();
296 254464 : const T dfYZ = dfXZ1 * dfDelta1.y() + dfXZ2 * dfDelta.y();
297 :
298 254464 : pdfRes = dfYZ;
299 254464 : return TRUE;
300 : };
301 :
302 256081 : auto apply4x4Kernel = [&](gdal::Vector2d dfDelta, T *adfValues,
303 : T &pdfRes) -> bool
304 : {
305 41 : T dfSumH = 0.0;
306 41 : T dfSumWeight = 0.0;
307 205 : for (int k_i = 0; k_i < 4; k_i++)
308 : {
309 : // Loop across the X axis.
310 820 : for (int k_j = 0; k_j < 4; k_j++)
311 : {
312 : // Calculate the weight for the specified pixel according
313 : // to the bicubic b-spline kernel we're using for
314 : // interpolation.
315 656 : const gdal::Vector2i dKernInd = {k_j - 1, k_i - 1};
316 656 : const gdal::Vector2d fPoint = dKernInd.cast<double>() - dfDelta;
317 656 : const double dfPixelWeight =
318 656 : eResampleAlg == GDALRIOResampleAlg::GRIORA_CubicSpline
319 320 : ? CubicSplineKernel(fPoint.x()) *
320 320 : CubicSplineKernel(fPoint.y())
321 336 : : CubicKernel(fPoint.x()) * CubicKernel(fPoint.y());
322 :
323 : // Create a sum of all values
324 : // adjusted for the pixel's calculated weight.
325 656 : const T dfElev = adfValues[k_j + k_i * 4];
326 656 : if (bGotNoDataValue && areEqualReal(dfNoDataValue, dfElev))
327 64 : continue;
328 :
329 112 : dfSumH += dfElev * dfPixelWeight;
330 592 : dfSumWeight += dfPixelWeight;
331 : }
332 : }
333 41 : if (dfSumWeight == 0.0)
334 : {
335 4 : return FALSE;
336 : }
337 :
338 7 : pdfRes = dfSumH / dfSumWeight;
339 :
340 37 : return TRUE;
341 : };
342 :
343 254600 : if (eResampleAlg == GDALRIOResampleAlg::GRIORA_CubicSpline ||
344 254580 : eResampleAlg == GDALRIOResampleAlg::GRIORA_Cubic)
345 : {
346 : // Convert from upper left corner of pixel coordinates to center of
347 : // pixel coordinates:
348 41 : const gdal::Vector2d df = inLoc - 0.5;
349 41 : const gdal::Vector2i d = df.floor().template cast<int>();
350 41 : const gdal::Vector2d delta = df - d.cast<double>();
351 41 : const gdal::Vector2i dNew = d - 1;
352 41 : const int nKernelSize = 4;
353 : const gdal::Vector2i dOutOfBorder =
354 41 : outOfBorderCorrection(dNew, nKernelSize);
355 :
356 : // CubicSpline interpolation.
357 41 : T adfReadData[16] = {0.0};
358 41 : if (!GDALInterpExtractValuesWindow(pBand, cache, dNew - dOutOfBorder,
359 : {nKernelSize, nKernelSize},
360 : adfReadData))
361 : {
362 0 : return FALSE;
363 : }
364 41 : dragReadDataInBorder(adfReadData, dOutOfBorder, nKernelSize);
365 41 : if (!apply4x4Kernel(delta, adfReadData, out))
366 4 : return FALSE;
367 :
368 37 : return TRUE;
369 : }
370 254559 : else if (eResampleAlg == GDALRIOResampleAlg::GRIORA_Bilinear)
371 : {
372 : // Convert from upper left corner of pixel coordinates to center of
373 : // pixel coordinates:
374 254468 : const gdal::Vector2d df = inLoc - 0.5;
375 254468 : const gdal::Vector2i d = df.floor().template cast<int>();
376 254468 : const gdal::Vector2d delta = df - d.cast<double>();
377 254468 : const int nKernelSize = 2;
378 : const gdal::Vector2i dOutOfBorder =
379 254468 : outOfBorderCorrection(d, nKernelSize);
380 :
381 : // Bilinear interpolation.
382 254468 : T adfReadData[4] = {0.0};
383 254468 : if (!GDALInterpExtractValuesWindow(pBand, cache, d - dOutOfBorder,
384 : {nKernelSize, nKernelSize},
385 : adfReadData))
386 : {
387 0 : return FALSE;
388 : }
389 254468 : dragReadDataInBorder(adfReadData, dOutOfBorder, nKernelSize);
390 254468 : if (!applyBilinearKernel(delta, adfReadData, out))
391 4 : return FALSE;
392 :
393 254464 : return TRUE;
394 : }
395 : else
396 : {
397 91 : const gdal::Vector2i d = inLoc.cast<int>();
398 91 : T adfOut[1] = {};
399 182 : if (!GDALInterpExtractValuesWindow(pBand, cache, d, {1, 1}, adfOut) ||
400 91 : (bGotNoDataValue && areEqualReal(dfNoDataValue, adfOut[0])))
401 : {
402 5 : return FALSE;
403 : }
404 :
405 86 : out = adfOut[0];
406 :
407 86 : return TRUE;
408 : }
409 : }
410 :
411 : /************************************************************************/
412 : /* GDALInterpolateAtPoint() */
413 : /************************************************************************/
414 :
415 307910 : bool GDALInterpolateAtPoint(GDALRasterBand *pBand,
416 : GDALRIOResampleAlg eResampleAlg,
417 : std::unique_ptr<DoublePointsCache> &cache,
418 : const double dfXIn, const double dfYIn,
419 : double *pdfOutputReal, double *pdfOutputImag)
420 : {
421 : const bool bIsComplex =
422 307910 : CPL_TO_BOOL(GDALDataTypeIsComplex(pBand->GetRasterDataType()));
423 307910 : bool res = TRUE;
424 307910 : if (bIsComplex)
425 : {
426 23 : std::complex<double> out{};
427 23 : res = GDALInterpolateAtPointImpl(pBand, eResampleAlg, cache, dfXIn,
428 : dfYIn, out);
429 23 : *pdfOutputReal = out.real();
430 23 : if (pdfOutputImag)
431 23 : *pdfOutputImag = out.imag();
432 : }
433 : else
434 : {
435 307887 : double out{};
436 307887 : res = GDALInterpolateAtPointImpl(pBand, eResampleAlg, cache, dfXIn,
437 : dfYIn, out);
438 307887 : *pdfOutputReal = out;
439 307887 : if (pdfOutputImag)
440 119 : *pdfOutputImag = 0;
441 : }
442 307910 : return res;
443 : }
|
