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 16749 : template <> bool areEqualReal(double dfNoDataValue, double dfOut)
28 : {
29 16749 : 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 253673 : bool GDALInterpExtractValuesWindow(GDALRasterBand *pBand,
40 : std::unique_ptr<DoublePointsCache> &cache,
41 : gdal::Vector2i point,
42 : gdal::Vector2i dimensions, T *padfOut)
43 : {
44 253673 : constexpr int BLOCK_SIZE = 64;
45 :
46 253673 : const int nX = point.x();
47 253673 : const int nY = point.y();
48 253673 : const int nWidth = dimensions.x();
49 253673 : const int nHeight = dimensions.y();
50 :
51 : // Request the DEM by blocks of BLOCK_SIZE * BLOCK_SIZE and put them
52 : // in cache
53 253673 : if (!cache)
54 81 : cache.reset(new DoublePointsCache{});
55 :
56 253673 : const int nXIters = (nX + nWidth - 1) / BLOCK_SIZE - nX / BLOCK_SIZE + 1;
57 253673 : const int nYIters = (nY + nHeight - 1) / BLOCK_SIZE - nY / BLOCK_SIZE + 1;
58 253673 : const int nRasterXSize = pBand->GetXSize();
59 253673 : const int nRasterYSize = pBand->GetYSize();
60 : const bool bIsComplex =
61 253673 : CPL_TO_BOOL(GDALDataTypeIsComplex(pBand->GetRasterDataType()));
62 :
63 507782 : for (int iY = 0; iY < nYIters; iY++)
64 : {
65 254109 : const int nBlockY = nY / BLOCK_SIZE + iY;
66 254109 : const int nReqYSize =
67 254109 : std::min(nRasterYSize - nBlockY * BLOCK_SIZE, BLOCK_SIZE);
68 254109 : const int nFirstLineInCachedBlock = (iY == 0) ? nY % BLOCK_SIZE : 0;
69 254109 : const int nFirstLineInOutput =
70 : (iY == 0) ? 0
71 436 : : BLOCK_SIZE - (nY % BLOCK_SIZE) + (iY - 1) * BLOCK_SIZE;
72 254981 : 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 512267 : for (int iX = 0; iX < nXIters; iX++)
78 : {
79 258158 : const int nBlockX = nX / BLOCK_SIZE + iX;
80 258158 : const int nReqXSize =
81 258158 : std::min(nRasterXSize - nBlockX * BLOCK_SIZE, BLOCK_SIZE);
82 258158 : const uint64_t nKey =
83 258158 : (static_cast<uint64_t>(nBlockY) << 32) | nBlockX;
84 258158 : const int nFirstColInCachedBlock = (iX == 0) ? nX % BLOCK_SIZE : 0;
85 258158 : const int nFirstColInOutput =
86 : (iX == 0)
87 : ? 0
88 4049 : : BLOCK_SIZE - (nX % BLOCK_SIZE) + (iX - 1) * BLOCK_SIZE;
89 266256 : 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 258158 : constexpr int nTypeFactor = sizeof(T) / sizeof(double);
104 0 : std::shared_ptr<std::vector<double>> poValue;
105 258158 : if (!cache->tryGet(nKey, poValue))
106 : {
107 720 : const GDALDataType eDataType =
108 : bIsComplex ? GDT_CFloat64 : GDT_Float64;
109 720 : const size_t nVectorSize =
110 720 : size_t(nReqXSize) * nReqYSize * nTypeFactor;
111 720 : poValue = std::make_shared<std::vector<double>>(nVectorSize);
112 720 : CPLErr eErr = pBand->RasterIO(
113 : GF_Read, nBlockX * BLOCK_SIZE, nBlockY * BLOCK_SIZE,
114 720 : nReqXSize, nReqYSize, poValue->data(), nReqXSize, nReqYSize,
115 : eDataType, 0, 0, nullptr);
116 720 : if (eErr != CE_None)
117 : {
118 0 : return false;
119 : }
120 720 : cache->insert(nKey, poValue);
121 : }
122 :
123 258158 : double *padfAsDouble = reinterpret_cast<double *>(padfOut);
124 : // Compose the cached block to the final buffer
125 773594 : for (int j = 0; j < nLinesToCopy; j++)
126 : {
127 1030868 : memcpy(padfAsDouble + ((nFirstLineInOutput + j) * nWidth +
128 515372 : nFirstColInOutput) *
129 : nTypeFactor,
130 515436 : poValue->data() +
131 515436 : ((nFirstLineInCachedBlock + j) * nReqXSize +
132 515372 : nFirstColInCachedBlock) *
133 : nTypeFactor,
134 515436 : nColsToCopy * sizeof(T));
135 : }
136 : }
137 : }
138 :
139 : #if 0
140 : CPLDebug("RPC_DEM", "DEM for %d,%d,%d,%d", nX, nY, nWidth, nHeight);
141 : for(int j = 0; j < nHeight; j++)
142 : {
143 : std::string osLine;
144 : for(int i = 0; i < nWidth; ++i )
145 : {
146 : if( !osLine.empty() )
147 : osLine += ", ";
148 : osLine += std::to_string(padfOut[j * nWidth + i]);
149 : }
150 : CPLDebug("RPC_DEM", "%s", osLine.c_str());
151 : }
152 : #endif
153 :
154 253673 : return true;
155 : }
156 :
157 : template <typename T>
158 306645 : bool GDALInterpolateAtPointImpl(GDALRasterBand *pBand,
159 : GDALRIOResampleAlg eResampleAlg,
160 : std::unique_ptr<DoublePointsCache> &cache,
161 : const double dfXIn, const double dfYIn, T &out)
162 : {
163 306645 : const gdal::Vector2i rasterSize{pBand->GetXSize(), pBand->GetYSize()};
164 306645 : const gdal::Vector2d inLoc{dfXIn, dfYIn};
165 :
166 306645 : int bGotNoDataValue = FALSE;
167 306645 : const double dfNoDataValue = pBand->GetNoDataValue(&bGotNoDataValue);
168 :
169 560338 : if (inLoc.x() < 0 || inLoc.x() > rasterSize.x() || inLoc.y() < 0 ||
170 253693 : inLoc.y() > rasterSize.y())
171 : {
172 52972 : return FALSE;
173 : }
174 :
175 : // Downgrade the interpolation algorithm if the image is too small
176 253681 : if ((rasterSize.x() < 4 || rasterSize.y() < 4) &&
177 8 : (eResampleAlg == GRIORA_CubicSpline || eResampleAlg == GRIORA_Cubic))
178 : {
179 0 : eResampleAlg = GRIORA_Bilinear;
180 : }
181 253673 : if ((rasterSize.x() < 2 || rasterSize.y() < 2) &&
182 0 : eResampleAlg == GRIORA_Bilinear)
183 : {
184 0 : eResampleAlg = GRIORA_NearestNeighbour;
185 : }
186 :
187 : auto outOfBorderCorrectionSimple =
188 507222 : [](int dNew, int nRasterSize, int nKernelsize)
189 : {
190 507222 : int dOutOfBorder = 0;
191 507222 : if (dNew < 0)
192 : {
193 34586 : dOutOfBorder = dNew;
194 : }
195 507222 : if (dNew + nKernelsize >= nRasterSize)
196 : {
197 45557 : dOutOfBorder = dNew + nKernelsize - nRasterSize;
198 : }
199 507222 : return dOutOfBorder;
200 : };
201 :
202 253673 : auto outOfBorderCorrection =
203 507222 : [&outOfBorderCorrectionSimple,
204 : &rasterSize](gdal::Vector2i input, int nKernelsize) -> gdal::Vector2i
205 : {
206 : return {
207 253611 : outOfBorderCorrectionSimple(input.x(), rasterSize.x(), nKernelsize),
208 253611 : outOfBorderCorrectionSimple(input.y(), rasterSize.y(),
209 760833 : nKernelsize)};
210 : };
211 :
212 : auto dragReadDataInBorderSimple =
213 1049038 : [](T *adfElevData, int dOutOfBorder, int nKernelSize, bool bIsX)
214 : {
215 541820 : while (dOutOfBorder < 0)
216 : {
217 103848 : for (int j = 0; j < nKernelSize; j++)
218 138716 : for (int ii = 0; ii < nKernelSize - 1; ii++)
219 : {
220 69466 : const int i = nKernelSize - ii - 2; // iterate in reverse
221 69466 : const int row_src = bIsX ? j : i;
222 69466 : const int row_dst = bIsX ? j : i + 1;
223 69466 : const int col_src = bIsX ? i : j;
224 69466 : const int col_dst = bIsX ? i + 1 : j;
225 69466 : adfElevData[nKernelSize * row_dst + col_dst] =
226 69466 : adfElevData[nKernelSize * row_src + col_src];
227 : }
228 34598 : dOutOfBorder++;
229 : }
230 541411 : while (dOutOfBorder > 0)
231 : {
232 102585 : for (int j = 0; j < nKernelSize; j++)
233 136864 : for (int i = 0; i < nKernelSize - 1; i++)
234 : {
235 68468 : const int row_src = bIsX ? j : i + 1;
236 68468 : const int row_dst = bIsX ? j : i;
237 68468 : const int col_src = bIsX ? i + 1 : j;
238 68468 : const int col_dst = bIsX ? i : j;
239 68468 : adfElevData[nKernelSize * row_dst + col_dst] =
240 68468 : adfElevData[nKernelSize * row_src + col_src];
241 : }
242 34189 : dOutOfBorder--;
243 : }
244 : };
245 1014506 : auto dragReadDataInBorder = [&dragReadDataInBorderSimple](
246 : T *adfElevData, gdal::Vector2i dOutOfBorder,
247 : int nKernelSize) -> void
248 : {
249 253611 : dragReadDataInBorderSimple(adfElevData, dOutOfBorder.x(), nKernelSize,
250 : true);
251 253611 : dragReadDataInBorderSimple(adfElevData, dOutOfBorder.y(), nKernelSize,
252 : false);
253 : };
254 :
255 523933 : auto applyBilinearKernel = [&](gdal::Vector2d dfDelta, T *adfValues,
256 : T &pdfRes) -> bool
257 : {
258 253572 : if (bGotNoDataValue)
259 : {
260 : // TODO: We could perhaps use a valid sample if there's one.
261 4172 : bool bFoundNoDataElev = false;
262 20860 : for (int k_i = 0; k_i < 4; k_i++)
263 : {
264 16688 : if (areEqualReal(dfNoDataValue, adfValues[k_i]))
265 16 : bFoundNoDataElev = true;
266 : }
267 4172 : if (bFoundNoDataElev)
268 : {
269 4 : return FALSE;
270 : }
271 : }
272 253568 : const gdal::Vector2d dfDelta1 = 1.0 - dfDelta;
273 :
274 253563 : const T dfXZ1 =
275 253568 : adfValues[0] * dfDelta1.x() + adfValues[1] * dfDelta.x();
276 253563 : const T dfXZ2 =
277 253568 : adfValues[2] * dfDelta1.x() + adfValues[3] * dfDelta.x();
278 253568 : const T dfYZ = dfXZ1 * dfDelta1.y() + dfXZ2 * dfDelta.y();
279 :
280 253568 : pdfRes = dfYZ;
281 253568 : return TRUE;
282 : };
283 :
284 255088 : auto apply4x4Kernel = [&](gdal::Vector2d dfDelta, T *adfValues,
285 : T &pdfRes) -> bool
286 : {
287 39 : T dfSumH = 0.0;
288 39 : T dfSumWeight = 0.0;
289 195 : for (int k_i = 0; k_i < 4; k_i++)
290 : {
291 : // Loop across the X axis.
292 780 : for (int k_j = 0; k_j < 4; k_j++)
293 : {
294 : // Calculate the weight for the specified pixel according
295 : // to the bicubic b-spline kernel we're using for
296 : // interpolation.
297 624 : const gdal::Vector2i dKernInd = {k_j - 1, k_i - 1};
298 624 : const gdal::Vector2d fPoint = dKernInd.cast<double>() - dfDelta;
299 624 : const double dfPixelWeight =
300 624 : eResampleAlg == GDALRIOResampleAlg::GRIORA_CubicSpline
301 304 : ? CubicSplineKernel(fPoint.x()) *
302 304 : CubicSplineKernel(fPoint.y())
303 320 : : CubicKernel(fPoint.x()) * CubicKernel(fPoint.y());
304 :
305 : // Create a sum of all values
306 : // adjusted for the pixel's calculated weight.
307 624 : const T dfElev = adfValues[k_j + k_i * 4];
308 624 : if (bGotNoDataValue && areEqualReal(dfNoDataValue, dfElev))
309 64 : continue;
310 :
311 112 : dfSumH += dfElev * dfPixelWeight;
312 560 : dfSumWeight += dfPixelWeight;
313 : }
314 : }
315 39 : if (dfSumWeight == 0.0)
316 : {
317 4 : return FALSE;
318 : }
319 :
320 7 : pdfRes = dfSumH / dfSumWeight;
321 :
322 35 : return TRUE;
323 : };
324 :
325 253673 : if (eResampleAlg == GDALRIOResampleAlg::GRIORA_CubicSpline ||
326 253654 : eResampleAlg == GDALRIOResampleAlg::GRIORA_Cubic)
327 : {
328 : // Convert from upper left corner of pixel coordinates to center of
329 : // pixel coordinates:
330 39 : const gdal::Vector2d df = inLoc - 0.5;
331 39 : const gdal::Vector2i d = df.floor().template cast<int>();
332 39 : const gdal::Vector2d delta = df - d.cast<double>();
333 39 : const gdal::Vector2i dNew = d - 1;
334 39 : const int nKernelSize = 4;
335 : const gdal::Vector2i dOutOfBorder =
336 39 : outOfBorderCorrection(dNew, nKernelSize);
337 :
338 : // CubicSpline interpolation.
339 39 : T adfReadData[16] = {0.0};
340 39 : if (!GDALInterpExtractValuesWindow(pBand, cache, dNew - dOutOfBorder,
341 : {nKernelSize, nKernelSize},
342 : adfReadData))
343 : {
344 0 : return FALSE;
345 : }
346 39 : dragReadDataInBorder(adfReadData, dOutOfBorder, nKernelSize);
347 39 : if (!apply4x4Kernel(delta, adfReadData, out))
348 4 : return FALSE;
349 :
350 35 : return TRUE;
351 : }
352 253634 : else if (eResampleAlg == GDALRIOResampleAlg::GRIORA_Bilinear)
353 : {
354 : // Convert from upper left corner of pixel coordinates to center of
355 : // pixel coordinates:
356 253572 : const gdal::Vector2d df = inLoc - 0.5;
357 253572 : const gdal::Vector2i d = df.floor().template cast<int>();
358 253572 : const gdal::Vector2d delta = df - d.cast<double>();
359 253572 : const int nKernelSize = 2;
360 : const gdal::Vector2i dOutOfBorder =
361 253572 : outOfBorderCorrection(d, nKernelSize);
362 :
363 : // Bilinear interpolation.
364 253572 : T adfReadData[4] = {0.0};
365 253572 : if (!GDALInterpExtractValuesWindow(pBand, cache, d - dOutOfBorder,
366 : {nKernelSize, nKernelSize},
367 : adfReadData))
368 : {
369 0 : return FALSE;
370 : }
371 253572 : dragReadDataInBorder(adfReadData, dOutOfBorder, nKernelSize);
372 253572 : if (!applyBilinearKernel(delta, adfReadData, out))
373 4 : return FALSE;
374 :
375 253568 : return TRUE;
376 : }
377 : else
378 : {
379 62 : const gdal::Vector2i d = inLoc.cast<int>();
380 62 : T adfOut[1] = {};
381 124 : if (!GDALInterpExtractValuesWindow(pBand, cache, d, {1, 1}, adfOut) ||
382 62 : (bGotNoDataValue && areEqualReal(dfNoDataValue, adfOut[0])))
383 : {
384 5 : return FALSE;
385 : }
386 :
387 57 : out = adfOut[0];
388 :
389 57 : return TRUE;
390 : }
391 : }
392 :
393 : /************************************************************************/
394 : /* GDALInterpolateAtPoint() */
395 : /************************************************************************/
396 :
397 306645 : bool GDALInterpolateAtPoint(GDALRasterBand *pBand,
398 : GDALRIOResampleAlg eResampleAlg,
399 : std::unique_ptr<DoublePointsCache> &cache,
400 : const double dfXIn, const double dfYIn,
401 : double *pdfOutputReal, double *pdfOutputImag)
402 : {
403 : const bool bIsComplex =
404 306645 : CPL_TO_BOOL(GDALDataTypeIsComplex(pBand->GetRasterDataType()));
405 306645 : bool res = TRUE;
406 306645 : if (bIsComplex)
407 : {
408 20 : std::complex<double> out{};
409 20 : res = GDALInterpolateAtPointImpl(pBand, eResampleAlg, cache, dfXIn,
410 : dfYIn, out);
411 20 : *pdfOutputReal = out.real();
412 20 : if (pdfOutputImag)
413 20 : *pdfOutputImag = out.imag();
414 : }
415 : else
416 : {
417 306625 : double out{};
418 306625 : res = GDALInterpolateAtPointImpl(pBand, eResampleAlg, cache, dfXIn,
419 : dfYIn, out);
420 306625 : *pdfOutputReal = out;
421 306625 : if (pdfOutputImag)
422 83 : *pdfOutputImag = 0;
423 : }
424 306645 : return res;
425 : }
|
