Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL Gridding API.
4 : * Purpose: Implementation of GDAL scattered data gridder.
5 : * Author: Andrey Kiselev, dron@ak4719.spb.edu
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2007, Andrey Kiselev <dron@ak4719.spb.edu>
9 : * Copyright (c) 2009-2013, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_port.h"
15 : #include "gdalgrid.h"
16 : #include "gdalgrid_priv.h"
17 :
18 : #include <cfloat>
19 : #include <climits>
20 : #include <cmath>
21 : #include <cstdlib>
22 : #include <cstring>
23 :
24 : #include <limits>
25 : #include <map>
26 : #include <utility>
27 : #include <algorithm>
28 :
29 : #include "cpl_conv.h"
30 : #include "cpl_cpu_features.h"
31 : #include "cpl_error.h"
32 : #include "cpl_multiproc.h"
33 : #include "cpl_progress.h"
34 : #include "cpl_quad_tree.h"
35 : #include "cpl_string.h"
36 : #include "cpl_vsi.h"
37 : #include "cpl_worker_thread_pool.h"
38 : #include "gdal.h"
39 :
40 : constexpr double TO_RADIANS = M_PI / 180.0;
41 :
42 : /************************************************************************/
43 : /* GDALGridGetPointBounds() */
44 : /************************************************************************/
45 :
46 5879260 : static void GDALGridGetPointBounds(const void *hFeature, CPLRectObj *pBounds)
47 : {
48 5879260 : const GDALGridPoint *psPoint = static_cast<const GDALGridPoint *>(hFeature);
49 5879260 : GDALGridXYArrays *psXYArrays = psPoint->psXYArrays;
50 5879260 : const int i = psPoint->i;
51 5879260 : const double dfX = psXYArrays->padfX[i];
52 5879260 : const double dfY = psXYArrays->padfY[i];
53 5879260 : pBounds->minx = dfX;
54 5879260 : pBounds->miny = dfY;
55 5879260 : pBounds->maxx = dfX;
56 5879260 : pBounds->maxy = dfY;
57 5879260 : }
58 :
59 : /************************************************************************/
60 : /* GDALGridInverseDistanceToAPower() */
61 : /************************************************************************/
62 :
63 : /**
64 : * Inverse distance to a power.
65 : *
66 : * The Inverse Distance to a Power gridding method is a weighted average
67 : * interpolator. You should supply the input arrays with the scattered data
68 : * values including coordinates of every data point and output grid geometry.
69 : * The function will compute interpolated value for the given position in
70 : * output grid.
71 : *
72 : * For every grid node the resulting value \f$Z\f$ will be calculated using
73 : * formula:
74 : *
75 : * \f[
76 : * Z=\frac{\sum_{i=1}^n{\frac{Z_i}{r_i^p}}}{\sum_{i=1}^n{\frac{1}{r_i^p}}}
77 : * \f]
78 : *
79 : * where
80 : * <ul>
81 : * <li> \f$Z_i\f$ is a known value at point \f$i\f$,
82 : * <li> \f$r_i\f$ is an Euclidean distance from the grid node
83 : * to point \f$i\f$,
84 : * <li> \f$p\f$ is a weighting power,
85 : * <li> \f$n\f$ is a total number of points in search ellipse.
86 : * </ul>
87 : *
88 : * In this method the weighting factor \f$w\f$ is
89 : *
90 : * \f[
91 : * w=\frac{1}{r^p}
92 : * \f]
93 : *
94 : * @param poOptionsIn Algorithm parameters. This should point to
95 : * GDALGridInverseDistanceToAPowerOptions object.
96 : * @param nPoints Number of elements in input arrays.
97 : * @param padfX Input array of X coordinates.
98 : * @param padfY Input array of Y coordinates.
99 : * @param padfZ Input array of Z values.
100 : * @param dfXPoint X coordinate of the point to compute.
101 : * @param dfYPoint Y coordinate of the point to compute.
102 : * @param pdfValue Pointer to variable where the computed grid node value
103 : * will be returned.
104 : * @param hExtraParamsIn extra parameters (unused)
105 : *
106 : * @return CE_None on success or CE_Failure if something goes wrong.
107 : */
108 :
109 471487 : CPLErr GDALGridInverseDistanceToAPower(const void *poOptionsIn, GUInt32 nPoints,
110 : const double *padfX, const double *padfY,
111 : const double *padfZ, double dfXPoint,
112 : double dfYPoint, double *pdfValue,
113 : CPL_UNUSED void *hExtraParamsIn)
114 : {
115 : // TODO: For optimization purposes pre-computed parameters should be moved
116 : // out of this routine to the calling function.
117 :
118 471487 : const GDALGridInverseDistanceToAPowerOptions *const poOptions =
119 : static_cast<const GDALGridInverseDistanceToAPowerOptions *>(
120 : poOptionsIn);
121 :
122 : // Pre-compute search ellipse parameters.
123 471487 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
124 471487 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
125 471487 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
126 :
127 : // Compute coefficients for coordinate system rotation.
128 471487 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
129 471487 : const bool bRotated = dfAngle != 0.0;
130 471487 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
131 471487 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
132 :
133 471487 : const double dfPowerDiv2 = poOptions->dfPower / 2;
134 471487 : const double dfSmoothing = poOptions->dfSmoothing;
135 471487 : const GUInt32 nMaxPoints = poOptions->nMaxPoints;
136 471487 : double dfNominator = 0.0;
137 471487 : double dfDenominator = 0.0;
138 471487 : GUInt32 n = 0;
139 :
140 2207140 : for (GUInt32 i = 0; i < nPoints; i++)
141 : {
142 1793720 : double dfRX = padfX[i] - dfXPoint;
143 1793720 : double dfRY = padfY[i] - dfYPoint;
144 1793720 : const double dfR2 =
145 1793720 : dfRX * dfRX + dfRY * dfRY + dfSmoothing * dfSmoothing;
146 :
147 1793720 : if (bRotated)
148 : {
149 248524 : const double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
150 248524 : const double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
151 :
152 248524 : dfRX = dfRXRotated;
153 248524 : dfRY = dfRYRotated;
154 : }
155 :
156 : // Is this point located inside the search ellipse?
157 1793720 : if (dfRadius2Square * dfRX * dfRX + dfRadius1Square * dfRY * dfRY <=
158 : dfR12Square)
159 : {
160 : // If the test point is close to the grid node, use the point
161 : // value directly as a node value to avoid singularity.
162 1266430 : if (dfR2 < 0.0000000000001)
163 : {
164 0 : *pdfValue = padfZ[i];
165 0 : return CE_None;
166 : }
167 :
168 1266430 : const double dfW = pow(dfR2, dfPowerDiv2);
169 1266430 : const double dfInvW = 1.0 / dfW;
170 1266430 : dfNominator += dfInvW * padfZ[i];
171 1266430 : dfDenominator += dfInvW;
172 1266430 : n++;
173 1266430 : if (nMaxPoints > 0 && n > nMaxPoints)
174 58061 : break;
175 : }
176 : }
177 :
178 471487 : if (n < poOptions->nMinPoints || dfDenominator == 0.0)
179 : {
180 18651 : *pdfValue = poOptions->dfNoDataValue;
181 : }
182 : else
183 : {
184 452836 : *pdfValue = dfNominator / dfDenominator;
185 : }
186 :
187 471487 : return CE_None;
188 : }
189 :
190 : /************************************************************************/
191 : /* GDALGridInverseDistanceToAPowerNearestNeighbor() */
192 : /************************************************************************/
193 :
194 : /**
195 : * Inverse distance to a power with nearest neighbor search, ideal when
196 : * max_points used.
197 : *
198 : * The Inverse Distance to a Power gridding method is a weighted average
199 : * interpolator. You should supply the input arrays with the scattered data
200 : * values including coordinates of every data point and output grid geometry.
201 : * The function will compute interpolated value for the given position in
202 : * output grid.
203 : *
204 : * For every grid node the resulting value \f$Z\f$ will be calculated using
205 : * formula for nearest matches:
206 : *
207 : * \f[
208 : * Z=\frac{\sum_{i=1}^n{\frac{Z_i}{r_i^p}}}{\sum_{i=1}^n{\frac{1}{r_i^p}}}
209 : * \f]
210 : *
211 : * where
212 : * <ul>
213 : * <li> \f$Z_i\f$ is a known value at point \f$i\f$,
214 : * <li> \f$r_i\f$ is an Euclidean distance from the grid node
215 : * to point \f$i\f$ (with an optional smoothing parameter \f$s\f$),
216 : * <li> \f$p\f$ is a weighting power,
217 : * <li> \f$n\f$ is a total number of points in search ellipse.
218 : * </ul>
219 : *
220 : * In this method the weighting factor \f$w\f$ is
221 : *
222 : * \f[
223 : * w=\frac{1}{r^p}
224 : * \f]
225 : *
226 : * @param poOptionsIn Algorithm parameters. This should point to
227 : * GDALGridInverseDistanceToAPowerNearestNeighborOptions object.
228 : * @param nPoints Number of elements in input arrays.
229 : * @param padfX Input array of X coordinates.
230 : * @param padfY Input array of Y coordinates.
231 : * @param padfZ Input array of Z values.
232 : * @param dfXPoint X coordinate of the point to compute.
233 : * @param dfYPoint Y coordinate of the point to compute.
234 : * @param pdfValue Pointer to variable where the computed grid node value
235 : * will be returned.
236 : * @param hExtraParamsIn extra parameters.
237 : *
238 : * @return CE_None on success or CE_Failure if something goes wrong.
239 : */
240 :
241 453923 : CPLErr GDALGridInverseDistanceToAPowerNearestNeighbor(
242 : const void *poOptionsIn, GUInt32 nPoints, const double *padfX,
243 : const double *padfY, const double *padfZ, double dfXPoint, double dfYPoint,
244 : double *pdfValue, void *hExtraParamsIn)
245 : {
246 453923 : CPL_IGNORE_RET_VAL(nPoints);
247 :
248 : const GDALGridInverseDistanceToAPowerNearestNeighborOptions
249 447203 : *const poOptions = static_cast<
250 : const GDALGridInverseDistanceToAPowerNearestNeighborOptions *>(
251 : poOptionsIn);
252 447203 : const double dfRadius = poOptions->dfRadius;
253 447203 : const double dfSmoothing = poOptions->dfSmoothing;
254 447203 : const double dfSmoothing2 = dfSmoothing * dfSmoothing;
255 :
256 447203 : const GUInt32 nMaxPoints = poOptions->nMaxPoints;
257 :
258 447203 : GDALGridExtraParameters *psExtraParams =
259 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
260 447203 : const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
261 447203 : CPLAssert(phQuadTree);
262 :
263 447203 : const double dfRPower2 = psExtraParams->dfRadiusPower2PreComp;
264 447203 : const double dfPowerDiv2 = psExtraParams->dfPowerDiv2PreComp;
265 :
266 896156 : std::multimap<double, double> oMapDistanceToZValues;
267 :
268 442299 : const double dfSearchRadius = dfRadius;
269 : CPLRectObj sAoi;
270 442299 : sAoi.minx = dfXPoint - dfSearchRadius;
271 442299 : sAoi.miny = dfYPoint - dfSearchRadius;
272 442299 : sAoi.maxx = dfXPoint + dfSearchRadius;
273 442299 : sAoi.maxy = dfYPoint + dfSearchRadius;
274 442299 : int nFeatureCount = 0;
275 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
276 442299 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
277 447675 : if (nFeatureCount != 0)
278 : {
279 2597810 : for (int k = 0; k < nFeatureCount; k++)
280 : {
281 2193380 : const int i = papsPoints[k]->i;
282 2193380 : const double dfRX = padfX[i] - dfXPoint;
283 2193380 : const double dfRY = padfY[i] - dfYPoint;
284 :
285 2193380 : const double dfR2 = dfRX * dfRX + dfRY * dfRY;
286 : // real distance + smoothing
287 2193380 : const double dfRsmoothed2 = dfR2 + dfSmoothing2;
288 2193380 : if (dfRsmoothed2 < 0.0000000000001)
289 : {
290 0 : *pdfValue = padfZ[i];
291 0 : CPLFree(papsPoints);
292 0 : return CE_None;
293 : }
294 : // is point within real distance?
295 2193380 : if (dfR2 <= dfRPower2)
296 : {
297 : oMapDistanceToZValues.insert(
298 2163700 : std::make_pair(dfRsmoothed2, padfZ[i]));
299 : }
300 : }
301 : }
302 405590 : CPLFree(papsPoints);
303 :
304 450906 : double dfNominator = 0.0;
305 450906 : double dfDenominator = 0.0;
306 450906 : GUInt32 n = 0;
307 :
308 : // Examine all "neighbors" within the radius (sorted by distance via the
309 : // multimap), and use the closest n points based on distance until the max
310 : // is reached.
311 1818730 : for (std::multimap<double, double>::iterator oMapDistanceToZValuesIter =
312 450906 : oMapDistanceToZValues.begin();
313 2257070 : oMapDistanceToZValuesIter != oMapDistanceToZValues.end();
314 1796770 : ++oMapDistanceToZValuesIter)
315 : {
316 1859440 : const double dfR2 = oMapDistanceToZValuesIter->first;
317 1836490 : const double dfZ = oMapDistanceToZValuesIter->second;
318 :
319 1860390 : const double dfW = pow(dfR2, dfPowerDiv2);
320 1860390 : const double dfInvW = 1.0 / dfW;
321 1860390 : dfNominator += dfInvW * dfZ;
322 1860390 : dfDenominator += dfInvW;
323 1860390 : n++;
324 1860390 : if (nMaxPoints > 0 && n >= nMaxPoints)
325 : {
326 63617 : break;
327 : }
328 : }
329 :
330 448953 : if (n < poOptions->nMinPoints || dfDenominator == 0.0)
331 : {
332 63983 : *pdfValue = poOptions->dfNoDataValue;
333 : }
334 : else
335 : {
336 384970 : *pdfValue = dfNominator / dfDenominator;
337 : }
338 :
339 448953 : return CE_None;
340 : }
341 :
342 : /************************************************************************/
343 : /* GDALGridInverseDistanceToAPowerNearestNeighborPerQuadrant() */
344 : /************************************************************************/
345 :
346 : /**
347 : * Inverse distance to a power with nearest neighbor search, with a per-quadrant
348 : * search logic.
349 : */
350 258425 : static CPLErr GDALGridInverseDistanceToAPowerNearestNeighborPerQuadrant(
351 : const void *poOptionsIn, GUInt32 /*nPoints*/, const double *padfX,
352 : const double *padfY, const double *padfZ, double dfXPoint, double dfYPoint,
353 : double *pdfValue, void *hExtraParamsIn)
354 : {
355 : const GDALGridInverseDistanceToAPowerNearestNeighborOptions
356 258425 : *const poOptions = static_cast<
357 : const GDALGridInverseDistanceToAPowerNearestNeighborOptions *>(
358 : poOptionsIn);
359 258425 : const double dfRadius = poOptions->dfRadius;
360 258425 : const double dfSmoothing = poOptions->dfSmoothing;
361 258425 : const double dfSmoothing2 = dfSmoothing * dfSmoothing;
362 :
363 258425 : const GUInt32 nMaxPoints = poOptions->nMaxPoints;
364 258425 : const GUInt32 nMinPointsPerQuadrant = poOptions->nMinPointsPerQuadrant;
365 258425 : const GUInt32 nMaxPointsPerQuadrant = poOptions->nMaxPointsPerQuadrant;
366 :
367 258425 : GDALGridExtraParameters *psExtraParams =
368 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
369 258425 : const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
370 258425 : CPLAssert(phQuadTree);
371 :
372 258425 : const double dfRPower2 = psExtraParams->dfRadiusPower2PreComp;
373 258425 : const double dfPowerDiv2 = psExtraParams->dfPowerDiv2PreComp;
374 2505980 : std::multimap<double, double> oMapDistanceToZValuesPerQuadrant[4];
375 :
376 255995 : const double dfSearchRadius = dfRadius;
377 : CPLRectObj sAoi;
378 255995 : sAoi.minx = dfXPoint - dfSearchRadius;
379 255995 : sAoi.miny = dfYPoint - dfSearchRadius;
380 255995 : sAoi.maxx = dfXPoint + dfSearchRadius;
381 255995 : sAoi.maxy = dfYPoint + dfSearchRadius;
382 255995 : int nFeatureCount = 0;
383 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
384 255995 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
385 256637 : if (nFeatureCount != 0)
386 : {
387 1466770 : for (int k = 0; k < nFeatureCount; k++)
388 : {
389 1226160 : const int i = papsPoints[k]->i;
390 1226160 : const double dfRX = padfX[i] - dfXPoint;
391 1226160 : const double dfRY = padfY[i] - dfYPoint;
392 :
393 1226160 : const double dfR2 = dfRX * dfRX + dfRY * dfRY;
394 : // real distance + smoothing
395 1226160 : const double dfRsmoothed2 = dfR2 + dfSmoothing2;
396 1226160 : if (dfRsmoothed2 < 0.0000000000001)
397 : {
398 0 : *pdfValue = padfZ[i];
399 0 : CPLFree(papsPoints);
400 0 : return CE_None;
401 : }
402 : // is point within real distance?
403 1226160 : if (dfR2 <= dfRPower2)
404 : {
405 1123460 : const int iQuadrant =
406 1123460 : ((dfRX >= 0) ? 1 : 0) | (((dfRY >= 0) ? 1 : 0) << 1);
407 : oMapDistanceToZValuesPerQuadrant[iQuadrant].insert(
408 1123460 : std::make_pair(dfRsmoothed2, padfZ[i]));
409 : }
410 : }
411 : }
412 239764 : CPLFree(papsPoints);
413 :
414 : std::multimap<double, double>::iterator aoIter[] = {
415 256426 : oMapDistanceToZValuesPerQuadrant[0].begin(),
416 250460 : oMapDistanceToZValuesPerQuadrant[1].begin(),
417 249295 : oMapDistanceToZValuesPerQuadrant[2].begin(),
418 251518 : oMapDistanceToZValuesPerQuadrant[3].begin(),
419 256426 : };
420 251780 : constexpr int ALL_QUADRANT_FLAGS = 1 + 2 + 4 + 8;
421 :
422 : // Examine all "neighbors" within the radius (sorted by distance via the
423 : // multimap), and use the closest n points based on distance until the max
424 : // is reached.
425 : // Do that by fetching the nearest point in quadrant 0, then the nearest
426 : // point in quadrant 1, 2 and 3, and starting again with the next nearest
427 : // point in quarant 0, etc.
428 251780 : int nQuadrantIterFinishedFlag = 0;
429 251780 : GUInt32 anPerQuadrant[4] = {0};
430 251780 : double dfNominator = 0.0;
431 251780 : double dfDenominator = 0.0;
432 251780 : GUInt32 n = 0;
433 251780 : for (int iQuadrant = 0; /* true */; iQuadrant = (iQuadrant + 1) % 4)
434 : {
435 2265130 : if (aoIter[iQuadrant] ==
436 5221970 : oMapDistanceToZValuesPerQuadrant[iQuadrant].end() ||
437 559521 : (nMaxPointsPerQuadrant > 0 &&
438 559521 : anPerQuadrant[iQuadrant] >= nMaxPointsPerQuadrant))
439 : {
440 1317990 : nQuadrantIterFinishedFlag |= 1 << iQuadrant;
441 1317990 : if (nQuadrantIterFinishedFlag == ALL_QUADRANT_FLAGS)
442 257751 : break;
443 1060240 : continue;
444 : }
445 :
446 1027340 : const double dfR2 = aoIter[iQuadrant]->first;
447 1009360 : const double dfZ = aoIter[iQuadrant]->second;
448 1014920 : ++aoIter[iQuadrant];
449 :
450 1005090 : const double dfW = pow(dfR2, dfPowerDiv2);
451 1005090 : const double dfInvW = 1.0 / dfW;
452 1005090 : dfNominator += dfInvW * dfZ;
453 1005090 : dfDenominator += dfInvW;
454 1005090 : n++;
455 1005090 : anPerQuadrant[iQuadrant]++;
456 1005090 : if (nMaxPoints > 0 && n >= nMaxPoints)
457 : {
458 1 : break;
459 : }
460 2065340 : }
461 :
462 257752 : if (nMinPointsPerQuadrant > 0 &&
463 128988 : (anPerQuadrant[0] < nMinPointsPerQuadrant ||
464 52349 : anPerQuadrant[1] < nMinPointsPerQuadrant ||
465 42502 : anPerQuadrant[2] < nMinPointsPerQuadrant ||
466 33357 : anPerQuadrant[3] < nMinPointsPerQuadrant))
467 : {
468 100181 : *pdfValue = poOptions->dfNoDataValue;
469 : }
470 157571 : else if (n < poOptions->nMinPoints || dfDenominator == 0.0)
471 : {
472 0 : *pdfValue = poOptions->dfNoDataValue;
473 : }
474 : else
475 : {
476 158331 : *pdfValue = dfNominator / dfDenominator;
477 : }
478 :
479 257752 : return CE_None;
480 : }
481 :
482 : /************************************************************************/
483 : /* GDALGridInverseDistanceToAPowerNoSearch() */
484 : /************************************************************************/
485 :
486 : /**
487 : * Inverse distance to a power for whole data set.
488 : *
489 : * This is somewhat optimized version of the Inverse Distance to a Power
490 : * method. It is used when the search ellips is not set. The algorithm and
491 : * parameters are the same as in GDALGridInverseDistanceToAPower(), but this
492 : * implementation works faster, because of no search.
493 : *
494 : * @see GDALGridInverseDistanceToAPower()
495 : */
496 :
497 116398 : CPLErr GDALGridInverseDistanceToAPowerNoSearch(
498 : const void *poOptionsIn, GUInt32 nPoints, const double *padfX,
499 : const double *padfY, const double *padfZ, double dfXPoint, double dfYPoint,
500 : double *pdfValue, void * /* hExtraParamsIn */)
501 : {
502 116398 : const GDALGridInverseDistanceToAPowerOptions *const poOptions =
503 : static_cast<const GDALGridInverseDistanceToAPowerOptions *>(
504 : poOptionsIn);
505 116398 : const double dfPowerDiv2 = poOptions->dfPower / 2.0;
506 116398 : const double dfSmoothing = poOptions->dfSmoothing;
507 116398 : const double dfSmoothing2 = dfSmoothing * dfSmoothing;
508 116398 : double dfNominator = 0.0;
509 116398 : double dfDenominator = 0.0;
510 116398 : const bool bPower2 = dfPowerDiv2 == 1.0;
511 :
512 116398 : GUInt32 i = 0; // Used after if.
513 116398 : if (bPower2)
514 : {
515 58297 : if (dfSmoothing2 > 0)
516 : {
517 237250 : for (i = 0; i < nPoints; i++)
518 : {
519 179986 : const double dfRX = padfX[i] - dfXPoint;
520 179986 : const double dfRY = padfY[i] - dfYPoint;
521 179986 : const double dfR2 = dfRX * dfRX + dfRY * dfRY + dfSmoothing2;
522 :
523 179986 : const double dfInvR2 = 1.0 / dfR2;
524 179986 : dfNominator += dfInvR2 * padfZ[i];
525 179986 : dfDenominator += dfInvR2;
526 : }
527 : }
528 : else
529 : {
530 27830 : for (i = 0; i < nPoints; i++)
531 : {
532 27195 : const double dfRX = padfX[i] - dfXPoint;
533 27195 : const double dfRY = padfY[i] - dfYPoint;
534 27195 : const double dfR2 = dfRX * dfRX + dfRY * dfRY;
535 :
536 : // If the test point is close to the grid node, use the point
537 : // value directly as a node value to avoid singularity.
538 27195 : if (dfR2 < 0.0000000000001)
539 : {
540 398 : break;
541 : }
542 :
543 26797 : const double dfInvR2 = 1.0 / dfR2;
544 26797 : dfNominator += dfInvR2 * padfZ[i];
545 26797 : dfDenominator += dfInvR2;
546 : }
547 : }
548 : }
549 : else
550 : {
551 302888 : for (i = 0; i < nPoints; i++)
552 : {
553 244787 : const double dfRX = padfX[i] - dfXPoint;
554 244787 : const double dfRY = padfY[i] - dfYPoint;
555 244787 : const double dfR2 = dfRX * dfRX + dfRY * dfRY + dfSmoothing2;
556 :
557 : // If the test point is close to the grid node, use the point
558 : // value directly as a node value to avoid singularity.
559 244787 : if (dfR2 < 0.0000000000001)
560 : {
561 0 : break;
562 : }
563 :
564 244787 : const double dfW = pow(dfR2, dfPowerDiv2);
565 244787 : const double dfInvW = 1.0 / dfW;
566 244787 : dfNominator += dfInvW * padfZ[i];
567 244787 : dfDenominator += dfInvW;
568 : }
569 : }
570 :
571 116398 : if (i != nPoints)
572 : {
573 397 : *pdfValue = padfZ[i];
574 : }
575 116001 : else if (dfDenominator == 0.0)
576 : {
577 0 : *pdfValue = poOptions->dfNoDataValue;
578 : }
579 : else
580 : {
581 116001 : *pdfValue = dfNominator / dfDenominator;
582 : }
583 :
584 116398 : return CE_None;
585 : }
586 :
587 : /************************************************************************/
588 : /* GDALGridMovingAverage() */
589 : /************************************************************************/
590 :
591 : /**
592 : * Moving average.
593 : *
594 : * The Moving Average is a simple data averaging algorithm. It uses a moving
595 : * window of elliptic form to search values and averages all data points
596 : * within the window. Search ellipse can be rotated by specified angle, the
597 : * center of ellipse located at the grid node. Also the minimum number of data
598 : * points to average can be set, if there are not enough points in window, the
599 : * grid node considered empty and will be filled with specified NODATA value.
600 : *
601 : * Mathematically it can be expressed with the formula:
602 : *
603 : * \f[
604 : * Z=\frac{\sum_{i=1}^n{Z_i}}{n}
605 : * \f]
606 : *
607 : * where
608 : * <ul>
609 : * <li> \f$Z\f$ is a resulting value at the grid node,
610 : * <li> \f$Z_i\f$ is a known value at point \f$i\f$,
611 : * <li> \f$n\f$ is a total number of points in search ellipse.
612 : * </ul>
613 : *
614 : * @param poOptionsIn Algorithm parameters. This should point to
615 : * GDALGridMovingAverageOptions object.
616 : * @param nPoints Number of elements in input arrays.
617 : * @param padfX Input array of X coordinates.
618 : * @param padfY Input array of Y coordinates.
619 : * @param padfZ Input array of Z values.
620 : * @param dfXPoint X coordinate of the point to compute.
621 : * @param dfYPoint Y coordinate of the point to compute.
622 : * @param pdfValue Pointer to variable where the computed grid node value
623 : * will be returned.
624 : * @param hExtraParamsIn extra parameters (unused)
625 : *
626 : * @return CE_None on success or CE_Failure if something goes wrong.
627 : */
628 :
629 625513 : CPLErr GDALGridMovingAverage(const void *poOptionsIn, GUInt32 nPoints,
630 : const double *padfX, const double *padfY,
631 : const double *padfZ, double dfXPoint,
632 : double dfYPoint, double *pdfValue,
633 : CPL_UNUSED void *hExtraParamsIn)
634 : {
635 : // TODO: For optimization purposes pre-computed parameters should be moved
636 : // out of this routine to the calling function.
637 :
638 625513 : const GDALGridMovingAverageOptions *const poOptions =
639 : static_cast<const GDALGridMovingAverageOptions *>(poOptionsIn);
640 : // Pre-compute search ellipse parameters.
641 625513 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
642 625513 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
643 : const double dfSearchRadius =
644 625513 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
645 601182 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
646 :
647 601182 : GDALGridExtraParameters *psExtraParams =
648 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
649 601182 : const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
650 :
651 : // Compute coefficients for coordinate system rotation.
652 601182 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
653 601182 : const bool bRotated = dfAngle != 0.0;
654 :
655 601182 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
656 601182 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
657 :
658 601182 : double dfAccumulator = 0.0;
659 :
660 601182 : GUInt32 n = 0; // Used after for.
661 601182 : if (phQuadTree != nullptr)
662 : {
663 : CPLRectObj sAoi;
664 1596 : sAoi.minx = dfXPoint - dfSearchRadius;
665 1596 : sAoi.miny = dfYPoint - dfSearchRadius;
666 1596 : sAoi.maxx = dfXPoint + dfSearchRadius;
667 1596 : sAoi.maxy = dfYPoint + dfSearchRadius;
668 1596 : int nFeatureCount = 0;
669 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
670 1596 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
671 1593 : if (nFeatureCount != 0)
672 : {
673 39443 : for (int k = 0; k < nFeatureCount; k++)
674 : {
675 37850 : const int i = papsPoints[k]->i;
676 37850 : const double dfRX = padfX[i] - dfXPoint;
677 37850 : const double dfRY = padfY[i] - dfYPoint;
678 :
679 37850 : if (dfRadius2Square * dfRX * dfRX +
680 37850 : dfRadius1Square * dfRY * dfRY <=
681 : dfR12Square)
682 : {
683 31259 : dfAccumulator += padfZ[i];
684 31259 : n++;
685 : }
686 : }
687 : }
688 1593 : CPLFree(papsPoints);
689 : }
690 : else
691 : {
692 3316680 : for (GUInt32 i = 0; i < nPoints; i++)
693 : {
694 2717100 : double dfRX = padfX[i] - dfXPoint;
695 2717100 : double dfRY = padfY[i] - dfYPoint;
696 :
697 2717100 : if (bRotated)
698 : {
699 333583 : const double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
700 333583 : const double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
701 :
702 333583 : dfRX = dfRXRotated;
703 333583 : dfRY = dfRYRotated;
704 : }
705 :
706 : // Is this point located inside the search ellipse?
707 2717100 : if (dfRadius2Square * dfRX * dfRX + dfRadius1Square * dfRY * dfRY <=
708 : dfR12Square)
709 : {
710 2060040 : dfAccumulator += padfZ[i];
711 2060040 : n++;
712 : }
713 : }
714 : }
715 :
716 601186 : if (n < poOptions->nMinPoints || n == 0)
717 : {
718 63465 : *pdfValue = poOptions->dfNoDataValue;
719 : }
720 : else
721 : {
722 537721 : *pdfValue = dfAccumulator / n;
723 : }
724 :
725 601186 : return CE_None;
726 : }
727 :
728 : /************************************************************************/
729 : /* GDALGridMovingAveragePerQuadrant() */
730 : /************************************************************************/
731 :
732 : /**
733 : * Moving average, with a per-quadrant search logic.
734 : */
735 192561 : static CPLErr GDALGridMovingAveragePerQuadrant(
736 : const void *poOptionsIn, GUInt32 /*nPoints*/, const double *padfX,
737 : const double *padfY, const double *padfZ, double dfXPoint, double dfYPoint,
738 : double *pdfValue, void *hExtraParamsIn)
739 : {
740 192561 : const GDALGridMovingAverageOptions *const poOptions =
741 : static_cast<const GDALGridMovingAverageOptions *>(poOptionsIn);
742 192561 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
743 192561 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
744 192561 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
745 :
746 192561 : const GUInt32 nMaxPoints = poOptions->nMaxPoints;
747 192561 : const GUInt32 nMinPointsPerQuadrant = poOptions->nMinPointsPerQuadrant;
748 192561 : const GUInt32 nMaxPointsPerQuadrant = poOptions->nMaxPointsPerQuadrant;
749 :
750 192561 : GDALGridExtraParameters *psExtraParams =
751 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
752 192561 : const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
753 192561 : CPLAssert(phQuadTree);
754 :
755 1859550 : std::multimap<double, double> oMapDistanceToZValuesPerQuadrant[4];
756 :
757 : const double dfSearchRadius =
758 190547 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
759 : CPLRectObj sAoi;
760 187990 : sAoi.minx = dfXPoint - dfSearchRadius;
761 187990 : sAoi.miny = dfYPoint - dfSearchRadius;
762 187990 : sAoi.maxx = dfXPoint + dfSearchRadius;
763 187990 : sAoi.maxy = dfYPoint + dfSearchRadius;
764 187990 : int nFeatureCount = 0;
765 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
766 187990 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
767 189687 : if (nFeatureCount != 0)
768 : {
769 1094060 : for (int k = 0; k < nFeatureCount; k++)
770 : {
771 924972 : const int i = papsPoints[k]->i;
772 924972 : const double dfRX = padfX[i] - dfXPoint;
773 924972 : const double dfRY = padfY[i] - dfYPoint;
774 924972 : const double dfRXSquare = dfRX * dfRX;
775 924972 : const double dfRYSquare = dfRY * dfRY;
776 :
777 924972 : if (dfRadius2Square * dfRXSquare + dfRadius1Square * dfRYSquare <=
778 : dfR12Square)
779 : {
780 935394 : const int iQuadrant =
781 935394 : ((dfRX >= 0) ? 1 : 0) | (((dfRY >= 0) ? 1 : 0) << 1);
782 : oMapDistanceToZValuesPerQuadrant[iQuadrant].insert(
783 935394 : std::make_pair(dfRXSquare + dfRYSquare, padfZ[i]));
784 : }
785 : }
786 : }
787 170302 : CPLFree(papsPoints);
788 :
789 : std::multimap<double, double>::iterator aoIter[] = {
790 191446 : oMapDistanceToZValuesPerQuadrant[0].begin(),
791 184017 : oMapDistanceToZValuesPerQuadrant[1].begin(),
792 184144 : oMapDistanceToZValuesPerQuadrant[2].begin(),
793 185234 : oMapDistanceToZValuesPerQuadrant[3].begin(),
794 191446 : };
795 182029 : constexpr int ALL_QUADRANT_FLAGS = 1 + 2 + 4 + 8;
796 :
797 : // Examine all "neighbors" within the radius (sorted by distance via the
798 : // multimap), and use the closest n points based on distance until the max
799 : // is reached.
800 : // Do that by fetching the nearest point in quadrant 0, then the nearest
801 : // point in quadrant 1, 2 and 3, and starting again with the next nearest
802 : // point in quarant 0, etc.
803 182029 : int nQuadrantIterFinishedFlag = 0;
804 182029 : GUInt32 anPerQuadrant[4] = {0};
805 182029 : double dfNominator = 0.0;
806 182029 : GUInt32 n = 0;
807 182029 : for (int iQuadrant = 0; /* true */; iQuadrant = (iQuadrant + 1) % 4)
808 : {
809 900173 : if (aoIter[iQuadrant] ==
810 2122500 : oMapDistanceToZValuesPerQuadrant[iQuadrant].end() ||
811 302643 : (nMaxPointsPerQuadrant > 0 &&
812 302643 : anPerQuadrant[iQuadrant] >= nMaxPointsPerQuadrant))
813 : {
814 243331 : nQuadrantIterFinishedFlag |= 1 << iQuadrant;
815 243331 : if (nQuadrantIterFinishedFlag == ALL_QUADRANT_FLAGS)
816 64036 : break;
817 179295 : continue;
818 : }
819 :
820 669008 : const double dfZ = aoIter[iQuadrant]->second;
821 660890 : ++aoIter[iQuadrant];
822 :
823 670885 : dfNominator += dfZ;
824 670885 : n++;
825 670885 : anPerQuadrant[iQuadrant]++;
826 670885 : if (nMaxPoints > 0 && n >= nMaxPoints)
827 : {
828 124687 : break;
829 : }
830 725493 : }
831 :
832 188723 : if (nMinPointsPerQuadrant > 0 &&
833 124996 : (anPerQuadrant[0] < nMinPointsPerQuadrant ||
834 125373 : anPerQuadrant[1] < nMinPointsPerQuadrant ||
835 125229 : anPerQuadrant[2] < nMinPointsPerQuadrant ||
836 125018 : anPerQuadrant[3] < nMinPointsPerQuadrant))
837 : {
838 61443 : *pdfValue = poOptions->dfNoDataValue;
839 : }
840 127280 : else if (n < poOptions->nMinPoints || n == 0)
841 : {
842 0 : *pdfValue = poOptions->dfNoDataValue;
843 : }
844 : else
845 : {
846 128999 : *pdfValue = dfNominator / n;
847 : }
848 :
849 378818 : return CE_None;
850 : }
851 :
852 : /************************************************************************/
853 : /* GDALGridNearestNeighbor() */
854 : /************************************************************************/
855 :
856 : /**
857 : * Nearest neighbor.
858 : *
859 : * The Nearest Neighbor method doesn't perform any interpolation or smoothing,
860 : * it just takes the value of nearest point found in grid node search ellipse
861 : * and returns it as a result. If there are no points found, the specified
862 : * NODATA value will be returned.
863 : *
864 : * @param poOptionsIn Algorithm parameters. This should point to
865 : * GDALGridNearestNeighborOptions object.
866 : * @param nPoints Number of elements in input arrays.
867 : * @param padfX Input array of X coordinates.
868 : * @param padfY Input array of Y coordinates.
869 : * @param padfZ Input array of Z values.
870 : * @param dfXPoint X coordinate of the point to compute.
871 : * @param dfYPoint Y coordinate of the point to compute.
872 : * @param pdfValue Pointer to variable where the computed grid node value
873 : * will be returned.
874 : * @param hExtraParamsIn extra parameters.
875 : *
876 : * @return CE_None on success or CE_Failure if something goes wrong.
877 : */
878 :
879 474457 : CPLErr GDALGridNearestNeighbor(const void *poOptionsIn, GUInt32 nPoints,
880 : const double *padfX, const double *padfY,
881 : const double *padfZ, double dfXPoint,
882 : double dfYPoint, double *pdfValue,
883 : void *hExtraParamsIn)
884 : {
885 : // TODO: For optimization purposes pre-computed parameters should be moved
886 : // out of this routine to the calling function.
887 :
888 474457 : const GDALGridNearestNeighborOptions *const poOptions =
889 : static_cast<const GDALGridNearestNeighborOptions *>(poOptionsIn);
890 : // Pre-compute search ellipse parameters.
891 474457 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
892 474457 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
893 474457 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
894 474457 : GDALGridExtraParameters *psExtraParams =
895 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
896 474457 : CPLQuadTree *hQuadTree = psExtraParams->hQuadTree;
897 :
898 : // Compute coefficients for coordinate system rotation.
899 474457 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
900 474457 : const bool bRotated = dfAngle != 0.0;
901 474457 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
902 474457 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
903 :
904 : // If the nearest point will not be found, its value remains as NODATA.
905 474457 : double dfNearestValue = poOptions->dfNoDataValue;
906 474457 : GUInt32 i = 0;
907 :
908 474457 : double dfSearchRadius = psExtraParams->dfInitialSearchRadius;
909 474457 : if (hQuadTree != nullptr)
910 : {
911 141829 : if (poOptions->dfRadius1 > 0 || poOptions->dfRadius2 > 0)
912 66912 : dfSearchRadius =
913 66960 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
914 : CPLRectObj sAoi;
915 350457 : while (dfSearchRadius > 0)
916 : {
917 344336 : sAoi.minx = dfXPoint - dfSearchRadius;
918 344336 : sAoi.miny = dfYPoint - dfSearchRadius;
919 344336 : sAoi.maxx = dfXPoint + dfSearchRadius;
920 344336 : sAoi.maxy = dfYPoint + dfSearchRadius;
921 344336 : int nFeatureCount = 0;
922 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
923 344336 : CPLQuadTreeSearch(hQuadTree, &sAoi, &nFeatureCount));
924 347950 : if (nFeatureCount != 0)
925 : {
926 : // Nearest distance will be initialized with the distance to the
927 : // first point in array.
928 77257 : double dfNearestRSquare = std::numeric_limits<double>::max();
929 371918 : for (int k = 0; k < nFeatureCount; k++)
930 : {
931 294661 : const int idx = papsPoints[k]->i;
932 294661 : const double dfRX = padfX[idx] - dfXPoint;
933 294661 : const double dfRY = padfY[idx] - dfYPoint;
934 :
935 294661 : const double dfR2 = dfRX * dfRX + dfRY * dfRY;
936 294661 : if (dfR2 <= dfNearestRSquare)
937 : {
938 174213 : dfNearestRSquare = dfR2;
939 174213 : dfNearestValue = padfZ[idx];
940 : }
941 : }
942 :
943 77257 : CPLFree(papsPoints);
944 131280 : break;
945 : }
946 :
947 270693 : CPLFree(papsPoints);
948 272672 : if (poOptions->dfRadius1 > 0 || poOptions->dfRadius2 > 0)
949 : break;
950 208676 : dfSearchRadius *= 2;
951 : #if DEBUG_VERBOSE
952 : CPLDebug("GDAL_GRID", "Increasing search radius to %.16g",
953 : dfSearchRadius);
954 : #endif
955 : }
956 : }
957 : else
958 : {
959 332628 : double dfNearestRSquare = std::numeric_limits<double>::max();
960 155751000 : while (i < nPoints)
961 : {
962 155419000 : double dfRX = padfX[i] - dfXPoint;
963 155419000 : double dfRY = padfY[i] - dfYPoint;
964 :
965 155419000 : if (bRotated)
966 : {
967 0 : const double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
968 0 : const double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
969 :
970 0 : dfRX = dfRXRotated;
971 0 : dfRY = dfRYRotated;
972 : }
973 :
974 : // Is this point located inside the search ellipse?
975 155419000 : const double dfRXSquare = dfRX * dfRX;
976 155419000 : const double dfRYSquare = dfRY * dfRY;
977 155419000 : if (dfRadius2Square * dfRXSquare + dfRadius1Square * dfRYSquare <=
978 : dfR12Square)
979 : {
980 156270000 : const double dfR2 = dfRXSquare + dfRYSquare;
981 156270000 : if (dfR2 <= dfNearestRSquare)
982 : {
983 16001300 : dfNearestRSquare = dfR2;
984 16001300 : dfNearestValue = padfZ[i];
985 : }
986 : }
987 :
988 155419000 : i++;
989 : }
990 : }
991 :
992 470029 : *pdfValue = dfNearestValue;
993 :
994 470029 : return CE_None;
995 : }
996 :
997 : /************************************************************************/
998 : /* GDALGridDataMetricMinimum() */
999 : /************************************************************************/
1000 :
1001 : /**
1002 : * Minimum data value (data metric).
1003 : *
1004 : * Minimum value found in grid node search ellipse. If there are no points
1005 : * found, the specified NODATA value will be returned.
1006 : *
1007 : * \f[
1008 : * Z=\min{(Z_1,Z_2,\ldots,Z_n)}
1009 : * \f]
1010 : *
1011 : * where
1012 : * <ul>
1013 : * <li> \f$Z\f$ is a resulting value at the grid node,
1014 : * <li> \f$Z_i\f$ is a known value at point \f$i\f$,
1015 : * <li> \f$n\f$ is a total number of points in search ellipse.
1016 : * </ul>
1017 : *
1018 : * @param poOptionsIn Algorithm parameters. This should point to
1019 : * GDALGridDataMetricsOptions object.
1020 : * @param nPoints Number of elements in input arrays.
1021 : * @param padfX Input array of X coordinates.
1022 : * @param padfY Input array of Y coordinates.
1023 : * @param padfZ Input array of Z values.
1024 : * @param dfXPoint X coordinate of the point to compute.
1025 : * @param dfYPoint Y coordinate of the point to compute.
1026 : * @param pdfValue Pointer to variable where the computed grid node value
1027 : * will be returned.
1028 : * @param hExtraParamsIn unused.
1029 : *
1030 : * @return CE_None on success or CE_Failure if something goes wrong.
1031 : */
1032 :
1033 311388 : CPLErr GDALGridDataMetricMinimum(const void *poOptionsIn, GUInt32 nPoints,
1034 : const double *padfX, const double *padfY,
1035 : const double *padfZ, double dfXPoint,
1036 : double dfYPoint, double *pdfValue,
1037 : void *hExtraParamsIn)
1038 : {
1039 : // TODO: For optimization purposes pre-computed parameters should be moved
1040 : // out of this routine to the calling function.
1041 :
1042 311388 : const GDALGridDataMetricsOptions *const poOptions =
1043 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
1044 :
1045 : // Pre-compute search ellipse parameters.
1046 311388 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
1047 311388 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
1048 : const double dfSearchRadius =
1049 311388 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
1050 313858 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
1051 :
1052 313858 : GDALGridExtraParameters *psExtraParams =
1053 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
1054 313858 : CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
1055 :
1056 : // Compute coefficients for coordinate system rotation.
1057 313858 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
1058 313858 : const bool bRotated = dfAngle != 0.0;
1059 313858 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
1060 313858 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
1061 :
1062 313858 : double dfMinimumValue = std::numeric_limits<double>::max();
1063 313858 : GUInt32 n = 0;
1064 313858 : if (phQuadTree != nullptr)
1065 : {
1066 : CPLRectObj sAoi;
1067 1592 : sAoi.minx = dfXPoint - dfSearchRadius;
1068 1592 : sAoi.miny = dfYPoint - dfSearchRadius;
1069 1592 : sAoi.maxx = dfXPoint + dfSearchRadius;
1070 1592 : sAoi.maxy = dfYPoint + dfSearchRadius;
1071 1592 : int nFeatureCount = 0;
1072 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
1073 1592 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
1074 1597 : if (nFeatureCount != 0)
1075 : {
1076 33717 : for (int k = 0; k < nFeatureCount; k++)
1077 : {
1078 32120 : const int i = papsPoints[k]->i;
1079 32120 : const double dfRX = padfX[i] - dfXPoint;
1080 32120 : const double dfRY = padfY[i] - dfYPoint;
1081 :
1082 32120 : if (dfRadius2Square * dfRX * dfRX +
1083 32120 : dfRadius1Square * dfRY * dfRY <=
1084 : dfR12Square)
1085 : {
1086 20573 : if (dfMinimumValue > padfZ[i])
1087 3268 : dfMinimumValue = padfZ[i];
1088 20573 : n++;
1089 : }
1090 : }
1091 : }
1092 1597 : CPLFree(papsPoints);
1093 : }
1094 : else
1095 : {
1096 312266 : GUInt32 i = 0;
1097 1697660 : while (i < nPoints)
1098 : {
1099 1385390 : double dfRX = padfX[i] - dfXPoint;
1100 1385390 : double dfRY = padfY[i] - dfYPoint;
1101 :
1102 1385390 : if (bRotated)
1103 : {
1104 70671 : const double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
1105 70671 : const double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
1106 :
1107 70671 : dfRX = dfRXRotated;
1108 70671 : dfRY = dfRYRotated;
1109 : }
1110 :
1111 : // Is this point located inside the search ellipse?
1112 1385390 : if (dfRadius2Square * dfRX * dfRX + dfRadius1Square * dfRY * dfRY <=
1113 : dfR12Square)
1114 : {
1115 576082 : if (dfMinimumValue > padfZ[i])
1116 368864 : dfMinimumValue = padfZ[i];
1117 576082 : n++;
1118 : }
1119 :
1120 1385390 : i++;
1121 : }
1122 : }
1123 :
1124 313863 : if (n < poOptions->nMinPoints || n == 0)
1125 : {
1126 76428 : *pdfValue = poOptions->dfNoDataValue;
1127 : }
1128 : else
1129 : {
1130 237435 : *pdfValue = dfMinimumValue;
1131 : }
1132 :
1133 313863 : return CE_None;
1134 : }
1135 :
1136 : /************************************************************************/
1137 : /* GDALGridDataMetricMinimumOrMaximumPerQuadrant() */
1138 : /************************************************************************/
1139 :
1140 : /**
1141 : * Minimum or maximum data value (data metric), with a per-quadrant search
1142 : * logic.
1143 : */
1144 : template <bool IS_MIN>
1145 129160 : static CPLErr GDALGridDataMetricMinimumOrMaximumPerQuadrant(
1146 : const void *poOptionsIn, const double *padfX, const double *padfY,
1147 : const double *padfZ, double dfXPoint, double dfYPoint, double *pdfValue,
1148 : void *hExtraParamsIn)
1149 : {
1150 129160 : const GDALGridDataMetricsOptions *const poOptions =
1151 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
1152 :
1153 : // Pre-compute search ellipse parameters.
1154 129160 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
1155 129160 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
1156 127509 : const double dfSearchRadius =
1157 129160 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
1158 127509 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
1159 :
1160 : // const GUInt32 nMaxPoints = poOptions->nMaxPoints;
1161 127509 : const GUInt32 nMinPointsPerQuadrant = poOptions->nMinPointsPerQuadrant;
1162 127509 : const GUInt32 nMaxPointsPerQuadrant = poOptions->nMaxPointsPerQuadrant;
1163 :
1164 127509 : GDALGridExtraParameters *psExtraParams =
1165 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
1166 127509 : const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
1167 127509 : CPLAssert(phQuadTree);
1168 :
1169 : CPLRectObj sAoi;
1170 127509 : sAoi.minx = dfXPoint - dfSearchRadius;
1171 127509 : sAoi.miny = dfYPoint - dfSearchRadius;
1172 127509 : sAoi.maxx = dfXPoint + dfSearchRadius;
1173 127509 : sAoi.maxy = dfYPoint + dfSearchRadius;
1174 127509 : int nFeatureCount = 0;
1175 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
1176 127509 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
1177 1255950 : std::multimap<double, double> oMapDistanceToZValuesPerQuadrant[4];
1178 :
1179 128172 : if (nFeatureCount != 0)
1180 : {
1181 524027 : for (int k = 0; k < nFeatureCount; k++)
1182 : {
1183 401838 : const int i = papsPoints[k]->i;
1184 401838 : const double dfRX = padfX[i] - dfXPoint;
1185 401838 : const double dfRY = padfY[i] - dfYPoint;
1186 401838 : const double dfRXSquare = dfRX * dfRX;
1187 401838 : const double dfRYSquare = dfRY * dfRY;
1188 :
1189 401838 : if (dfRadius2Square * dfRXSquare + dfRadius1Square * dfRYSquare <=
1190 : dfR12Square)
1191 : {
1192 387818 : const int iQuadrant =
1193 387818 : ((dfRX >= 0) ? 1 : 0) | (((dfRY >= 0) ? 1 : 0) << 1);
1194 381341 : oMapDistanceToZValuesPerQuadrant[iQuadrant].insert(
1195 769190 : std::make_pair(dfRXSquare + dfRYSquare, padfZ[i]));
1196 : }
1197 : }
1198 : }
1199 121726 : CPLFree(papsPoints);
1200 :
1201 626351 : std::multimap<double, double>::iterator aoIter[] = {
1202 127501 : oMapDistanceToZValuesPerQuadrant[0].begin(),
1203 124361 : oMapDistanceToZValuesPerQuadrant[1].begin(),
1204 124221 : oMapDistanceToZValuesPerQuadrant[2].begin(),
1205 124948 : oMapDistanceToZValuesPerQuadrant[3].begin(),
1206 : };
1207 125320 : constexpr int ALL_QUADRANT_FLAGS = 1 + 2 + 4 + 8;
1208 :
1209 : // Examine all "neighbors" within the radius (sorted by distance via the
1210 : // multimap), and use the closest n points based on distance until the max
1211 : // is reached.
1212 : // Do that by fetching the nearest point in quadrant 0, then the nearest
1213 : // point in quadrant 1, 2 and 3, and starting again with the next nearest
1214 : // point in quarant 0, etc.
1215 125320 : int nQuadrantIterFinishedFlag = 0;
1216 125320 : GUInt32 anPerQuadrant[4] = {0};
1217 125320 : double dfExtremum = IS_MIN ? std::numeric_limits<double>::max()
1218 : : -std::numeric_limits<double>::max();
1219 125320 : GUInt32 n = 0;
1220 903815 : for (int iQuadrant = 0; /* true */; iQuadrant = (iQuadrant + 1) % 4)
1221 : {
1222 861975 : if (aoIter[iQuadrant] ==
1223 2128850 : oMapDistanceToZValuesPerQuadrant[iQuadrant].end() ||
1224 314152 : (nMaxPointsPerQuadrant > 0 &&
1225 314152 : anPerQuadrant[iQuadrant] >= nMaxPointsPerQuadrant))
1226 : {
1227 590948 : nQuadrantIterFinishedFlag |= 1 << iQuadrant;
1228 590948 : if (nQuadrantIterFinishedFlag == ALL_QUADRANT_FLAGS)
1229 128410 : break;
1230 462538 : continue;
1231 : }
1232 :
1233 319938 : const double dfZ = aoIter[iQuadrant]->second;
1234 316785 : ++aoIter[iQuadrant];
1235 :
1236 : if (IS_MIN)
1237 : {
1238 315941 : if (dfExtremum > dfZ)
1239 243288 : dfExtremum = dfZ;
1240 : }
1241 : else
1242 : {
1243 16 : if (dfExtremum < dfZ)
1244 6 : dfExtremum = dfZ;
1245 : }
1246 315957 : n++;
1247 315957 : anPerQuadrant[iQuadrant]++;
1248 : /*if( nMaxPoints > 0 && n >= nMaxPoints )
1249 : {
1250 : break;
1251 : }*/
1252 : }
1253 :
1254 128410 : if (nMinPointsPerQuadrant > 0 &&
1255 63619 : (anPerQuadrant[0] < nMinPointsPerQuadrant ||
1256 8 : anPerQuadrant[1] < nMinPointsPerQuadrant ||
1257 6 : anPerQuadrant[2] < nMinPointsPerQuadrant ||
1258 6 : anPerQuadrant[3] < nMinPointsPerQuadrant))
1259 : {
1260 63953 : *pdfValue = poOptions->dfNoDataValue;
1261 : }
1262 64457 : else if (n < poOptions->nMinPoints || n == 0)
1263 : {
1264 1 : *pdfValue = poOptions->dfNoDataValue;
1265 : }
1266 : else
1267 : {
1268 64941 : *pdfValue = dfExtremum;
1269 : }
1270 :
1271 256732 : return CE_None;
1272 : }
1273 :
1274 : /************************************************************************/
1275 : /* GDALGridDataMetricMinimumPerQuadrant() */
1276 : /************************************************************************/
1277 :
1278 : /**
1279 : * Minimum data value (data metric), with a per-quadrant search logic.
1280 : */
1281 128857 : static CPLErr GDALGridDataMetricMinimumPerQuadrant(
1282 : const void *poOptionsIn, GUInt32 /* nPoints */, const double *padfX,
1283 : const double *padfY, const double *padfZ, double dfXPoint, double dfYPoint,
1284 : double *pdfValue, void *hExtraParamsIn)
1285 : {
1286 128857 : return GDALGridDataMetricMinimumOrMaximumPerQuadrant</*IS_MIN=*/true>(
1287 : poOptionsIn, padfX, padfY, padfZ, dfXPoint, dfYPoint, pdfValue,
1288 128965 : hExtraParamsIn);
1289 : }
1290 :
1291 : /************************************************************************/
1292 : /* GDALGridDataMetricMaximum() */
1293 : /************************************************************************/
1294 :
1295 : /**
1296 : * Maximum data value (data metric).
1297 : *
1298 : * Maximum value found in grid node search ellipse. If there are no points
1299 : * found, the specified NODATA value will be returned.
1300 : *
1301 : * \f[
1302 : * Z=\max{(Z_1,Z_2,\ldots,Z_n)}
1303 : * \f]
1304 : *
1305 : * where
1306 : * <ul>
1307 : * <li> \f$Z\f$ is a resulting value at the grid node,
1308 : * <li> \f$Z_i\f$ is a known value at point \f$i\f$,
1309 : * <li> \f$n\f$ is a total number of points in search ellipse.
1310 : * </ul>
1311 : *
1312 : * @param poOptionsIn Algorithm parameters. This should point to
1313 : * GDALGridDataMetricsOptions object.
1314 : * @param nPoints Number of elements in input arrays.
1315 : * @param padfX Input array of X coordinates.
1316 : * @param padfY Input array of Y coordinates.
1317 : * @param padfZ Input array of Z values.
1318 : * @param dfXPoint X coordinate of the point to compute.
1319 : * @param dfYPoint Y coordinate of the point to compute.
1320 : * @param pdfValue Pointer to variable where the computed grid node value
1321 : * will be returned.
1322 : * @param hExtraParamsIn extra parameters (unused)
1323 : *
1324 : * @return CE_None on success or CE_Failure if something goes wrong.
1325 : */
1326 :
1327 66605 : CPLErr GDALGridDataMetricMaximum(const void *poOptionsIn, GUInt32 nPoints,
1328 : const double *padfX, const double *padfY,
1329 : const double *padfZ, double dfXPoint,
1330 : double dfYPoint, double *pdfValue,
1331 : void *hExtraParamsIn)
1332 : {
1333 : // TODO: For optimization purposes pre-computed parameters should be moved
1334 : // out of this routine to the calling function.
1335 :
1336 66605 : const GDALGridDataMetricsOptions *const poOptions =
1337 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
1338 :
1339 : // Pre-compute search ellipse parameters.
1340 66605 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
1341 66605 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
1342 : const double dfSearchRadius =
1343 66605 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
1344 66429 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
1345 :
1346 66429 : GDALGridExtraParameters *psExtraParams =
1347 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
1348 66429 : CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
1349 :
1350 : // Compute coefficients for coordinate system rotation.
1351 66429 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
1352 66429 : const bool bRotated = dfAngle != 0.0;
1353 66429 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
1354 66429 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
1355 :
1356 66429 : double dfMaximumValue = -std::numeric_limits<double>::max();
1357 66429 : GUInt32 n = 0;
1358 66429 : if (phQuadTree != nullptr)
1359 : {
1360 : CPLRectObj sAoi;
1361 1197 : sAoi.minx = dfXPoint - dfSearchRadius;
1362 1197 : sAoi.miny = dfYPoint - dfSearchRadius;
1363 1197 : sAoi.maxx = dfXPoint + dfSearchRadius;
1364 1197 : sAoi.maxy = dfYPoint + dfSearchRadius;
1365 1197 : int nFeatureCount = 0;
1366 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
1367 1197 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
1368 1198 : if (nFeatureCount != 0)
1369 : {
1370 35970 : for (int k = 0; k < nFeatureCount; k++)
1371 : {
1372 34772 : const int i = papsPoints[k]->i;
1373 34772 : const double dfRX = padfX[i] - dfXPoint;
1374 34772 : const double dfRY = padfY[i] - dfYPoint;
1375 :
1376 34772 : if (dfRadius2Square * dfRX * dfRX +
1377 34772 : dfRadius1Square * dfRY * dfRY <=
1378 : dfR12Square)
1379 : {
1380 23139 : if (dfMaximumValue < padfZ[i])
1381 3320 : dfMaximumValue = padfZ[i];
1382 23139 : n++;
1383 : }
1384 : }
1385 : }
1386 1198 : CPLFree(papsPoints);
1387 : }
1388 : else
1389 : {
1390 65232 : GUInt32 i = 0;
1391 653956 : while (i < nPoints)
1392 : {
1393 588724 : double dfRX = padfX[i] - dfXPoint;
1394 588724 : double dfRY = padfY[i] - dfYPoint;
1395 :
1396 588724 : if (bRotated)
1397 : {
1398 213354 : const double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
1399 213354 : const double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
1400 :
1401 213354 : dfRX = dfRXRotated;
1402 213354 : dfRY = dfRYRotated;
1403 : }
1404 :
1405 : // Is this point located inside the search ellipse?
1406 588724 : if (dfRadius2Square * dfRX * dfRX + dfRadius1Square * dfRY * dfRY <=
1407 : dfR12Square)
1408 : {
1409 365123 : if (dfMaximumValue < padfZ[i])
1410 122226 : dfMaximumValue = padfZ[i];
1411 365123 : n++;
1412 : }
1413 :
1414 588724 : i++;
1415 : }
1416 : }
1417 :
1418 66431 : if (n < poOptions->nMinPoints || n == 0)
1419 : {
1420 615 : *pdfValue = poOptions->dfNoDataValue;
1421 : }
1422 : else
1423 : {
1424 65816 : *pdfValue = dfMaximumValue;
1425 : }
1426 :
1427 66431 : return CE_None;
1428 : }
1429 :
1430 : /************************************************************************/
1431 : /* GDALGridDataMetricMaximumPerQuadrant() */
1432 : /************************************************************************/
1433 :
1434 : /**
1435 : * Maximum data value (data metric), with a per-quadrant search logic.
1436 : */
1437 5 : static CPLErr GDALGridDataMetricMaximumPerQuadrant(
1438 : const void *poOptionsIn, GUInt32 /* nPoints */, const double *padfX,
1439 : const double *padfY, const double *padfZ, double dfXPoint, double dfYPoint,
1440 : double *pdfValue, void *hExtraParamsIn)
1441 : {
1442 5 : return GDALGridDataMetricMinimumOrMaximumPerQuadrant</*IS_MIN=*/false>(
1443 : poOptionsIn, padfX, padfY, padfZ, dfXPoint, dfYPoint, pdfValue,
1444 5 : hExtraParamsIn);
1445 : }
1446 :
1447 : /************************************************************************/
1448 : /* GDALGridDataMetricRange() */
1449 : /************************************************************************/
1450 :
1451 : /**
1452 : * Data range (data metric).
1453 : *
1454 : * A difference between the minimum and maximum values found in grid node
1455 : * search ellipse. If there are no points found, the specified NODATA
1456 : * value will be returned.
1457 : *
1458 : * \f[
1459 : * Z=\max{(Z_1,Z_2,\ldots,Z_n)}-\min{(Z_1,Z_2,\ldots,Z_n)}
1460 : * \f]
1461 : *
1462 : * where
1463 : * <ul>
1464 : * <li> \f$Z\f$ is a resulting value at the grid node,
1465 : * <li> \f$Z_i\f$ is a known value at point \f$i\f$,
1466 : * <li> \f$n\f$ is a total number of points in search ellipse.
1467 : * </ul>
1468 : *
1469 : * @param poOptionsIn Algorithm parameters. This should point to
1470 : * GDALGridDataMetricsOptions object.
1471 : * @param nPoints Number of elements in input arrays.
1472 : * @param padfX Input array of X coordinates.
1473 : * @param padfY Input array of Y coordinates.
1474 : * @param padfZ Input array of Z values.
1475 : * @param dfXPoint X coordinate of the point to compute.
1476 : * @param dfYPoint Y coordinate of the point to compute.
1477 : * @param pdfValue Pointer to variable where the computed grid node value
1478 : * will be returned.
1479 : * @param hExtraParamsIn extra parameters (unused)
1480 : *
1481 : * @return CE_None on success or CE_Failure if something goes wrong.
1482 : */
1483 :
1484 65841 : CPLErr GDALGridDataMetricRange(const void *poOptionsIn, GUInt32 nPoints,
1485 : const double *padfX, const double *padfY,
1486 : const double *padfZ, double dfXPoint,
1487 : double dfYPoint, double *pdfValue,
1488 : void *hExtraParamsIn)
1489 : {
1490 : // TODO: For optimization purposes pre-computed parameters should be moved
1491 : // out of this routine to the calling function.
1492 :
1493 65841 : const GDALGridDataMetricsOptions *const poOptions =
1494 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
1495 : // Pre-compute search ellipse parameters.
1496 65841 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
1497 65841 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
1498 : const double dfSearchRadius =
1499 65841 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
1500 65619 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
1501 :
1502 65619 : GDALGridExtraParameters *psExtraParams =
1503 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
1504 65619 : CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
1505 :
1506 : // Compute coefficients for coordinate system rotation.
1507 65619 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
1508 65619 : const bool bRotated = dfAngle != 0.0;
1509 65619 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
1510 65619 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
1511 :
1512 65619 : double dfMaximumValue = -std::numeric_limits<double>::max();
1513 65619 : double dfMinimumValue = std::numeric_limits<double>::max();
1514 65619 : GUInt32 n = 0;
1515 65619 : if (phQuadTree != nullptr)
1516 : {
1517 : CPLRectObj sAoi;
1518 800 : sAoi.minx = dfXPoint - dfSearchRadius;
1519 800 : sAoi.miny = dfYPoint - dfSearchRadius;
1520 800 : sAoi.maxx = dfXPoint + dfSearchRadius;
1521 800 : sAoi.maxy = dfYPoint + dfSearchRadius;
1522 800 : int nFeatureCount = 0;
1523 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
1524 800 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
1525 800 : if (nFeatureCount != 0)
1526 : {
1527 7458 : for (int k = 0; k < nFeatureCount; k++)
1528 : {
1529 6658 : const int i = papsPoints[k]->i;
1530 6658 : const double dfRX = padfX[i] - dfXPoint;
1531 6658 : const double dfRY = padfY[i] - dfYPoint;
1532 :
1533 6658 : if (dfRadius2Square * dfRX * dfRX +
1534 6658 : dfRadius1Square * dfRY * dfRY <=
1535 : dfR12Square)
1536 : {
1537 6667 : if (dfMinimumValue > padfZ[i])
1538 1945 : dfMinimumValue = padfZ[i];
1539 6667 : if (dfMaximumValue < padfZ[i])
1540 1829 : dfMaximumValue = padfZ[i];
1541 6667 : n++;
1542 : }
1543 : }
1544 : }
1545 800 : CPLFree(papsPoints);
1546 : }
1547 : else
1548 : {
1549 64819 : GUInt32 i = 0;
1550 440936 : while (i < nPoints)
1551 : {
1552 376117 : double dfRX = padfX[i] - dfXPoint;
1553 376117 : double dfRY = padfY[i] - dfYPoint;
1554 :
1555 376117 : if (bRotated)
1556 : {
1557 0 : const double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
1558 0 : const double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
1559 :
1560 0 : dfRX = dfRXRotated;
1561 0 : dfRY = dfRYRotated;
1562 : }
1563 :
1564 : // Is this point located inside the search ellipse?
1565 376117 : if (dfRadius2Square * dfRX * dfRX + dfRadius1Square * dfRY * dfRY <=
1566 : dfR12Square)
1567 : {
1568 366787 : if (dfMinimumValue > padfZ[i])
1569 189478 : dfMinimumValue = padfZ[i];
1570 366787 : if (dfMaximumValue < padfZ[i])
1571 182190 : dfMaximumValue = padfZ[i];
1572 366787 : n++;
1573 : }
1574 :
1575 376117 : i++;
1576 : }
1577 : }
1578 :
1579 65619 : if (n < poOptions->nMinPoints || n == 0)
1580 : {
1581 0 : *pdfValue = poOptions->dfNoDataValue;
1582 : }
1583 : else
1584 : {
1585 65786 : *pdfValue = dfMaximumValue - dfMinimumValue;
1586 : }
1587 :
1588 65619 : return CE_None;
1589 : }
1590 :
1591 : /************************************************************************/
1592 : /* GDALGridDataMetricRangePerQuadrant() */
1593 : /************************************************************************/
1594 :
1595 : /**
1596 : * Data range (data metric), with a per-quadrant search logic.
1597 : */
1598 5 : static CPLErr GDALGridDataMetricRangePerQuadrant(
1599 : const void *poOptionsIn, GUInt32 /* nPoints */, const double *padfX,
1600 : const double *padfY, const double *padfZ, double dfXPoint, double dfYPoint,
1601 : double *pdfValue, void *hExtraParamsIn)
1602 : {
1603 5 : const GDALGridDataMetricsOptions *const poOptions =
1604 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
1605 :
1606 : // Pre-compute search ellipse parameters.
1607 5 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
1608 5 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
1609 : const double dfSearchRadius =
1610 5 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
1611 5 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
1612 :
1613 : // const GUInt32 nMaxPoints = poOptions->nMaxPoints;
1614 5 : const GUInt32 nMinPointsPerQuadrant = poOptions->nMinPointsPerQuadrant;
1615 5 : const GUInt32 nMaxPointsPerQuadrant = poOptions->nMaxPointsPerQuadrant;
1616 :
1617 5 : GDALGridExtraParameters *psExtraParams =
1618 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
1619 5 : const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
1620 5 : CPLAssert(phQuadTree);
1621 :
1622 : CPLRectObj sAoi;
1623 5 : sAoi.minx = dfXPoint - dfSearchRadius;
1624 5 : sAoi.miny = dfYPoint - dfSearchRadius;
1625 5 : sAoi.maxx = dfXPoint + dfSearchRadius;
1626 5 : sAoi.maxy = dfYPoint + dfSearchRadius;
1627 5 : int nFeatureCount = 0;
1628 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
1629 5 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
1630 50 : std::multimap<double, double> oMapDistanceToZValuesPerQuadrant[4];
1631 :
1632 5 : if (nFeatureCount != 0)
1633 : {
1634 26 : for (int k = 0; k < nFeatureCount; k++)
1635 : {
1636 21 : const int i = papsPoints[k]->i;
1637 21 : const double dfRX = padfX[i] - dfXPoint;
1638 21 : const double dfRY = padfY[i] - dfYPoint;
1639 21 : const double dfRXSquare = dfRX * dfRX;
1640 21 : const double dfRYSquare = dfRY * dfRY;
1641 :
1642 21 : if (dfRadius2Square * dfRXSquare + dfRadius1Square * dfRYSquare <=
1643 : dfR12Square)
1644 : {
1645 17 : const int iQuadrant =
1646 17 : ((dfRX >= 0) ? 1 : 0) | (((dfRY >= 0) ? 1 : 0) << 1);
1647 : oMapDistanceToZValuesPerQuadrant[iQuadrant].insert(
1648 17 : std::make_pair(dfRXSquare + dfRYSquare, padfZ[i]));
1649 : }
1650 : }
1651 : }
1652 5 : CPLFree(papsPoints);
1653 :
1654 : std::multimap<double, double>::iterator aoIter[] = {
1655 5 : oMapDistanceToZValuesPerQuadrant[0].begin(),
1656 5 : oMapDistanceToZValuesPerQuadrant[1].begin(),
1657 5 : oMapDistanceToZValuesPerQuadrant[2].begin(),
1658 5 : oMapDistanceToZValuesPerQuadrant[3].begin(),
1659 5 : };
1660 5 : constexpr int ALL_QUADRANT_FLAGS = 1 + 2 + 4 + 8;
1661 :
1662 : // Examine all "neighbors" within the radius (sorted by distance via the
1663 : // multimap), and use the closest n points based on distance until the max
1664 : // is reached.
1665 : // Do that by fetching the nearest point in quadrant 0, then the nearest
1666 : // point in quadrant 1, 2 and 3, and starting again with the next nearest
1667 : // point in quarant 0, etc.
1668 5 : int nQuadrantIterFinishedFlag = 0;
1669 5 : GUInt32 anPerQuadrant[4] = {0};
1670 5 : double dfMaximumValue = -std::numeric_limits<double>::max();
1671 5 : double dfMinimumValue = std::numeric_limits<double>::max();
1672 5 : GUInt32 n = 0;
1673 5 : for (int iQuadrant = 0; /* true */; iQuadrant = (iQuadrant + 1) % 4)
1674 : {
1675 40 : if (aoIter[iQuadrant] ==
1676 90 : oMapDistanceToZValuesPerQuadrant[iQuadrant].end() ||
1677 10 : (nMaxPointsPerQuadrant > 0 &&
1678 10 : anPerQuadrant[iQuadrant] >= nMaxPointsPerQuadrant))
1679 : {
1680 24 : nQuadrantIterFinishedFlag |= 1 << iQuadrant;
1681 24 : if (nQuadrantIterFinishedFlag == ALL_QUADRANT_FLAGS)
1682 5 : break;
1683 19 : continue;
1684 : }
1685 :
1686 16 : const double dfZ = aoIter[iQuadrant]->second;
1687 16 : ++aoIter[iQuadrant];
1688 :
1689 16 : if (dfMinimumValue > dfZ)
1690 7 : dfMinimumValue = dfZ;
1691 16 : if (dfMaximumValue < dfZ)
1692 6 : dfMaximumValue = dfZ;
1693 16 : n++;
1694 16 : anPerQuadrant[iQuadrant]++;
1695 : /*if( nMaxPoints > 0 && n >= nMaxPoints )
1696 : {
1697 : break;
1698 : }*/
1699 35 : }
1700 :
1701 5 : if (nMinPointsPerQuadrant > 0 &&
1702 5 : (anPerQuadrant[0] < nMinPointsPerQuadrant ||
1703 4 : anPerQuadrant[1] < nMinPointsPerQuadrant ||
1704 3 : anPerQuadrant[2] < nMinPointsPerQuadrant ||
1705 3 : anPerQuadrant[3] < nMinPointsPerQuadrant))
1706 : {
1707 2 : *pdfValue = poOptions->dfNoDataValue;
1708 : }
1709 3 : else if (n < poOptions->nMinPoints || n == 0)
1710 : {
1711 1 : *pdfValue = poOptions->dfNoDataValue;
1712 : }
1713 : else
1714 : {
1715 2 : *pdfValue = dfMaximumValue - dfMinimumValue;
1716 : }
1717 :
1718 10 : return CE_None;
1719 : }
1720 :
1721 : /************************************************************************/
1722 : /* GDALGridDataMetricCount() */
1723 : /************************************************************************/
1724 :
1725 : /**
1726 : * Number of data points (data metric).
1727 : *
1728 : * A number of data points found in grid node search ellipse.
1729 : *
1730 : * \f[
1731 : * Z=n
1732 : * \f]
1733 : *
1734 : * where
1735 : * <ul>
1736 : * <li> \f$Z\f$ is a resulting value at the grid node,
1737 : * <li> \f$n\f$ is a total number of points in search ellipse.
1738 : * </ul>
1739 : *
1740 : * @param poOptionsIn Algorithm parameters. This should point to
1741 : * GDALGridDataMetricsOptions object.
1742 : * @param nPoints Number of elements in input arrays.
1743 : * @param padfX Input array of X coordinates.
1744 : * @param padfY Input array of Y coordinates.
1745 : * @param padfZ Input array of Z values.
1746 : * @param dfXPoint X coordinate of the point to compute.
1747 : * @param dfYPoint Y coordinate of the point to compute.
1748 : * @param pdfValue Pointer to variable where the computed grid node value
1749 : * will be returned.
1750 : * @param hExtraParamsIn extra parameters (unused)
1751 : *
1752 : * @return CE_None on success or CE_Failure if something goes wrong.
1753 : */
1754 :
1755 67174 : CPLErr GDALGridDataMetricCount(const void *poOptionsIn, GUInt32 nPoints,
1756 : const double *padfX, const double *padfY,
1757 : CPL_UNUSED const double *padfZ, double dfXPoint,
1758 : double dfYPoint, double *pdfValue,
1759 : void *hExtraParamsIn)
1760 : {
1761 : // TODO: For optimization purposes pre-computed parameters should be moved
1762 : // out of this routine to the calling function.
1763 :
1764 67174 : const GDALGridDataMetricsOptions *const poOptions =
1765 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
1766 :
1767 : // Pre-compute search ellipse parameters.
1768 67174 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
1769 67174 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
1770 : const double dfSearchRadius =
1771 67174 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
1772 67058 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
1773 :
1774 67058 : GDALGridExtraParameters *psExtraParams =
1775 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
1776 67058 : CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
1777 :
1778 : // Compute coefficients for coordinate system rotation.
1779 67058 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
1780 67058 : const bool bRotated = dfAngle != 0.0;
1781 67058 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
1782 67058 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
1783 :
1784 67058 : GUInt32 n = 0;
1785 67058 : if (phQuadTree != nullptr)
1786 : {
1787 : CPLRectObj sAoi;
1788 1999 : sAoi.minx = dfXPoint - dfSearchRadius;
1789 1999 : sAoi.miny = dfYPoint - dfSearchRadius;
1790 1999 : sAoi.maxx = dfXPoint + dfSearchRadius;
1791 1999 : sAoi.maxy = dfYPoint + dfSearchRadius;
1792 1999 : int nFeatureCount = 0;
1793 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
1794 1999 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
1795 1992 : if (nFeatureCount != 0)
1796 : {
1797 79151 : for (int k = 0; k < nFeatureCount; k++)
1798 : {
1799 77159 : const int i = papsPoints[k]->i;
1800 77159 : const double dfRX = padfX[i] - dfXPoint;
1801 77159 : const double dfRY = padfY[i] - dfYPoint;
1802 :
1803 77159 : if (dfRadius2Square * dfRX * dfRX +
1804 77159 : dfRadius1Square * dfRY * dfRY <=
1805 : dfR12Square)
1806 : {
1807 54713 : n++;
1808 : }
1809 : }
1810 : }
1811 1992 : CPLFree(papsPoints);
1812 : }
1813 : else
1814 : {
1815 65059 : GUInt32 i = 0;
1816 385678 : while (i < nPoints)
1817 : {
1818 320619 : double dfRX = padfX[i] - dfXPoint;
1819 320619 : double dfRY = padfY[i] - dfYPoint;
1820 :
1821 320619 : if (bRotated)
1822 : {
1823 0 : const double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
1824 0 : const double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
1825 :
1826 0 : dfRX = dfRXRotated;
1827 0 : dfRY = dfRYRotated;
1828 : }
1829 :
1830 : // Is this point located inside the search ellipse?
1831 320619 : if (dfRadius2Square * dfRX * dfRX + dfRadius1Square * dfRY * dfRY <=
1832 : dfR12Square)
1833 : {
1834 70199 : n++;
1835 : }
1836 :
1837 320619 : i++;
1838 : }
1839 : }
1840 :
1841 67059 : if (n < poOptions->nMinPoints)
1842 : {
1843 0 : *pdfValue = poOptions->dfNoDataValue;
1844 : }
1845 : else
1846 : {
1847 67059 : *pdfValue = static_cast<double>(n);
1848 : }
1849 :
1850 67059 : return CE_None;
1851 : }
1852 :
1853 : /************************************************************************/
1854 : /* GDALGridDataMetricCountPerQuadrant() */
1855 : /************************************************************************/
1856 :
1857 : /**
1858 : * Number of data points (data metric), with a per-quadrant search logic.
1859 : */
1860 5 : static CPLErr GDALGridDataMetricCountPerQuadrant(
1861 : const void *poOptionsIn, GUInt32 /* nPoints */, const double *padfX,
1862 : const double *padfY, const double *padfZ, double dfXPoint, double dfYPoint,
1863 : double *pdfValue, void *hExtraParamsIn)
1864 : {
1865 5 : const GDALGridDataMetricsOptions *const poOptions =
1866 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
1867 :
1868 : // Pre-compute search ellipse parameters.
1869 5 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
1870 5 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
1871 : const double dfSearchRadius =
1872 5 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
1873 5 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
1874 :
1875 : // const GUInt32 nMaxPoints = poOptions->nMaxPoints;
1876 5 : const GUInt32 nMinPointsPerQuadrant = poOptions->nMinPointsPerQuadrant;
1877 5 : const GUInt32 nMaxPointsPerQuadrant = poOptions->nMaxPointsPerQuadrant;
1878 :
1879 5 : GDALGridExtraParameters *psExtraParams =
1880 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
1881 5 : const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
1882 5 : CPLAssert(phQuadTree);
1883 :
1884 : CPLRectObj sAoi;
1885 5 : sAoi.minx = dfXPoint - dfSearchRadius;
1886 5 : sAoi.miny = dfYPoint - dfSearchRadius;
1887 5 : sAoi.maxx = dfXPoint + dfSearchRadius;
1888 5 : sAoi.maxy = dfYPoint + dfSearchRadius;
1889 5 : int nFeatureCount = 0;
1890 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
1891 5 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
1892 50 : std::multimap<double, double> oMapDistanceToZValuesPerQuadrant[4];
1893 :
1894 5 : if (nFeatureCount != 0)
1895 : {
1896 26 : for (int k = 0; k < nFeatureCount; k++)
1897 : {
1898 21 : const int i = papsPoints[k]->i;
1899 21 : const double dfRX = padfX[i] - dfXPoint;
1900 21 : const double dfRY = padfY[i] - dfYPoint;
1901 21 : const double dfRXSquare = dfRX * dfRX;
1902 21 : const double dfRYSquare = dfRY * dfRY;
1903 :
1904 21 : if (dfRadius2Square * dfRXSquare + dfRadius1Square * dfRYSquare <=
1905 : dfR12Square)
1906 : {
1907 17 : const int iQuadrant =
1908 17 : ((dfRX >= 0) ? 1 : 0) | (((dfRY >= 0) ? 1 : 0) << 1);
1909 : oMapDistanceToZValuesPerQuadrant[iQuadrant].insert(
1910 17 : std::make_pair(dfRXSquare + dfRYSquare, padfZ[i]));
1911 : }
1912 : }
1913 : }
1914 5 : CPLFree(papsPoints);
1915 :
1916 : std::multimap<double, double>::iterator aoIter[] = {
1917 5 : oMapDistanceToZValuesPerQuadrant[0].begin(),
1918 5 : oMapDistanceToZValuesPerQuadrant[1].begin(),
1919 5 : oMapDistanceToZValuesPerQuadrant[2].begin(),
1920 5 : oMapDistanceToZValuesPerQuadrant[3].begin(),
1921 5 : };
1922 5 : constexpr int ALL_QUADRANT_FLAGS = 1 + 2 + 4 + 8;
1923 :
1924 : // Examine all "neighbors" within the radius (sorted by distance via the
1925 : // multimap), and use the closest n points based on distance until the max
1926 : // is reached.
1927 : // Do that by fetching the nearest point in quadrant 0, then the nearest
1928 : // point in quadrant 1, 2 and 3, and starting again with the next nearest
1929 : // point in quarant 0, etc.
1930 5 : int nQuadrantIterFinishedFlag = 0;
1931 5 : GUInt32 anPerQuadrant[4] = {0};
1932 5 : GUInt32 n = 0;
1933 40 : for (int iQuadrant = 0; /* true */; iQuadrant = (iQuadrant + 1) % 4)
1934 : {
1935 40 : if (aoIter[iQuadrant] ==
1936 90 : oMapDistanceToZValuesPerQuadrant[iQuadrant].end() ||
1937 10 : (nMaxPointsPerQuadrant > 0 &&
1938 10 : anPerQuadrant[iQuadrant] >= nMaxPointsPerQuadrant))
1939 : {
1940 24 : nQuadrantIterFinishedFlag |= 1 << iQuadrant;
1941 24 : if (nQuadrantIterFinishedFlag == ALL_QUADRANT_FLAGS)
1942 5 : break;
1943 19 : continue;
1944 : }
1945 :
1946 16 : ++aoIter[iQuadrant];
1947 :
1948 16 : n++;
1949 16 : anPerQuadrant[iQuadrant]++;
1950 : /*if( nMaxPoints > 0 && n >= nMaxPoints )
1951 : {
1952 : break;
1953 : }*/
1954 : }
1955 :
1956 5 : if (nMinPointsPerQuadrant > 0 &&
1957 5 : (anPerQuadrant[0] < nMinPointsPerQuadrant ||
1958 4 : anPerQuadrant[1] < nMinPointsPerQuadrant ||
1959 3 : anPerQuadrant[2] < nMinPointsPerQuadrant ||
1960 3 : anPerQuadrant[3] < nMinPointsPerQuadrant))
1961 : {
1962 2 : *pdfValue = poOptions->dfNoDataValue;
1963 : }
1964 3 : else if (n < poOptions->nMinPoints)
1965 : {
1966 1 : *pdfValue = poOptions->dfNoDataValue;
1967 : }
1968 : else
1969 : {
1970 2 : *pdfValue = static_cast<double>(n);
1971 : }
1972 :
1973 10 : return CE_None;
1974 : }
1975 :
1976 : /************************************************************************/
1977 : /* GDALGridDataMetricAverageDistance() */
1978 : /************************************************************************/
1979 :
1980 : /**
1981 : * Average distance (data metric).
1982 : *
1983 : * An average distance between the grid node (center of the search ellipse)
1984 : * and all of the data points found in grid node search ellipse. If there are
1985 : * no points found, the specified NODATA value will be returned.
1986 : *
1987 : * \f[
1988 : * Z=\frac{\sum_{i = 1}^n r_i}{n}
1989 : * \f]
1990 : *
1991 : * where
1992 : * <ul>
1993 : * <li> \f$Z\f$ is a resulting value at the grid node,
1994 : * <li> \f$r_i\f$ is an Euclidean distance from the grid node
1995 : * to point \f$i\f$,
1996 : * <li> \f$n\f$ is a total number of points in search ellipse.
1997 : * </ul>
1998 : *
1999 : * @param poOptionsIn Algorithm parameters. This should point to
2000 : * GDALGridDataMetricsOptions object.
2001 : * @param nPoints Number of elements in input arrays.
2002 : * @param padfX Input array of X coordinates.
2003 : * @param padfY Input array of Y coordinates.
2004 : * @param padfZ Input array of Z values (unused)
2005 : * @param dfXPoint X coordinate of the point to compute.
2006 : * @param dfYPoint Y coordinate of the point to compute.
2007 : * @param pdfValue Pointer to variable where the computed grid node value
2008 : * will be returned.
2009 : * @param hExtraParamsIn extra parameters (unused)
2010 : *
2011 : * @return CE_None on success or CE_Failure if something goes wrong.
2012 : */
2013 :
2014 66545 : CPLErr GDALGridDataMetricAverageDistance(const void *poOptionsIn,
2015 : GUInt32 nPoints, const double *padfX,
2016 : const double *padfY,
2017 : CPL_UNUSED const double *padfZ,
2018 : double dfXPoint, double dfYPoint,
2019 : double *pdfValue, void *hExtraParamsIn)
2020 : {
2021 : // TODO: For optimization purposes pre-computed parameters should be moved
2022 : // out of this routine to the calling function.
2023 :
2024 66545 : const GDALGridDataMetricsOptions *const poOptions =
2025 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
2026 :
2027 : // Pre-compute search ellipse parameters.
2028 66545 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
2029 66545 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
2030 : const double dfSearchRadius =
2031 66545 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
2032 66617 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
2033 :
2034 66617 : GDALGridExtraParameters *psExtraParams =
2035 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
2036 66617 : CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
2037 :
2038 : // Compute coefficients for coordinate system rotation.
2039 66617 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
2040 66617 : const bool bRotated = dfAngle != 0.0;
2041 66617 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
2042 66617 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
2043 :
2044 66617 : double dfAccumulator = 0.0;
2045 66617 : GUInt32 n = 0;
2046 66617 : if (phQuadTree != nullptr)
2047 : {
2048 : CPLRectObj sAoi;
2049 798 : sAoi.minx = dfXPoint - dfSearchRadius;
2050 798 : sAoi.miny = dfYPoint - dfSearchRadius;
2051 798 : sAoi.maxx = dfXPoint + dfSearchRadius;
2052 798 : sAoi.maxy = dfYPoint + dfSearchRadius;
2053 798 : int nFeatureCount = 0;
2054 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
2055 798 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
2056 800 : if (nFeatureCount != 0)
2057 : {
2058 17918 : for (int k = 0; k < nFeatureCount; k++)
2059 : {
2060 17118 : const int i = papsPoints[k]->i;
2061 17118 : const double dfRX = padfX[i] - dfXPoint;
2062 17118 : const double dfRY = padfY[i] - dfYPoint;
2063 :
2064 17118 : if (dfRadius2Square * dfRX * dfRX +
2065 17118 : dfRadius1Square * dfRY * dfRY <=
2066 : dfR12Square)
2067 : {
2068 14614 : dfAccumulator += sqrt(dfRX * dfRX + dfRY * dfRY);
2069 14614 : n++;
2070 : }
2071 : }
2072 : }
2073 800 : CPLFree(papsPoints);
2074 : }
2075 : else
2076 : {
2077 65819 : GUInt32 i = 0;
2078 :
2079 376962 : while (i < nPoints)
2080 : {
2081 311143 : double dfRX = padfX[i] - dfXPoint;
2082 311143 : double dfRY = padfY[i] - dfYPoint;
2083 :
2084 311143 : if (bRotated)
2085 : {
2086 0 : const double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
2087 0 : const double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
2088 :
2089 0 : dfRX = dfRXRotated;
2090 0 : dfRY = dfRYRotated;
2091 : }
2092 :
2093 : // Is this point located inside the search ellipse?
2094 311143 : if (dfRadius2Square * dfRX * dfRX + dfRadius1Square * dfRY * dfRY <=
2095 : dfR12Square)
2096 : {
2097 426786 : dfAccumulator += sqrt(dfRX * dfRX + dfRY * dfRY);
2098 426786 : n++;
2099 : }
2100 :
2101 311143 : i++;
2102 : }
2103 : }
2104 :
2105 66617 : if (n < poOptions->nMinPoints || n == 0)
2106 : {
2107 2 : *pdfValue = poOptions->dfNoDataValue;
2108 : }
2109 : else
2110 : {
2111 66615 : *pdfValue = dfAccumulator / n;
2112 : }
2113 :
2114 66617 : return CE_None;
2115 : }
2116 :
2117 : /************************************************************************/
2118 : /* GDALGridDataMetricAverageDistancePerQuadrant() */
2119 : /************************************************************************/
2120 :
2121 : /**
2122 : * Average distance (data metric), with a per-quadrant search logic.
2123 : */
2124 5 : static CPLErr GDALGridDataMetricAverageDistancePerQuadrant(
2125 : const void *poOptionsIn, GUInt32 /* nPoints */, const double *padfX,
2126 : const double *padfY, const double *padfZ, double dfXPoint, double dfYPoint,
2127 : double *pdfValue, void *hExtraParamsIn)
2128 : {
2129 5 : const GDALGridDataMetricsOptions *const poOptions =
2130 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
2131 :
2132 : // Pre-compute search ellipse parameters.
2133 5 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
2134 5 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
2135 : const double dfSearchRadius =
2136 5 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
2137 5 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
2138 :
2139 : // const GUInt32 nMaxPoints = poOptions->nMaxPoints;
2140 5 : const GUInt32 nMinPointsPerQuadrant = poOptions->nMinPointsPerQuadrant;
2141 5 : const GUInt32 nMaxPointsPerQuadrant = poOptions->nMaxPointsPerQuadrant;
2142 :
2143 5 : GDALGridExtraParameters *psExtraParams =
2144 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
2145 5 : const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
2146 5 : CPLAssert(phQuadTree);
2147 :
2148 : CPLRectObj sAoi;
2149 5 : sAoi.minx = dfXPoint - dfSearchRadius;
2150 5 : sAoi.miny = dfYPoint - dfSearchRadius;
2151 5 : sAoi.maxx = dfXPoint + dfSearchRadius;
2152 5 : sAoi.maxy = dfYPoint + dfSearchRadius;
2153 5 : int nFeatureCount = 0;
2154 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
2155 5 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
2156 50 : std::multimap<double, double> oMapDistanceToZValuesPerQuadrant[4];
2157 :
2158 5 : if (nFeatureCount != 0)
2159 : {
2160 26 : for (int k = 0; k < nFeatureCount; k++)
2161 : {
2162 21 : const int i = papsPoints[k]->i;
2163 21 : const double dfRX = padfX[i] - dfXPoint;
2164 21 : const double dfRY = padfY[i] - dfYPoint;
2165 21 : const double dfRXSquare = dfRX * dfRX;
2166 21 : const double dfRYSquare = dfRY * dfRY;
2167 :
2168 21 : if (dfRadius2Square * dfRXSquare + dfRadius1Square * dfRYSquare <=
2169 : dfR12Square)
2170 : {
2171 17 : const int iQuadrant =
2172 17 : ((dfRX >= 0) ? 1 : 0) | (((dfRY >= 0) ? 1 : 0) << 1);
2173 : oMapDistanceToZValuesPerQuadrant[iQuadrant].insert(
2174 17 : std::make_pair(dfRXSquare + dfRYSquare, padfZ[i]));
2175 : }
2176 : }
2177 : }
2178 5 : CPLFree(papsPoints);
2179 :
2180 : std::multimap<double, double>::iterator aoIter[] = {
2181 5 : oMapDistanceToZValuesPerQuadrant[0].begin(),
2182 5 : oMapDistanceToZValuesPerQuadrant[1].begin(),
2183 5 : oMapDistanceToZValuesPerQuadrant[2].begin(),
2184 5 : oMapDistanceToZValuesPerQuadrant[3].begin(),
2185 5 : };
2186 5 : constexpr int ALL_QUADRANT_FLAGS = 1 + 2 + 4 + 8;
2187 :
2188 : // Examine all "neighbors" within the radius (sorted by distance via the
2189 : // multimap), and use the closest n points based on distance until the max
2190 : // is reached.
2191 : // Do that by fetching the nearest point in quadrant 0, then the nearest
2192 : // point in quadrant 1, 2 and 3, and starting again with the next nearest
2193 : // point in quarant 0, etc.
2194 5 : int nQuadrantIterFinishedFlag = 0;
2195 5 : GUInt32 anPerQuadrant[4] = {0};
2196 5 : GUInt32 n = 0;
2197 5 : double dfAccumulator = 0;
2198 40 : for (int iQuadrant = 0; /* true */; iQuadrant = (iQuadrant + 1) % 4)
2199 : {
2200 40 : if (aoIter[iQuadrant] ==
2201 90 : oMapDistanceToZValuesPerQuadrant[iQuadrant].end() ||
2202 10 : (nMaxPointsPerQuadrant > 0 &&
2203 10 : anPerQuadrant[iQuadrant] >= nMaxPointsPerQuadrant))
2204 : {
2205 24 : nQuadrantIterFinishedFlag |= 1 << iQuadrant;
2206 24 : if (nQuadrantIterFinishedFlag == ALL_QUADRANT_FLAGS)
2207 5 : break;
2208 19 : continue;
2209 : }
2210 :
2211 16 : dfAccumulator += sqrt(aoIter[iQuadrant]->first);
2212 16 : ++aoIter[iQuadrant];
2213 :
2214 16 : n++;
2215 16 : anPerQuadrant[iQuadrant]++;
2216 : /*if( nMaxPoints > 0 && n >= nMaxPoints )
2217 : {
2218 : break;
2219 : }*/
2220 : }
2221 :
2222 5 : if (nMinPointsPerQuadrant > 0 &&
2223 5 : (anPerQuadrant[0] < nMinPointsPerQuadrant ||
2224 4 : anPerQuadrant[1] < nMinPointsPerQuadrant ||
2225 3 : anPerQuadrant[2] < nMinPointsPerQuadrant ||
2226 3 : anPerQuadrant[3] < nMinPointsPerQuadrant))
2227 : {
2228 2 : *pdfValue = poOptions->dfNoDataValue;
2229 : }
2230 3 : else if (n < poOptions->nMinPoints || n == 0)
2231 : {
2232 1 : *pdfValue = poOptions->dfNoDataValue;
2233 : }
2234 : else
2235 : {
2236 2 : *pdfValue = dfAccumulator / n;
2237 : }
2238 :
2239 10 : return CE_None;
2240 : }
2241 :
2242 : /************************************************************************/
2243 : /* GDALGridDataMetricAverageDistancePts() */
2244 : /************************************************************************/
2245 :
2246 : /**
2247 : * Average distance between points (data metric).
2248 : *
2249 : * An average distance between the data points found in grid node search
2250 : * ellipse. The distance between each pair of points within ellipse is
2251 : * calculated and average of all distances is set as a grid node value. If
2252 : * there are no points found, the specified NODATA value will be returned.
2253 :
2254 : *
2255 : * \f[
2256 : * Z=\frac{\sum_{i = 1}^{n-1}\sum_{j=i+1}^{n}
2257 : r_{ij}}{\left(n-1\right)\,n-\frac{n+{\left(n-1\right)}^{2}-1}{2}}
2258 : * \f]
2259 : *
2260 : * where
2261 : * <ul>
2262 : * <li> \f$Z\f$ is a resulting value at the grid node,
2263 : * <li> \f$r_{ij}\f$ is an Euclidean distance between points
2264 : * \f$i\f$ and \f$j\f$,
2265 : * <li> \f$n\f$ is a total number of points in search ellipse.
2266 : * </ul>
2267 : *
2268 : * @param poOptionsIn Algorithm parameters. This should point to
2269 : * GDALGridDataMetricsOptions object.
2270 : * @param nPoints Number of elements in input arrays.
2271 : * @param padfX Input array of X coordinates.
2272 : * @param padfY Input array of Y coordinates.
2273 : * @param padfZ Input array of Z values (unused)
2274 : * @param dfXPoint X coordinate of the point to compute.
2275 : * @param dfYPoint Y coordinate of the point to compute.
2276 : * @param pdfValue Pointer to variable where the computed grid node value
2277 : * will be returned.
2278 : * @param hExtraParamsIn extra parameters (unused)
2279 : *
2280 : * @return CE_None on success or CE_Failure if something goes wrong.
2281 : */
2282 :
2283 66598 : CPLErr GDALGridDataMetricAverageDistancePts(
2284 : const void *poOptionsIn, GUInt32 nPoints, const double *padfX,
2285 : const double *padfY, CPL_UNUSED const double *padfZ, double dfXPoint,
2286 : double dfYPoint, double *pdfValue, void *hExtraParamsIn)
2287 : {
2288 : // TODO: For optimization purposes pre-computed parameters should be moved
2289 : // out of this routine to the calling function.
2290 :
2291 66598 : const GDALGridDataMetricsOptions *const poOptions =
2292 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
2293 : // Pre-compute search ellipse parameters.
2294 66598 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
2295 66598 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
2296 : const double dfSearchRadius =
2297 66598 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
2298 66517 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
2299 :
2300 66517 : GDALGridExtraParameters *psExtraParams =
2301 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
2302 66517 : CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
2303 :
2304 : // Compute coefficients for coordinate system rotation.
2305 66517 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
2306 66517 : const bool bRotated = dfAngle != 0.0;
2307 66517 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
2308 66517 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
2309 :
2310 66517 : double dfAccumulator = 0.0;
2311 66517 : GUInt32 n = 0;
2312 66517 : if (phQuadTree != nullptr)
2313 : {
2314 : CPLRectObj sAoi;
2315 800 : sAoi.minx = dfXPoint - dfSearchRadius;
2316 800 : sAoi.miny = dfYPoint - dfSearchRadius;
2317 800 : sAoi.maxx = dfXPoint + dfSearchRadius;
2318 800 : sAoi.maxy = dfYPoint + dfSearchRadius;
2319 800 : int nFeatureCount = 0;
2320 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
2321 800 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
2322 800 : if (nFeatureCount != 0)
2323 : {
2324 17444 : for (int k = 0; k < nFeatureCount - 1; k++)
2325 : {
2326 16644 : const int i = papsPoints[k]->i;
2327 16644 : const double dfRX1 = padfX[i] - dfXPoint;
2328 16644 : const double dfRY1 = padfY[i] - dfYPoint;
2329 :
2330 16644 : if (dfRadius2Square * dfRX1 * dfRX1 +
2331 16644 : dfRadius1Square * dfRY1 * dfRY1 <=
2332 : dfR12Square)
2333 : {
2334 163975 : for (int j = k; j < nFeatureCount; j++)
2335 : // Search all the remaining points within the ellipse and
2336 : // compute distances between them and the first point.
2337 : {
2338 149496 : const int ji = papsPoints[j]->i;
2339 149496 : double dfRX2 = padfX[ji] - dfXPoint;
2340 149496 : double dfRY2 = padfY[ji] - dfYPoint;
2341 :
2342 149496 : if (dfRadius2Square * dfRX2 * dfRX2 +
2343 149496 : dfRadius1Square * dfRY2 * dfRY2 <=
2344 : dfR12Square)
2345 : {
2346 131826 : const double dfRX = padfX[ji] - padfX[i];
2347 131826 : const double dfRY = padfY[ji] - padfY[i];
2348 :
2349 131826 : dfAccumulator += sqrt(dfRX * dfRX + dfRY * dfRY);
2350 131826 : n++;
2351 : }
2352 : }
2353 : }
2354 : }
2355 : }
2356 800 : CPLFree(papsPoints);
2357 : }
2358 : else
2359 : {
2360 65717 : GUInt32 i = 0;
2361 446358 : while (i < nPoints - 1)
2362 : {
2363 380641 : double dfRX1 = padfX[i] - dfXPoint;
2364 380641 : double dfRY1 = padfY[i] - dfYPoint;
2365 :
2366 380641 : if (bRotated)
2367 : {
2368 148652 : const double dfRXRotated = dfRX1 * dfCoeff1 + dfRY1 * dfCoeff2;
2369 148652 : const double dfRYRotated = dfRY1 * dfCoeff1 - dfRX1 * dfCoeff2;
2370 :
2371 148652 : dfRX1 = dfRXRotated;
2372 148652 : dfRY1 = dfRYRotated;
2373 : }
2374 :
2375 : // Is this point located inside the search ellipse?
2376 380641 : if (dfRadius2Square * dfRX1 * dfRX1 +
2377 380641 : dfRadius1Square * dfRY1 * dfRY1 <=
2378 : dfR12Square)
2379 : {
2380 : // Search all the remaining points within the ellipse and
2381 : // compute distances between them and the first point.
2382 1107830 : for (GUInt32 j = i + 1; j < nPoints; j++)
2383 : {
2384 859434 : double dfRX2 = padfX[j] - dfXPoint;
2385 859434 : double dfRY2 = padfY[j] - dfYPoint;
2386 :
2387 859434 : if (bRotated)
2388 : {
2389 391456 : const double dfRXRotated =
2390 391456 : dfRX2 * dfCoeff1 + dfRY2 * dfCoeff2;
2391 391456 : const double dfRYRotated =
2392 391456 : dfRY2 * dfCoeff1 - dfRX2 * dfCoeff2;
2393 :
2394 391456 : dfRX2 = dfRXRotated;
2395 391456 : dfRY2 = dfRYRotated;
2396 : }
2397 :
2398 859434 : if (dfRadius2Square * dfRX2 * dfRX2 +
2399 859434 : dfRadius1Square * dfRY2 * dfRY2 <=
2400 : dfR12Square)
2401 : {
2402 555083 : const double dfRX = padfX[j] - padfX[i];
2403 555083 : const double dfRY = padfY[j] - padfY[i];
2404 :
2405 555083 : dfAccumulator += sqrt(dfRX * dfRX + dfRY * dfRY);
2406 555083 : n++;
2407 : }
2408 : }
2409 : }
2410 :
2411 380641 : i++;
2412 : }
2413 : }
2414 :
2415 : // Search for the first point within the search ellipse.
2416 66516 : if (n < poOptions->nMinPoints || n == 0)
2417 : {
2418 0 : *pdfValue = poOptions->dfNoDataValue;
2419 : }
2420 : else
2421 : {
2422 66593 : *pdfValue = dfAccumulator / n;
2423 : }
2424 :
2425 66516 : return CE_None;
2426 : }
2427 :
2428 : /************************************************************************/
2429 : /* GDALGridLinear() */
2430 : /************************************************************************/
2431 :
2432 : /**
2433 : * Linear interpolation
2434 : *
2435 : * The Linear method performs linear interpolation by finding in which triangle
2436 : * of a Delaunay triangulation the point is, and by doing interpolation from
2437 : * its barycentric coordinates within the triangle.
2438 : * If the point is not in any triangle, depending on the radius, the
2439 : * algorithm will use the value of the nearest point (radius != 0),
2440 : * or the nodata value (radius == 0)
2441 : *
2442 : * @param poOptionsIn Algorithm parameters. This should point to
2443 : * GDALGridLinearOptions object.
2444 : * @param nPoints Number of elements in input arrays.
2445 : * @param padfX Input array of X coordinates.
2446 : * @param padfY Input array of Y coordinates.
2447 : * @param padfZ Input array of Z values.
2448 : * @param dfXPoint X coordinate of the point to compute.
2449 : * @param dfYPoint Y coordinate of the point to compute.
2450 : * @param pdfValue Pointer to variable where the computed grid node value
2451 : * will be returned.
2452 : * @param hExtraParams extra parameters
2453 : *
2454 : * @return CE_None on success or CE_Failure if something goes wrong.
2455 : *
2456 : */
2457 :
2458 342151 : CPLErr GDALGridLinear(const void *poOptionsIn, GUInt32 nPoints,
2459 : const double *padfX, const double *padfY,
2460 : const double *padfZ, double dfXPoint, double dfYPoint,
2461 : double *pdfValue, void *hExtraParams)
2462 : {
2463 342151 : GDALGridExtraParameters *psExtraParams =
2464 : static_cast<GDALGridExtraParameters *>(hExtraParams);
2465 342151 : GDALTriangulation *psTriangulation = psExtraParams->psTriangulation;
2466 :
2467 342151 : int nOutputFacetIdx = -1;
2468 342151 : const bool bRet = CPL_TO_BOOL(GDALTriangulationFindFacetDirected(
2469 : psTriangulation, psExtraParams->nInitialFacetIdx, dfXPoint, dfYPoint,
2470 : &nOutputFacetIdx));
2471 :
2472 343758 : if (bRet)
2473 : {
2474 82197 : CPLAssert(nOutputFacetIdx >= 0);
2475 : // Reuse output facet idx as next initial index since we proceed line by
2476 : // line.
2477 82197 : psExtraParams->nInitialFacetIdx = nOutputFacetIdx;
2478 :
2479 82197 : double lambda1 = 0.0;
2480 82197 : double lambda2 = 0.0;
2481 82197 : double lambda3 = 0.0;
2482 82197 : GDALTriangulationComputeBarycentricCoordinates(
2483 : psTriangulation, nOutputFacetIdx, dfXPoint, dfYPoint, &lambda1,
2484 : &lambda2, &lambda3);
2485 80777 : const int i1 =
2486 80777 : psTriangulation->pasFacets[nOutputFacetIdx].anVertexIdx[0];
2487 80777 : const int i2 =
2488 80777 : psTriangulation->pasFacets[nOutputFacetIdx].anVertexIdx[1];
2489 80777 : const int i3 =
2490 80777 : psTriangulation->pasFacets[nOutputFacetIdx].anVertexIdx[2];
2491 80777 : *pdfValue =
2492 80777 : lambda1 * padfZ[i1] + lambda2 * padfZ[i2] + lambda3 * padfZ[i3];
2493 : }
2494 : else
2495 : {
2496 261561 : if (nOutputFacetIdx >= 0)
2497 : {
2498 : // Also reuse this failed output facet, when valid, as seed for
2499 : // next search.
2500 251996 : psExtraParams->nInitialFacetIdx = nOutputFacetIdx;
2501 : }
2502 :
2503 261561 : const GDALGridLinearOptions *const poOptions =
2504 : static_cast<const GDALGridLinearOptions *>(poOptionsIn);
2505 261561 : const double dfRadius = poOptions->dfRadius;
2506 261561 : if (dfRadius == 0.0)
2507 : {
2508 124510 : *pdfValue = poOptions->dfNoDataValue;
2509 : }
2510 : else
2511 : {
2512 : GDALGridNearestNeighborOptions sNeighbourOptions;
2513 137051 : sNeighbourOptions.nSizeOfStructure = sizeof(sNeighbourOptions);
2514 137019 : sNeighbourOptions.dfRadius1 =
2515 126990 : dfRadius < 0.0 || dfRadius >= std::numeric_limits<double>::max()
2516 264009 : ? 0.0
2517 : : dfRadius;
2518 136963 : sNeighbourOptions.dfRadius2 =
2519 127175 : dfRadius < 0.0 || dfRadius >= std::numeric_limits<double>::max()
2520 264138 : ? 0.0
2521 : : dfRadius;
2522 136963 : sNeighbourOptions.dfAngle = 0.0;
2523 136963 : sNeighbourOptions.dfNoDataValue = poOptions->dfNoDataValue;
2524 136963 : return GDALGridNearestNeighbor(&sNeighbourOptions, nPoints, padfX,
2525 : padfY, padfZ, dfXPoint, dfYPoint,
2526 136586 : pdfValue, hExtraParams);
2527 : }
2528 : }
2529 :
2530 205287 : return CE_None;
2531 : }
2532 :
2533 : /************************************************************************/
2534 : /* GDALGridJob */
2535 : /************************************************************************/
2536 :
2537 : typedef struct _GDALGridJob GDALGridJob;
2538 :
2539 : struct _GDALGridJob
2540 : {
2541 : GUInt32 nYStart;
2542 :
2543 : GByte *pabyData;
2544 : GUInt32 nYStep;
2545 : GUInt32 nXSize;
2546 : GUInt32 nYSize;
2547 : double dfXMin;
2548 : double dfYMin;
2549 : double dfDeltaX;
2550 : double dfDeltaY;
2551 : GUInt32 nPoints;
2552 : const double *padfX;
2553 : const double *padfY;
2554 : const double *padfZ;
2555 : const void *poOptions;
2556 : GDALGridFunction pfnGDALGridMethod;
2557 : GDALGridExtraParameters *psExtraParameters;
2558 : int (*pfnProgress)(GDALGridJob *psJob);
2559 : GDALDataType eType;
2560 :
2561 : int *pnCounter;
2562 : int nCounterSingleThreaded;
2563 : volatile int *pbStop;
2564 : CPLCond *hCond;
2565 : CPLMutex *hCondMutex;
2566 :
2567 : GDALProgressFunc pfnRealProgress;
2568 : void *pRealProgressArg;
2569 : };
2570 :
2571 : /************************************************************************/
2572 : /* GDALGridProgressMultiThread() */
2573 : /************************************************************************/
2574 :
2575 : // Return TRUE if the computation must be interrupted.
2576 20816 : static int GDALGridProgressMultiThread(GDALGridJob *psJob)
2577 : {
2578 20816 : CPLAcquireMutex(psJob->hCondMutex, 1.0);
2579 20852 : ++(*psJob->pnCounter);
2580 20852 : CPLCondSignal(psJob->hCond);
2581 20852 : const int bStop = *psJob->pbStop;
2582 20852 : CPLReleaseMutex(psJob->hCondMutex);
2583 :
2584 20852 : return bStop;
2585 : }
2586 :
2587 : /************************************************************************/
2588 : /* GDALGridProgressMonoThread() */
2589 : /************************************************************************/
2590 :
2591 : // Return TRUE if the computation must be interrupted.
2592 40 : static int GDALGridProgressMonoThread(GDALGridJob *psJob)
2593 : {
2594 40 : const int nCounter = ++(psJob->nCounterSingleThreaded);
2595 40 : if (!psJob->pfnRealProgress(nCounter / static_cast<double>(psJob->nYSize),
2596 : "", psJob->pRealProgressArg))
2597 : {
2598 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
2599 0 : *psJob->pbStop = TRUE;
2600 0 : return TRUE;
2601 : }
2602 40 : return FALSE;
2603 : }
2604 :
2605 : /************************************************************************/
2606 : /* GDALGridJobProcess() */
2607 : /************************************************************************/
2608 :
2609 726 : static void GDALGridJobProcess(void *user_data)
2610 : {
2611 726 : GDALGridJob *const psJob = static_cast<GDALGridJob *>(user_data);
2612 726 : int (*pfnProgress)(GDALGridJob * psJob) = psJob->pfnProgress;
2613 726 : const GUInt32 nXSize = psJob->nXSize;
2614 :
2615 : /* -------------------------------------------------------------------- */
2616 : /* Allocate a buffer of scanline size, fill it with gridded values */
2617 : /* and use GDALCopyWords() to copy values into output data array with */
2618 : /* appropriate data type conversion. */
2619 : /* -------------------------------------------------------------------- */
2620 : double *padfValues =
2621 726 : static_cast<double *>(VSI_MALLOC2_VERBOSE(sizeof(double), nXSize));
2622 717 : if (padfValues == nullptr)
2623 : {
2624 0 : *(psJob->pbStop) = TRUE;
2625 0 : if (pfnProgress != nullptr)
2626 0 : pfnProgress(psJob); // To notify the main thread.
2627 0 : return;
2628 : }
2629 :
2630 717 : const GUInt32 nYStart = psJob->nYStart;
2631 717 : const GUInt32 nYStep = psJob->nYStep;
2632 717 : GByte *pabyData = psJob->pabyData;
2633 :
2634 717 : const GUInt32 nYSize = psJob->nYSize;
2635 717 : const double dfXMin = psJob->dfXMin;
2636 717 : const double dfYMin = psJob->dfYMin;
2637 717 : const double dfDeltaX = psJob->dfDeltaX;
2638 717 : const double dfDeltaY = psJob->dfDeltaY;
2639 717 : const GUInt32 nPoints = psJob->nPoints;
2640 717 : const double *padfX = psJob->padfX;
2641 717 : const double *padfY = psJob->padfY;
2642 717 : const double *padfZ = psJob->padfZ;
2643 717 : const void *poOptions = psJob->poOptions;
2644 717 : GDALGridFunction pfnGDALGridMethod = psJob->pfnGDALGridMethod;
2645 : // Have a local copy of sExtraParameters since we want to modify
2646 : // nInitialFacetIdx.
2647 717 : GDALGridExtraParameters sExtraParameters = *psJob->psExtraParameters;
2648 717 : const GDALDataType eType = psJob->eType;
2649 :
2650 717 : const int nDataTypeSize = GDALGetDataTypeSizeBytes(eType);
2651 716 : const int nLineSpace = nXSize * nDataTypeSize;
2652 :
2653 21594 : for (GUInt32 nYPoint = nYStart; nYPoint < nYSize; nYPoint += nYStep)
2654 : {
2655 20880 : const double dfYPoint = dfYMin + (nYPoint + 0.5) * dfDeltaY;
2656 :
2657 4746500 : for (GUInt32 nXPoint = 0; nXPoint < nXSize; nXPoint++)
2658 : {
2659 4730030 : const double dfXPoint = dfXMin + (nXPoint + 0.5) * dfDeltaX;
2660 :
2661 9455650 : if ((*pfnGDALGridMethod)(poOptions, nPoints, padfX, padfY, padfZ,
2662 4730030 : dfXPoint, dfYPoint, padfValues + nXPoint,
2663 4725620 : &sExtraParameters) != CE_None)
2664 : {
2665 0 : CPLError(CE_Failure, CPLE_AppDefined,
2666 : "Gridding failed at X position %lu, Y position %lu",
2667 : static_cast<long unsigned int>(nXPoint),
2668 : static_cast<long unsigned int>(nYPoint));
2669 0 : *psJob->pbStop = TRUE;
2670 0 : if (pfnProgress != nullptr)
2671 0 : pfnProgress(psJob); // To notify the main thread.
2672 0 : break;
2673 : }
2674 : }
2675 :
2676 16463 : GDALCopyWords(padfValues, GDT_Float64, sizeof(double),
2677 16463 : pabyData + nYPoint * nLineSpace, eType, nDataTypeSize,
2678 : nXSize);
2679 :
2680 20837 : if (*psJob->pbStop || (pfnProgress != nullptr && pfnProgress(psJob)))
2681 0 : break;
2682 : }
2683 :
2684 714 : CPLFree(padfValues);
2685 : }
2686 :
2687 : /************************************************************************/
2688 : /* GDALGridContextCreate() */
2689 : /************************************************************************/
2690 :
2691 : struct GDALGridContext
2692 : {
2693 : GDALGridAlgorithm eAlgorithm;
2694 : void *poOptions;
2695 : GDALGridFunction pfnGDALGridMethod;
2696 :
2697 : GUInt32 nPoints;
2698 : GDALGridPoint *pasGridPoints;
2699 : GDALGridXYArrays sXYArrays;
2700 :
2701 : GDALGridExtraParameters sExtraParameters;
2702 : double *padfX;
2703 : double *padfY;
2704 : double *padfZ;
2705 : bool bFreePadfXYZArrays;
2706 :
2707 : CPLWorkerThreadPool *poWorkerThreadPool;
2708 : };
2709 :
2710 : static void GDALGridContextCreateQuadTree(GDALGridContext *psContext);
2711 :
2712 : /**
2713 : * Creates a context to do regular gridding from the scattered data.
2714 : *
2715 : * This function takes the arrays of X and Y coordinates and corresponding Z
2716 : * values as input to prepare computation of regular grid (or call it a raster)
2717 : * from these scattered data.
2718 : *
2719 : * On Intel/AMD i386/x86_64 architectures, some
2720 : * gridding methods will be optimized with SSE instructions (provided GDAL
2721 : * has been compiled with such support, and it is available at runtime).
2722 : * Currently, only 'invdist' algorithm with default parameters has an optimized
2723 : * implementation.
2724 : * This can provide substantial speed-up, but sometimes at the expense of
2725 : * reduced floating point precision. This can be disabled by setting the
2726 : * GDAL_USE_SSE configuration option to NO.
2727 : * A further optimized version can use the AVX
2728 : * instruction set. This can be disabled by setting the GDAL_USE_AVX
2729 : * configuration option to NO.
2730 : *
2731 : * It is possible to set the GDAL_NUM_THREADS
2732 : * configuration option to parallelize the processing. The value to set is
2733 : * the number of worker threads, or ALL_CPUS to use all the cores/CPUs of the
2734 : * computer (default value).
2735 : *
2736 : * @param eAlgorithm Gridding method.
2737 : * @param poOptions Options to control chosen gridding method.
2738 : * @param nPoints Number of elements in input arrays.
2739 : * @param padfX Input array of X coordinates.
2740 : * @param padfY Input array of Y coordinates.
2741 : * @param padfZ Input array of Z values.
2742 : * @param bCallerWillKeepPointArraysAlive Whether the provided padfX, padfY,
2743 : * padfZ arrays will still be "alive" during the calls to
2744 : * GDALGridContextProcess(). Setting to TRUE prevent them from being
2745 : * duplicated in the context. If unsure, set to FALSE.
2746 : *
2747 : * @return the context (to be freed with GDALGridContextFree()) or NULL in case
2748 : * or error.
2749 : *
2750 : */
2751 :
2752 183 : GDALGridContext *GDALGridContextCreate(GDALGridAlgorithm eAlgorithm,
2753 : const void *poOptions, GUInt32 nPoints,
2754 : const double *padfX, const double *padfY,
2755 : const double *padfZ,
2756 : int bCallerWillKeepPointArraysAlive)
2757 : {
2758 183 : CPLAssert(poOptions);
2759 183 : CPLAssert(padfX);
2760 183 : CPLAssert(padfY);
2761 183 : CPLAssert(padfZ);
2762 183 : bool bCreateQuadTree = false;
2763 :
2764 : const unsigned int nPointCountThreshold =
2765 183 : atoi(CPLGetConfigOption("GDAL_GRID_POINT_COUNT_THRESHOLD", "100"));
2766 :
2767 : // Starting address aligned on 32-byte boundary for AVX.
2768 183 : float *pafXAligned = nullptr;
2769 183 : float *pafYAligned = nullptr;
2770 183 : float *pafZAligned = nullptr;
2771 :
2772 183 : void *poOptionsNew = nullptr;
2773 :
2774 183 : GDALGridFunction pfnGDALGridMethod = nullptr;
2775 :
2776 183 : switch (eAlgorithm)
2777 : {
2778 44 : case GGA_InverseDistanceToAPower:
2779 : {
2780 44 : const auto poOptionsOld =
2781 : static_cast<const GDALGridInverseDistanceToAPowerOptions *>(
2782 : poOptions);
2783 44 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
2784 : {
2785 0 : CPLError(CE_Failure, CPLE_AppDefined,
2786 : "Wrong value of nSizeOfStructure member");
2787 0 : return nullptr;
2788 : }
2789 : poOptionsNew =
2790 44 : CPLMalloc(sizeof(GDALGridInverseDistanceToAPowerOptions));
2791 44 : memcpy(poOptionsNew, poOptions,
2792 : sizeof(GDALGridInverseDistanceToAPowerOptions));
2793 :
2794 44 : const GDALGridInverseDistanceToAPowerOptions *const poPower =
2795 : static_cast<const GDALGridInverseDistanceToAPowerOptions *>(
2796 : poOptions);
2797 44 : if (poPower->dfRadius1 == 0.0 && poPower->dfRadius2 == 0.0)
2798 : {
2799 35 : const double dfPower = poPower->dfPower;
2800 35 : const double dfSmoothing = poPower->dfSmoothing;
2801 :
2802 35 : pfnGDALGridMethod = GDALGridInverseDistanceToAPowerNoSearch;
2803 35 : if (dfPower == 2.0 && dfSmoothing == 0.0)
2804 : {
2805 : #ifdef HAVE_AVX_AT_COMPILE_TIME
2806 :
2807 33 : if (CPLTestBool(
2808 60 : CPLGetConfigOption("GDAL_USE_AVX", "YES")) &&
2809 27 : CPLHaveRuntimeAVX())
2810 : {
2811 : pafXAligned = static_cast<float *>(
2812 27 : VSI_MALLOC_ALIGNED_AUTO_VERBOSE(sizeof(float) *
2813 : nPoints));
2814 : pafYAligned = static_cast<float *>(
2815 27 : VSI_MALLOC_ALIGNED_AUTO_VERBOSE(sizeof(float) *
2816 : nPoints));
2817 : pafZAligned = static_cast<float *>(
2818 27 : VSI_MALLOC_ALIGNED_AUTO_VERBOSE(sizeof(float) *
2819 : nPoints));
2820 27 : if (pafXAligned != nullptr && pafYAligned != nullptr &&
2821 : pafZAligned != nullptr)
2822 : {
2823 27 : CPLDebug("GDAL_GRID",
2824 : "Using AVX optimized version");
2825 27 : pfnGDALGridMethod =
2826 : GDALGridInverseDistanceToAPower2NoSmoothingNoSearchAVX;
2827 1335 : for (GUInt32 i = 0; i < nPoints; i++)
2828 : {
2829 1308 : pafXAligned[i] = static_cast<float>(padfX[i]);
2830 1308 : pafYAligned[i] = static_cast<float>(padfY[i]);
2831 1308 : pafZAligned[i] = static_cast<float>(padfZ[i]);
2832 27 : }
2833 : }
2834 : else
2835 : {
2836 0 : VSIFree(pafXAligned);
2837 0 : VSIFree(pafYAligned);
2838 0 : VSIFree(pafZAligned);
2839 0 : pafXAligned = nullptr;
2840 0 : pafYAligned = nullptr;
2841 0 : pafZAligned = nullptr;
2842 : }
2843 : }
2844 : #endif
2845 :
2846 : #ifdef HAVE_SSE_AT_COMPILE_TIME
2847 :
2848 39 : if (pafXAligned == nullptr &&
2849 6 : CPLTestBool(CPLGetConfigOption("GDAL_USE_SSE", "YES"))
2850 : #if !defined(USE_NEON_OPTIMIZATIONS)
2851 39 : && CPLHaveRuntimeSSE()
2852 : #endif
2853 : )
2854 : {
2855 : pafXAligned = static_cast<float *>(
2856 3 : VSI_MALLOC_ALIGNED_AUTO_VERBOSE(sizeof(float) *
2857 : nPoints));
2858 : pafYAligned = static_cast<float *>(
2859 3 : VSI_MALLOC_ALIGNED_AUTO_VERBOSE(sizeof(float) *
2860 : nPoints));
2861 : pafZAligned = static_cast<float *>(
2862 3 : VSI_MALLOC_ALIGNED_AUTO_VERBOSE(sizeof(float) *
2863 : nPoints));
2864 3 : if (pafXAligned != nullptr && pafYAligned != nullptr &&
2865 : pafZAligned != nullptr)
2866 : {
2867 3 : CPLDebug("GDAL_GRID",
2868 : "Using SSE optimized version");
2869 3 : pfnGDALGridMethod =
2870 : GDALGridInverseDistanceToAPower2NoSmoothingNoSearchSSE;
2871 405 : for (GUInt32 i = 0; i < nPoints; i++)
2872 : {
2873 402 : pafXAligned[i] = static_cast<float>(padfX[i]);
2874 402 : pafYAligned[i] = static_cast<float>(padfY[i]);
2875 402 : pafZAligned[i] = static_cast<float>(padfZ[i]);
2876 3 : }
2877 : }
2878 : else
2879 : {
2880 0 : VSIFree(pafXAligned);
2881 0 : VSIFree(pafYAligned);
2882 0 : VSIFree(pafZAligned);
2883 0 : pafXAligned = nullptr;
2884 0 : pafYAligned = nullptr;
2885 0 : pafZAligned = nullptr;
2886 : }
2887 : }
2888 : #endif // HAVE_SSE_AT_COMPILE_TIME
2889 35 : }
2890 : }
2891 : else
2892 : {
2893 9 : pfnGDALGridMethod = GDALGridInverseDistanceToAPower;
2894 : }
2895 44 : break;
2896 : }
2897 22 : case GGA_InverseDistanceToAPowerNearestNeighbor:
2898 : {
2899 22 : const auto poOptionsOld = static_cast<
2900 : const GDALGridInverseDistanceToAPowerNearestNeighborOptions *>(
2901 : poOptions);
2902 22 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
2903 : {
2904 0 : CPLError(CE_Failure, CPLE_AppDefined,
2905 : "Wrong value of nSizeOfStructure member");
2906 0 : return nullptr;
2907 : }
2908 22 : poOptionsNew = CPLMalloc(
2909 : sizeof(GDALGridInverseDistanceToAPowerNearestNeighborOptions));
2910 22 : memcpy(
2911 : poOptionsNew, poOptions,
2912 : sizeof(GDALGridInverseDistanceToAPowerNearestNeighborOptions));
2913 :
2914 22 : if (poOptionsOld->nMinPointsPerQuadrant != 0 ||
2915 12 : poOptionsOld->nMaxPointsPerQuadrant != 0)
2916 : {
2917 12 : pfnGDALGridMethod =
2918 : GDALGridInverseDistanceToAPowerNearestNeighborPerQuadrant;
2919 : }
2920 : else
2921 : {
2922 10 : pfnGDALGridMethod =
2923 : GDALGridInverseDistanceToAPowerNearestNeighbor;
2924 : }
2925 22 : bCreateQuadTree = true;
2926 22 : break;
2927 : }
2928 26 : case GGA_MovingAverage:
2929 : {
2930 26 : const auto poOptionsOld =
2931 : static_cast<const GDALGridMovingAverageOptions *>(poOptions);
2932 26 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
2933 : {
2934 0 : CPLError(CE_Failure, CPLE_AppDefined,
2935 : "Wrong value of nSizeOfStructure member");
2936 0 : return nullptr;
2937 : }
2938 26 : poOptionsNew = CPLMalloc(sizeof(GDALGridMovingAverageOptions));
2939 26 : memcpy(poOptionsNew, poOptions,
2940 : sizeof(GDALGridMovingAverageOptions));
2941 :
2942 26 : if (poOptionsOld->nMinPointsPerQuadrant != 0 ||
2943 18 : poOptionsOld->nMaxPointsPerQuadrant != 0)
2944 : {
2945 9 : pfnGDALGridMethod = GDALGridMovingAveragePerQuadrant;
2946 9 : bCreateQuadTree = true;
2947 : }
2948 : else
2949 : {
2950 17 : pfnGDALGridMethod = GDALGridMovingAverage;
2951 17 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
2952 22 : poOptionsOld->dfAngle == 0.0 &&
2953 5 : (poOptionsOld->dfRadius1 > 0.0 ||
2954 1 : poOptionsOld->dfRadius2 > 0.0));
2955 : }
2956 26 : break;
2957 : }
2958 21 : case GGA_NearestNeighbor:
2959 : {
2960 21 : const auto poOptionsOld =
2961 : static_cast<const GDALGridNearestNeighborOptions *>(poOptions);
2962 21 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
2963 : {
2964 0 : CPLError(CE_Failure, CPLE_AppDefined,
2965 : "Wrong value of nSizeOfStructure member");
2966 0 : return nullptr;
2967 : }
2968 21 : poOptionsNew = CPLMalloc(sizeof(GDALGridNearestNeighborOptions));
2969 21 : memcpy(poOptionsNew, poOptions,
2970 : sizeof(GDALGridNearestNeighborOptions));
2971 :
2972 21 : pfnGDALGridMethod = GDALGridNearestNeighbor;
2973 34 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
2974 34 : poOptionsOld->dfAngle == 0.0 &&
2975 13 : (poOptionsOld->dfRadius1 > 0.0 ||
2976 5 : poOptionsOld->dfRadius2 > 0.0));
2977 21 : break;
2978 : }
2979 18 : case GGA_MetricMinimum:
2980 : {
2981 18 : const auto poOptionsOld =
2982 : static_cast<const GDALGridDataMetricsOptions *>(poOptions);
2983 18 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
2984 : {
2985 0 : CPLError(CE_Failure, CPLE_AppDefined,
2986 : "Wrong value of nSizeOfStructure member");
2987 0 : return nullptr;
2988 : }
2989 18 : poOptionsNew = CPLMalloc(sizeof(GDALGridDataMetricsOptions));
2990 18 : memcpy(poOptionsNew, poOptions, sizeof(GDALGridDataMetricsOptions));
2991 :
2992 18 : if (poOptionsOld->nMinPointsPerQuadrant != 0 ||
2993 12 : poOptionsOld->nMaxPointsPerQuadrant != 0)
2994 : {
2995 7 : pfnGDALGridMethod = GDALGridDataMetricMinimumPerQuadrant;
2996 7 : bCreateQuadTree = true;
2997 : }
2998 : else
2999 : {
3000 11 : pfnGDALGridMethod = GDALGridDataMetricMinimum;
3001 11 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
3002 16 : poOptionsOld->dfAngle == 0.0 &&
3003 5 : (poOptionsOld->dfRadius1 > 0.0 ||
3004 1 : poOptionsOld->dfRadius2 > 0.0));
3005 : }
3006 18 : break;
3007 : }
3008 12 : case GGA_MetricMaximum:
3009 : {
3010 12 : const auto poOptionsOld =
3011 : static_cast<const GDALGridDataMetricsOptions *>(poOptions);
3012 12 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
3013 : {
3014 0 : CPLError(CE_Failure, CPLE_AppDefined,
3015 : "Wrong value of nSizeOfStructure member");
3016 0 : return nullptr;
3017 : }
3018 12 : poOptionsNew = CPLMalloc(sizeof(GDALGridDataMetricsOptions));
3019 12 : memcpy(poOptionsNew, poOptions, sizeof(GDALGridDataMetricsOptions));
3020 :
3021 12 : if (poOptionsOld->nMinPointsPerQuadrant != 0 ||
3022 7 : poOptionsOld->nMaxPointsPerQuadrant != 0)
3023 : {
3024 5 : pfnGDALGridMethod = GDALGridDataMetricMaximumPerQuadrant;
3025 5 : bCreateQuadTree = true;
3026 : }
3027 : else
3028 : {
3029 7 : pfnGDALGridMethod = GDALGridDataMetricMaximum;
3030 7 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
3031 11 : poOptionsOld->dfAngle == 0.0 &&
3032 4 : (poOptionsOld->dfRadius1 > 0.0 ||
3033 1 : poOptionsOld->dfRadius2 > 0.0));
3034 : }
3035 :
3036 12 : break;
3037 : }
3038 9 : case GGA_MetricRange:
3039 : {
3040 9 : const auto poOptionsOld =
3041 : static_cast<const GDALGridDataMetricsOptions *>(poOptions);
3042 9 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
3043 : {
3044 0 : CPLError(CE_Failure, CPLE_AppDefined,
3045 : "Wrong value of nSizeOfStructure member");
3046 0 : return nullptr;
3047 : }
3048 9 : poOptionsNew = CPLMalloc(sizeof(GDALGridDataMetricsOptions));
3049 9 : memcpy(poOptionsNew, poOptions, sizeof(GDALGridDataMetricsOptions));
3050 :
3051 9 : if (poOptionsOld->nMinPointsPerQuadrant != 0 ||
3052 4 : poOptionsOld->nMaxPointsPerQuadrant != 0)
3053 : {
3054 5 : pfnGDALGridMethod = GDALGridDataMetricRangePerQuadrant;
3055 5 : bCreateQuadTree = true;
3056 : }
3057 : else
3058 : {
3059 4 : pfnGDALGridMethod = GDALGridDataMetricRange;
3060 4 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
3061 7 : poOptionsOld->dfAngle == 0.0 &&
3062 3 : (poOptionsOld->dfRadius1 > 0.0 ||
3063 1 : poOptionsOld->dfRadius2 > 0.0));
3064 : }
3065 :
3066 9 : break;
3067 : }
3068 11 : case GGA_MetricCount:
3069 : {
3070 11 : const auto poOptionsOld =
3071 : static_cast<const GDALGridDataMetricsOptions *>(poOptions);
3072 11 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
3073 : {
3074 0 : CPLError(CE_Failure, CPLE_AppDefined,
3075 : "Wrong value of nSizeOfStructure member");
3076 0 : return nullptr;
3077 : }
3078 11 : poOptionsNew = CPLMalloc(sizeof(GDALGridDataMetricsOptions));
3079 11 : memcpy(poOptionsNew, poOptions, sizeof(GDALGridDataMetricsOptions));
3080 :
3081 11 : if (poOptionsOld->nMinPointsPerQuadrant != 0 ||
3082 6 : poOptionsOld->nMaxPointsPerQuadrant != 0)
3083 : {
3084 5 : pfnGDALGridMethod = GDALGridDataMetricCountPerQuadrant;
3085 5 : bCreateQuadTree = true;
3086 : }
3087 : else
3088 : {
3089 6 : pfnGDALGridMethod = GDALGridDataMetricCount;
3090 6 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
3091 11 : poOptionsOld->dfAngle == 0.0 &&
3092 5 : (poOptionsOld->dfRadius1 > 0.0 ||
3093 0 : poOptionsOld->dfRadius2 > 0.0));
3094 : }
3095 :
3096 11 : break;
3097 : }
3098 9 : case GGA_MetricAverageDistance:
3099 : {
3100 9 : const auto poOptionsOld =
3101 : static_cast<const GDALGridDataMetricsOptions *>(poOptions);
3102 9 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
3103 : {
3104 0 : CPLError(CE_Failure, CPLE_AppDefined,
3105 : "Wrong value of nSizeOfStructure member");
3106 0 : return nullptr;
3107 : }
3108 9 : poOptionsNew = CPLMalloc(sizeof(GDALGridDataMetricsOptions));
3109 9 : memcpy(poOptionsNew, poOptions, sizeof(GDALGridDataMetricsOptions));
3110 :
3111 9 : if (poOptionsOld->nMinPointsPerQuadrant != 0 ||
3112 4 : poOptionsOld->nMaxPointsPerQuadrant != 0)
3113 : {
3114 5 : pfnGDALGridMethod =
3115 : GDALGridDataMetricAverageDistancePerQuadrant;
3116 5 : bCreateQuadTree = true;
3117 : }
3118 : else
3119 : {
3120 4 : pfnGDALGridMethod = GDALGridDataMetricAverageDistance;
3121 4 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
3122 7 : poOptionsOld->dfAngle == 0.0 &&
3123 3 : (poOptionsOld->dfRadius1 > 0.0 ||
3124 1 : poOptionsOld->dfRadius2 > 0.0));
3125 : }
3126 :
3127 9 : break;
3128 : }
3129 4 : case GGA_MetricAverageDistancePts:
3130 : {
3131 4 : const auto poOptionsOld =
3132 : static_cast<const GDALGridDataMetricsOptions *>(poOptions);
3133 4 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
3134 : {
3135 0 : CPLError(CE_Failure, CPLE_AppDefined,
3136 : "Wrong value of nSizeOfStructure member");
3137 0 : return nullptr;
3138 : }
3139 4 : poOptionsNew = CPLMalloc(sizeof(GDALGridDataMetricsOptions));
3140 4 : memcpy(poOptionsNew, poOptions, sizeof(GDALGridDataMetricsOptions));
3141 :
3142 4 : pfnGDALGridMethod = GDALGridDataMetricAverageDistancePts;
3143 7 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
3144 6 : poOptionsOld->dfAngle == 0.0 &&
3145 2 : (poOptionsOld->dfRadius1 > 0.0 ||
3146 0 : poOptionsOld->dfRadius2 > 0.0));
3147 :
3148 4 : break;
3149 : }
3150 7 : case GGA_Linear:
3151 : {
3152 7 : const auto poOptionsOld =
3153 : static_cast<const GDALGridLinearOptions *>(poOptions);
3154 7 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
3155 : {
3156 0 : CPLError(CE_Failure, CPLE_AppDefined,
3157 : "Wrong value of nSizeOfStructure member");
3158 0 : return nullptr;
3159 : }
3160 7 : poOptionsNew = CPLMalloc(sizeof(GDALGridLinearOptions));
3161 7 : memcpy(poOptionsNew, poOptions, sizeof(GDALGridLinearOptions));
3162 :
3163 7 : pfnGDALGridMethod = GDALGridLinear;
3164 7 : break;
3165 : }
3166 0 : default:
3167 0 : CPLError(CE_Failure, CPLE_IllegalArg,
3168 : "GDAL does not support gridding method %d", eAlgorithm);
3169 0 : return nullptr;
3170 : }
3171 :
3172 183 : if (pafXAligned == nullptr && !bCallerWillKeepPointArraysAlive)
3173 : {
3174 : double *padfXNew =
3175 0 : static_cast<double *>(VSI_MALLOC2_VERBOSE(nPoints, sizeof(double)));
3176 : double *padfYNew =
3177 0 : static_cast<double *>(VSI_MALLOC2_VERBOSE(nPoints, sizeof(double)));
3178 : double *padfZNew =
3179 0 : static_cast<double *>(VSI_MALLOC2_VERBOSE(nPoints, sizeof(double)));
3180 0 : if (padfXNew == nullptr || padfYNew == nullptr || padfZNew == nullptr)
3181 : {
3182 0 : VSIFree(padfXNew);
3183 0 : VSIFree(padfYNew);
3184 0 : VSIFree(padfZNew);
3185 0 : CPLFree(poOptionsNew);
3186 0 : return nullptr;
3187 : }
3188 0 : memcpy(padfXNew, padfX, nPoints * sizeof(double));
3189 0 : memcpy(padfYNew, padfY, nPoints * sizeof(double));
3190 0 : memcpy(padfZNew, padfZ, nPoints * sizeof(double));
3191 0 : padfX = padfXNew;
3192 0 : padfY = padfYNew;
3193 0 : padfZ = padfZNew;
3194 : }
3195 : GDALGridContext *psContext =
3196 183 : static_cast<GDALGridContext *>(CPLCalloc(1, sizeof(GDALGridContext)));
3197 183 : psContext->eAlgorithm = eAlgorithm;
3198 183 : psContext->poOptions = poOptionsNew;
3199 183 : psContext->pfnGDALGridMethod = pfnGDALGridMethod;
3200 183 : psContext->nPoints = nPoints;
3201 183 : psContext->pasGridPoints = nullptr;
3202 183 : psContext->sXYArrays.padfX = padfX;
3203 183 : psContext->sXYArrays.padfY = padfY;
3204 183 : psContext->sExtraParameters.hQuadTree = nullptr;
3205 183 : psContext->sExtraParameters.dfInitialSearchRadius = 0.0;
3206 183 : psContext->sExtraParameters.pafX = pafXAligned;
3207 183 : psContext->sExtraParameters.pafY = pafYAligned;
3208 183 : psContext->sExtraParameters.pafZ = pafZAligned;
3209 183 : psContext->sExtraParameters.psTriangulation = nullptr;
3210 183 : psContext->sExtraParameters.nInitialFacetIdx = 0;
3211 183 : psContext->padfX = pafXAligned ? nullptr : const_cast<double *>(padfX);
3212 183 : psContext->padfY = pafXAligned ? nullptr : const_cast<double *>(padfY);
3213 183 : psContext->padfZ = pafXAligned ? nullptr : const_cast<double *>(padfZ);
3214 183 : psContext->bFreePadfXYZArrays =
3215 183 : pafXAligned ? false : !bCallerWillKeepPointArraysAlive;
3216 :
3217 : /* -------------------------------------------------------------------- */
3218 : /* Create quadtree if requested and possible. */
3219 : /* -------------------------------------------------------------------- */
3220 183 : if (bCreateQuadTree)
3221 : {
3222 88 : GDALGridContextCreateQuadTree(psContext);
3223 88 : if (psContext->sExtraParameters.hQuadTree == nullptr &&
3224 0 : (eAlgorithm == GGA_InverseDistanceToAPowerNearestNeighbor ||
3225 : pfnGDALGridMethod == GDALGridMovingAveragePerQuadrant))
3226 : {
3227 : // shouldn't happen unless memory allocation failure occurs
3228 0 : GDALGridContextFree(psContext);
3229 0 : return nullptr;
3230 : }
3231 : }
3232 :
3233 : /* -------------------------------------------------------------------- */
3234 : /* Pre-compute extra parameters in GDALGridExtraParameters */
3235 : /* -------------------------------------------------------------------- */
3236 183 : if (eAlgorithm == GGA_InverseDistanceToAPowerNearestNeighbor)
3237 : {
3238 22 : const double dfPower =
3239 : static_cast<
3240 : const GDALGridInverseDistanceToAPowerNearestNeighborOptions *>(
3241 : poOptions)
3242 : ->dfPower;
3243 22 : psContext->sExtraParameters.dfPowerDiv2PreComp = dfPower / 2;
3244 :
3245 22 : const double dfRadius =
3246 : static_cast<
3247 : const GDALGridInverseDistanceToAPowerNearestNeighborOptions *>(
3248 : poOptions)
3249 : ->dfRadius;
3250 22 : psContext->sExtraParameters.dfRadiusPower2PreComp = pow(dfRadius, 2);
3251 : }
3252 :
3253 183 : if (eAlgorithm == GGA_Linear)
3254 : {
3255 7 : psContext->sExtraParameters.psTriangulation =
3256 7 : GDALTriangulationCreateDelaunay(nPoints, padfX, padfY);
3257 7 : if (psContext->sExtraParameters.psTriangulation == nullptr)
3258 : {
3259 0 : GDALGridContextFree(psContext);
3260 0 : return nullptr;
3261 : }
3262 7 : GDALTriangulationComputeBarycentricCoefficients(
3263 : psContext->sExtraParameters.psTriangulation, padfX, padfY);
3264 : }
3265 :
3266 : /* -------------------------------------------------------------------- */
3267 : /* Start thread pool. */
3268 : /* -------------------------------------------------------------------- */
3269 183 : const char *pszThreads = CPLGetConfigOption("GDAL_NUM_THREADS", "ALL_CPUS");
3270 183 : int nThreads = 0;
3271 183 : if (EQUAL(pszThreads, "ALL_CPUS"))
3272 181 : nThreads = CPLGetNumCPUs();
3273 : else
3274 2 : nThreads = atoi(pszThreads);
3275 183 : if (nThreads > 128)
3276 0 : nThreads = 128;
3277 183 : if (nThreads > 1)
3278 : {
3279 181 : psContext->poWorkerThreadPool = new CPLWorkerThreadPool();
3280 181 : if (!psContext->poWorkerThreadPool->Setup(nThreads, nullptr, nullptr))
3281 : {
3282 0 : delete psContext->poWorkerThreadPool;
3283 0 : psContext->poWorkerThreadPool = nullptr;
3284 : }
3285 : else
3286 : {
3287 181 : CPLDebug("GDAL_GRID", "Using %d threads", nThreads);
3288 : }
3289 : }
3290 : else
3291 2 : psContext->poWorkerThreadPool = nullptr;
3292 :
3293 183 : return psContext;
3294 : }
3295 :
3296 : /************************************************************************/
3297 : /* GDALGridContextCreateQuadTree() */
3298 : /************************************************************************/
3299 :
3300 94 : void GDALGridContextCreateQuadTree(GDALGridContext *psContext)
3301 : {
3302 94 : const GUInt32 nPoints = psContext->nPoints;
3303 94 : psContext->pasGridPoints = static_cast<GDALGridPoint *>(
3304 94 : VSI_MALLOC2_VERBOSE(nPoints, sizeof(GDALGridPoint)));
3305 94 : if (psContext->pasGridPoints != nullptr)
3306 : {
3307 94 : const double *const padfX = psContext->padfX;
3308 94 : const double *const padfY = psContext->padfY;
3309 :
3310 : // Determine point extents.
3311 : CPLRectObj sRect;
3312 94 : sRect.minx = padfX[0];
3313 94 : sRect.miny = padfY[0];
3314 94 : sRect.maxx = padfX[0];
3315 94 : sRect.maxy = padfY[0];
3316 28131 : for (GUInt32 i = 1; i < nPoints; i++)
3317 : {
3318 28037 : if (padfX[i] < sRect.minx)
3319 40 : sRect.minx = padfX[i];
3320 28037 : if (padfY[i] < sRect.miny)
3321 788 : sRect.miny = padfY[i];
3322 28037 : if (padfX[i] > sRect.maxx)
3323 794 : sRect.maxx = padfX[i];
3324 28037 : if (padfY[i] > sRect.maxy)
3325 30 : sRect.maxy = padfY[i];
3326 : }
3327 :
3328 : // Initial value for search radius is the typical dimension of a
3329 : // "pixel" of the point array (assuming rather uniform distribution).
3330 94 : psContext->sExtraParameters.dfInitialSearchRadius = sqrt(
3331 94 : (sRect.maxx - sRect.minx) * (sRect.maxy - sRect.miny) / nPoints);
3332 :
3333 94 : psContext->sExtraParameters.hQuadTree =
3334 94 : CPLQuadTreeCreate(&sRect, GDALGridGetPointBounds);
3335 :
3336 28225 : for (GUInt32 i = 0; i < nPoints; i++)
3337 : {
3338 28131 : psContext->pasGridPoints[i].psXYArrays = &(psContext->sXYArrays);
3339 28131 : psContext->pasGridPoints[i].i = i;
3340 28131 : CPLQuadTreeInsert(psContext->sExtraParameters.hQuadTree,
3341 28131 : psContext->pasGridPoints + i);
3342 : }
3343 : }
3344 94 : }
3345 :
3346 : /************************************************************************/
3347 : /* GDALGridContextFree() */
3348 : /************************************************************************/
3349 :
3350 : /**
3351 : * Free a context used created by GDALGridContextCreate()
3352 : *
3353 : * @param psContext the context.
3354 : *
3355 : */
3356 183 : void GDALGridContextFree(GDALGridContext *psContext)
3357 : {
3358 183 : if (psContext)
3359 : {
3360 183 : CPLFree(psContext->poOptions);
3361 183 : CPLFree(psContext->pasGridPoints);
3362 183 : if (psContext->sExtraParameters.hQuadTree != nullptr)
3363 94 : CPLQuadTreeDestroy(psContext->sExtraParameters.hQuadTree);
3364 183 : if (psContext->bFreePadfXYZArrays)
3365 : {
3366 0 : CPLFree(psContext->padfX);
3367 0 : CPLFree(psContext->padfY);
3368 0 : CPLFree(psContext->padfZ);
3369 : }
3370 183 : VSIFreeAligned(psContext->sExtraParameters.pafX);
3371 183 : VSIFreeAligned(psContext->sExtraParameters.pafY);
3372 183 : VSIFreeAligned(psContext->sExtraParameters.pafZ);
3373 183 : if (psContext->sExtraParameters.psTriangulation)
3374 7 : GDALTriangulationFree(psContext->sExtraParameters.psTriangulation);
3375 183 : delete psContext->poWorkerThreadPool;
3376 183 : CPLFree(psContext);
3377 : }
3378 183 : }
3379 :
3380 : /************************************************************************/
3381 : /* GDALGridContextProcess() */
3382 : /************************************************************************/
3383 :
3384 : /**
3385 : * Do the gridding of a window of a raster.
3386 : *
3387 : * This function takes the gridding context as input to prepare computation
3388 : * of regular grid (or call it a raster) from these scattered data.
3389 : * You should supply the extent of the output grid and allocate array
3390 : * sufficient to hold such a grid.
3391 : *
3392 : * @param psContext Gridding context.
3393 : * @param dfXMin Lowest X border of output grid.
3394 : * @param dfXMax Highest X border of output grid.
3395 : * @param dfYMin Lowest Y border of output grid.
3396 : * @param dfYMax Highest Y border of output grid.
3397 : * @param nXSize Number of columns in output grid.
3398 : * @param nYSize Number of rows in output grid.
3399 : * @param eType Data type of output array.
3400 : * @param pData Pointer to array where the computed grid will be stored.
3401 : * @param pfnProgress a GDALProgressFunc() compatible callback function for
3402 : * reporting progress or NULL.
3403 : * @param pProgressArg argument to be passed to pfnProgress. May be NULL.
3404 : *
3405 : * @return CE_None on success or CE_Failure if something goes wrong.
3406 : *
3407 : */
3408 :
3409 183 : CPLErr GDALGridContextProcess(GDALGridContext *psContext, double dfXMin,
3410 : double dfXMax, double dfYMin, double dfYMax,
3411 : GUInt32 nXSize, GUInt32 nYSize,
3412 : GDALDataType eType, void *pData,
3413 : GDALProgressFunc pfnProgress, void *pProgressArg)
3414 : {
3415 183 : CPLAssert(psContext);
3416 183 : CPLAssert(pData);
3417 :
3418 183 : if (nXSize == 0 || nYSize == 0)
3419 : {
3420 0 : CPLError(CE_Failure, CPLE_IllegalArg,
3421 : "Output raster dimensions should have non-zero size.");
3422 0 : return CE_Failure;
3423 : }
3424 :
3425 183 : const double dfDeltaX = (dfXMax - dfXMin) / nXSize;
3426 183 : const double dfDeltaY = (dfYMax - dfYMin) / nYSize;
3427 :
3428 : // For linear, check if we will need to fallback to nearest neighbour
3429 : // by sampling along the edges. If all points on edges are within
3430 : // triangles, then interior points will also be.
3431 183 : if (psContext->eAlgorithm == GGA_Linear &&
3432 7 : psContext->sExtraParameters.hQuadTree == nullptr)
3433 : {
3434 7 : bool bNeedNearest = false;
3435 7 : int nStartLeft = 0;
3436 7 : int nStartRight = 0;
3437 7 : const double dfXPointMin = dfXMin + (0 + 0.5) * dfDeltaX;
3438 7 : const double dfXPointMax = dfXMin + (nXSize - 1 + 0.5) * dfDeltaX;
3439 269 : for (GUInt32 nYPoint = 0; !bNeedNearest && nYPoint < nYSize; nYPoint++)
3440 : {
3441 262 : const double dfYPoint = dfYMin + (nYPoint + 0.5) * dfDeltaY;
3442 :
3443 262 : if (!GDALTriangulationFindFacetDirected(
3444 262 : psContext->sExtraParameters.psTriangulation, nStartLeft,
3445 : dfXPointMin, dfYPoint, &nStartLeft))
3446 : {
3447 6 : bNeedNearest = true;
3448 : }
3449 262 : if (!GDALTriangulationFindFacetDirected(
3450 262 : psContext->sExtraParameters.psTriangulation, nStartRight,
3451 : dfXPointMax, dfYPoint, &nStartRight))
3452 : {
3453 6 : bNeedNearest = true;
3454 : }
3455 : }
3456 7 : int nStartTop = 0;
3457 7 : int nStartBottom = 0;
3458 7 : const double dfYPointMin = dfYMin + (0 + 0.5) * dfDeltaY;
3459 7 : const double dfYPointMax = dfYMin + (nYSize - 1 + 0.5) * dfDeltaY;
3460 261 : for (GUInt32 nXPoint = 1; !bNeedNearest && nXPoint + 1 < nXSize;
3461 : nXPoint++)
3462 : {
3463 254 : const double dfXPoint = dfXMin + (nXPoint + 0.5) * dfDeltaX;
3464 :
3465 254 : if (!GDALTriangulationFindFacetDirected(
3466 254 : psContext->sExtraParameters.psTriangulation, nStartTop,
3467 : dfXPoint, dfYPointMin, &nStartTop))
3468 : {
3469 0 : bNeedNearest = true;
3470 : }
3471 254 : if (!GDALTriangulationFindFacetDirected(
3472 254 : psContext->sExtraParameters.psTriangulation, nStartBottom,
3473 : dfXPoint, dfYPointMax, &nStartBottom))
3474 : {
3475 0 : bNeedNearest = true;
3476 : }
3477 : }
3478 7 : if (bNeedNearest)
3479 : {
3480 6 : CPLDebug("GDAL_GRID", "Will need nearest neighbour");
3481 6 : GDALGridContextCreateQuadTree(psContext);
3482 : }
3483 : }
3484 :
3485 183 : int nCounter = 0;
3486 183 : volatile int bStop = FALSE;
3487 : GDALGridJob sJob;
3488 183 : sJob.nYStart = 0;
3489 183 : sJob.pabyData = static_cast<GByte *>(pData);
3490 183 : sJob.nYStep = 1;
3491 183 : sJob.nXSize = nXSize;
3492 183 : sJob.nYSize = nYSize;
3493 183 : sJob.dfXMin = dfXMin;
3494 183 : sJob.dfYMin = dfYMin;
3495 183 : sJob.dfDeltaX = dfDeltaX;
3496 183 : sJob.dfDeltaY = dfDeltaY;
3497 183 : sJob.nPoints = psContext->nPoints;
3498 183 : sJob.padfX = psContext->padfX;
3499 183 : sJob.padfY = psContext->padfY;
3500 183 : sJob.padfZ = psContext->padfZ;
3501 183 : sJob.poOptions = psContext->poOptions;
3502 183 : sJob.pfnGDALGridMethod = psContext->pfnGDALGridMethod;
3503 183 : sJob.psExtraParameters = &psContext->sExtraParameters;
3504 183 : sJob.pfnProgress = nullptr;
3505 183 : sJob.eType = eType;
3506 183 : sJob.pfnRealProgress = pfnProgress;
3507 183 : sJob.pRealProgressArg = pProgressArg;
3508 183 : sJob.nCounterSingleThreaded = 0;
3509 183 : sJob.pnCounter = &nCounter;
3510 183 : sJob.pbStop = &bStop;
3511 183 : sJob.hCond = nullptr;
3512 183 : sJob.hCondMutex = nullptr;
3513 :
3514 183 : if (psContext->poWorkerThreadPool == nullptr)
3515 : {
3516 2 : if (sJob.pfnRealProgress != nullptr &&
3517 2 : sJob.pfnRealProgress != GDALDummyProgress)
3518 : {
3519 2 : sJob.pfnProgress = GDALGridProgressMonoThread;
3520 : }
3521 :
3522 2 : GDALGridJobProcess(&sJob);
3523 : }
3524 : else
3525 : {
3526 181 : int nThreads = psContext->poWorkerThreadPool->GetThreadCount();
3527 : GDALGridJob *pasJobs = static_cast<GDALGridJob *>(
3528 181 : CPLMalloc(sizeof(GDALGridJob) * nThreads));
3529 :
3530 181 : sJob.nYStep = nThreads;
3531 181 : sJob.hCondMutex = CPLCreateMutex(); /* and implicitly take the mutex */
3532 181 : sJob.hCond = CPLCreateCond();
3533 181 : sJob.pfnProgress = GDALGridProgressMultiThread;
3534 :
3535 : /* --------------------------------------------------------------------
3536 : */
3537 : /* Start threads. */
3538 : /* --------------------------------------------------------------------
3539 : */
3540 905 : for (int i = 0; i < nThreads && !bStop; i++)
3541 : {
3542 724 : memcpy(&pasJobs[i], &sJob, sizeof(GDALGridJob));
3543 724 : pasJobs[i].nYStart = i;
3544 724 : psContext->poWorkerThreadPool->SubmitJob(GDALGridJobProcess,
3545 724 : &pasJobs[i]);
3546 : }
3547 :
3548 : /* --------------------------------------------------------------------
3549 : */
3550 : /* Report progress. */
3551 : /* --------------------------------------------------------------------
3552 : */
3553 7692 : while (*(sJob.pnCounter) < static_cast<int>(nYSize) && !bStop)
3554 : {
3555 7511 : CPLCondWait(sJob.hCond, sJob.hCondMutex);
3556 :
3557 7511 : int nLocalCounter = *(sJob.pnCounter);
3558 7511 : CPLReleaseMutex(sJob.hCondMutex);
3559 :
3560 15016 : if (pfnProgress != nullptr &&
3561 7505 : !pfnProgress(nLocalCounter / static_cast<double>(nYSize), "",
3562 : pProgressArg))
3563 : {
3564 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
3565 0 : bStop = TRUE;
3566 : }
3567 :
3568 7511 : CPLAcquireMutex(sJob.hCondMutex, 1.0);
3569 : }
3570 :
3571 : // Release mutex before joining threads, otherwise they will dead-lock
3572 : // forever in GDALGridProgressMultiThread().
3573 181 : CPLReleaseMutex(sJob.hCondMutex);
3574 :
3575 : /* --------------------------------------------------------------------
3576 : */
3577 : /* Wait for all threads to complete and finish. */
3578 : /* --------------------------------------------------------------------
3579 : */
3580 181 : psContext->poWorkerThreadPool->WaitCompletion();
3581 :
3582 181 : CPLFree(pasJobs);
3583 181 : CPLDestroyCond(sJob.hCond);
3584 181 : CPLDestroyMutex(sJob.hCondMutex);
3585 : }
3586 :
3587 183 : return bStop ? CE_Failure : CE_None;
3588 : }
3589 :
3590 : /************************************************************************/
3591 : /* GDALGridCreate() */
3592 : /************************************************************************/
3593 :
3594 : /**
3595 : * Create regular grid from the scattered data.
3596 : *
3597 : * This function takes the arrays of X and Y coordinates and corresponding Z
3598 : * values as input and computes regular grid (or call it a raster) from these
3599 : * scattered data. You should supply geometry and extent of the output grid
3600 : * and allocate array sufficient to hold such a grid.
3601 : *
3602 : * It is possible to set the GDAL_NUM_THREADS
3603 : * configuration option to parallelize the processing. The value to set is
3604 : * the number of worker threads, or ALL_CPUS to use all the cores/CPUs of the
3605 : * computer (default value).
3606 : *
3607 : * On Intel/AMD i386/x86_64 architectures, some
3608 : * gridding methods will be optimized with SSE instructions (provided GDAL
3609 : * has been compiled with such support, and it is available at runtime).
3610 : * Currently, only 'invdist' algorithm with default parameters has an optimized
3611 : * implementation.
3612 : * This can provide substantial speed-up, but sometimes at the expense of
3613 : * reduced floating point precision. This can be disabled by setting the
3614 : * GDAL_USE_SSE configuration option to NO.
3615 : * A further optimized version can use the AVX
3616 : * instruction set. This can be disabled by setting the GDAL_USE_AVX
3617 : * configuration option to NO.
3618 : *
3619 : * Note: it will be more efficient to use GDALGridContextCreate(),
3620 : * GDALGridContextProcess() and GDALGridContextFree() when doing repeated
3621 : * gridding operations with the same algorithm, parameters and points, and
3622 : * moving the window in the output grid.
3623 : *
3624 : * @param eAlgorithm Gridding method.
3625 : * @param poOptions Options to control chosen gridding method.
3626 : * @param nPoints Number of elements in input arrays.
3627 : * @param padfX Input array of X coordinates.
3628 : * @param padfY Input array of Y coordinates.
3629 : * @param padfZ Input array of Z values.
3630 : * @param dfXMin Lowest X border of output grid.
3631 : * @param dfXMax Highest X border of output grid.
3632 : * @param dfYMin Lowest Y border of output grid.
3633 : * @param dfYMax Highest Y border of output grid.
3634 : * @param nXSize Number of columns in output grid.
3635 : * @param nYSize Number of rows in output grid.
3636 : * @param eType Data type of output array.
3637 : * @param pData Pointer to array where the computed grid will be stored.
3638 : * @param pfnProgress a GDALProgressFunc() compatible callback function for
3639 : * reporting progress or NULL.
3640 : * @param pProgressArg argument to be passed to pfnProgress. May be NULL.
3641 : *
3642 : * @return CE_None on success or CE_Failure if something goes wrong.
3643 : */
3644 :
3645 5 : CPLErr GDALGridCreate(GDALGridAlgorithm eAlgorithm, const void *poOptions,
3646 : GUInt32 nPoints, const double *padfX, const double *padfY,
3647 : const double *padfZ, double dfXMin, double dfXMax,
3648 : double dfYMin, double dfYMax, GUInt32 nXSize,
3649 : GUInt32 nYSize, GDALDataType eType, void *pData,
3650 : GDALProgressFunc pfnProgress, void *pProgressArg)
3651 : {
3652 5 : GDALGridContext *psContext = GDALGridContextCreate(
3653 : eAlgorithm, poOptions, nPoints, padfX, padfY, padfZ, TRUE);
3654 5 : CPLErr eErr = CE_Failure;
3655 5 : if (psContext)
3656 : {
3657 5 : eErr = GDALGridContextProcess(psContext, dfXMin, dfXMax, dfYMin, dfYMax,
3658 : nXSize, nYSize, eType, pData, pfnProgress,
3659 : pProgressArg);
3660 : }
3661 :
3662 5 : GDALGridContextFree(psContext);
3663 5 : return eErr;
3664 : }
3665 :
3666 : /************************************************************************/
3667 : /* GDALGridParseAlgorithmAndOptions() */
3668 : /************************************************************************/
3669 :
3670 : /** Translates mnemonic gridding algorithm names into GDALGridAlgorithm code,
3671 : * parse control parameters and assign defaults.
3672 : */
3673 383 : CPLErr GDALGridParseAlgorithmAndOptions(const char *pszAlgorithm,
3674 : GDALGridAlgorithm *peAlgorithm,
3675 : void **ppOptions)
3676 : {
3677 383 : CPLAssert(pszAlgorithm);
3678 383 : CPLAssert(peAlgorithm);
3679 383 : CPLAssert(ppOptions);
3680 :
3681 383 : *ppOptions = nullptr;
3682 :
3683 383 : char **papszParams = CSLTokenizeString2(pszAlgorithm, ":", FALSE);
3684 :
3685 383 : if (CSLCount(papszParams) < 1)
3686 : {
3687 0 : CSLDestroy(papszParams);
3688 0 : return CE_Failure;
3689 : }
3690 :
3691 383 : if (EQUAL(papszParams[0], szAlgNameInvDist))
3692 : {
3693 495 : if (CSLFetchNameValue(papszParams, "min_points_per_quadrant") ||
3694 246 : CSLFetchNameValue(papszParams, "max_points_per_quadrant"))
3695 : {
3696 : // Remap invdist to invdistnn if per quadrant is required
3697 4 : if (CSLFetchNameValue(papszParams, "radius") == nullptr)
3698 : {
3699 : const double dfRadius1 =
3700 0 : CPLAtofM(CSLFetchNameValueDef(papszParams, "radius1", "1"));
3701 : const double dfRadius2 =
3702 0 : CPLAtofM(CSLFetchNameValueDef(papszParams, "radius2", "1"));
3703 0 : if (dfRadius1 != dfRadius2)
3704 : {
3705 0 : CPLError(CE_Failure, CPLE_NotSupported,
3706 : "radius1 != radius2 not supported when "
3707 : "min_points_per_quadrant and/or "
3708 : "max_points_per_quadrant is specified");
3709 0 : CSLDestroy(papszParams);
3710 0 : return CE_Failure;
3711 : }
3712 : }
3713 :
3714 4 : if (CPLAtofM(CSLFetchNameValueDef(papszParams, "angle", "0")) != 0)
3715 : {
3716 0 : CPLError(CE_Failure, CPLE_NotSupported,
3717 : "angle != 0 not supported when "
3718 : "min_points_per_quadrant and/or "
3719 : "max_points_per_quadrant is specified");
3720 0 : CSLDestroy(papszParams);
3721 0 : return CE_Failure;
3722 : }
3723 :
3724 4 : char **papszNewParams = CSLAddString(nullptr, "invdistnn");
3725 4 : if (CSLFetchNameValue(papszParams, "radius") == nullptr)
3726 : {
3727 0 : papszNewParams = CSLSetNameValue(
3728 : papszNewParams, "radius",
3729 : CSLFetchNameValueDef(papszParams, "radius1", "1"));
3730 : }
3731 32 : for (const char *pszItem :
3732 : {"radius", "power", "smoothing", "max_points", "min_points",
3733 : "nodata", "min_points_per_quadrant",
3734 36 : "max_points_per_quadrant"})
3735 : {
3736 32 : const char *pszValue = CSLFetchNameValue(papszParams, pszItem);
3737 32 : if (pszValue)
3738 : papszNewParams =
3739 19 : CSLSetNameValue(papszNewParams, pszItem, pszValue);
3740 : }
3741 4 : CSLDestroy(papszParams);
3742 4 : papszParams = papszNewParams;
3743 :
3744 4 : *peAlgorithm = GGA_InverseDistanceToAPowerNearestNeighbor;
3745 : }
3746 : else
3747 : {
3748 245 : *peAlgorithm = GGA_InverseDistanceToAPower;
3749 : }
3750 : }
3751 134 : else if (EQUAL(papszParams[0], szAlgNameInvDistNearestNeighbor))
3752 : {
3753 18 : *peAlgorithm = GGA_InverseDistanceToAPowerNearestNeighbor;
3754 : }
3755 116 : else if (EQUAL(papszParams[0], szAlgNameAverage))
3756 : {
3757 26 : *peAlgorithm = GGA_MovingAverage;
3758 : }
3759 90 : else if (EQUAL(papszParams[0], szAlgNameNearest))
3760 : {
3761 19 : *peAlgorithm = GGA_NearestNeighbor;
3762 : }
3763 71 : else if (EQUAL(papszParams[0], szAlgNameMinimum))
3764 : {
3765 18 : *peAlgorithm = GGA_MetricMinimum;
3766 : }
3767 53 : else if (EQUAL(papszParams[0], szAlgNameMaximum))
3768 : {
3769 12 : *peAlgorithm = GGA_MetricMaximum;
3770 : }
3771 41 : else if (EQUAL(papszParams[0], szAlgNameRange))
3772 : {
3773 9 : *peAlgorithm = GGA_MetricRange;
3774 : }
3775 32 : else if (EQUAL(papszParams[0], szAlgNameCount))
3776 : {
3777 11 : *peAlgorithm = GGA_MetricCount;
3778 : }
3779 21 : else if (EQUAL(papszParams[0], szAlgNameAverageDistance))
3780 : {
3781 9 : *peAlgorithm = GGA_MetricAverageDistance;
3782 : }
3783 12 : else if (EQUAL(papszParams[0], szAlgNameAverageDistancePts))
3784 : {
3785 4 : *peAlgorithm = GGA_MetricAverageDistancePts;
3786 : }
3787 8 : else if (EQUAL(papszParams[0], szAlgNameLinear))
3788 : {
3789 7 : *peAlgorithm = GGA_Linear;
3790 : }
3791 : else
3792 : {
3793 1 : CPLError(CE_Failure, CPLE_IllegalArg,
3794 : "Unsupported gridding method \"%s\"", papszParams[0]);
3795 1 : CSLDestroy(papszParams);
3796 1 : return CE_Failure;
3797 : }
3798 :
3799 : /* -------------------------------------------------------------------- */
3800 : /* Parse algorithm parameters and assign defaults. */
3801 : /* -------------------------------------------------------------------- */
3802 382 : const char *const *papszKnownOptions = nullptr;
3803 :
3804 382 : switch (*peAlgorithm)
3805 : {
3806 245 : case GGA_InverseDistanceToAPower:
3807 : default:
3808 : {
3809 245 : *ppOptions =
3810 245 : CPLMalloc(sizeof(GDALGridInverseDistanceToAPowerOptions));
3811 :
3812 245 : GDALGridInverseDistanceToAPowerOptions *const poPowerOpts =
3813 : static_cast<GDALGridInverseDistanceToAPowerOptions *>(
3814 : *ppOptions);
3815 :
3816 245 : poPowerOpts->nSizeOfStructure = sizeof(*poPowerOpts);
3817 :
3818 245 : const char *pszValue = CSLFetchNameValue(papszParams, "power");
3819 245 : poPowerOpts->dfPower = pszValue ? CPLAtofM(pszValue) : 2.0;
3820 :
3821 245 : pszValue = CSLFetchNameValue(papszParams, "smoothing");
3822 245 : poPowerOpts->dfSmoothing = pszValue ? CPLAtofM(pszValue) : 0.0;
3823 :
3824 245 : pszValue = CSLFetchNameValue(papszParams, "radius");
3825 245 : if (pszValue)
3826 : {
3827 5 : poPowerOpts->dfRadius1 = CPLAtofM(pszValue);
3828 5 : poPowerOpts->dfRadius2 = poPowerOpts->dfRadius1;
3829 : }
3830 : else
3831 : {
3832 240 : pszValue = CSLFetchNameValue(papszParams, "radius1");
3833 240 : poPowerOpts->dfRadius1 = pszValue ? CPLAtofM(pszValue) : 0.0;
3834 :
3835 240 : pszValue = CSLFetchNameValue(papszParams, "radius2");
3836 240 : poPowerOpts->dfRadius2 = pszValue ? CPLAtofM(pszValue) : 0.0;
3837 : }
3838 :
3839 245 : pszValue = CSLFetchNameValue(papszParams, "angle");
3840 245 : poPowerOpts->dfAngle = pszValue ? CPLAtofM(pszValue) : 0.0;
3841 :
3842 245 : pszValue = CSLFetchNameValue(papszParams, "max_points");
3843 245 : poPowerOpts->nMaxPoints =
3844 245 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3845 :
3846 245 : pszValue = CSLFetchNameValue(papszParams, "min_points");
3847 245 : poPowerOpts->nMinPoints =
3848 245 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3849 :
3850 245 : pszValue = CSLFetchNameValue(papszParams, "nodata");
3851 245 : poPowerOpts->dfNoDataValue = pszValue ? CPLAtofM(pszValue) : 0.0;
3852 :
3853 : static const char *const apszKnownOptions[] = {
3854 : "power", "smoothing", "radius", "radius1", "radius2",
3855 : "angle", "max_points", "min_points", "nodata", nullptr};
3856 245 : papszKnownOptions = apszKnownOptions;
3857 :
3858 245 : break;
3859 : }
3860 22 : case GGA_InverseDistanceToAPowerNearestNeighbor:
3861 : {
3862 22 : *ppOptions = CPLMalloc(
3863 : sizeof(GDALGridInverseDistanceToAPowerNearestNeighborOptions));
3864 :
3865 : GDALGridInverseDistanceToAPowerNearestNeighborOptions
3866 22 : *const poPowerOpts = static_cast<
3867 : GDALGridInverseDistanceToAPowerNearestNeighborOptions *>(
3868 : *ppOptions);
3869 :
3870 22 : poPowerOpts->nSizeOfStructure = sizeof(*poPowerOpts);
3871 :
3872 22 : const char *pszValue = CSLFetchNameValue(papszParams, "power");
3873 22 : poPowerOpts->dfPower = pszValue ? CPLAtofM(pszValue) : 2.0;
3874 :
3875 22 : pszValue = CSLFetchNameValue(papszParams, "smoothing");
3876 22 : poPowerOpts->dfSmoothing = pszValue ? CPLAtofM(pszValue) : 0.0;
3877 :
3878 22 : pszValue = CSLFetchNameValue(papszParams, "radius");
3879 22 : poPowerOpts->dfRadius = pszValue ? CPLAtofM(pszValue) : 1.0;
3880 22 : if (!(poPowerOpts->dfRadius > 0))
3881 : {
3882 0 : CPLError(CE_Failure, CPLE_IllegalArg,
3883 : "Radius value should be strictly positive");
3884 0 : CSLDestroy(papszParams);
3885 0 : return CE_Failure;
3886 : }
3887 :
3888 22 : pszValue = CSLFetchNameValue(papszParams, "max_points");
3889 22 : poPowerOpts->nMaxPoints =
3890 22 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 12);
3891 :
3892 22 : pszValue = CSLFetchNameValue(papszParams, "min_points");
3893 22 : poPowerOpts->nMinPoints =
3894 22 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3895 :
3896 22 : pszValue = CSLFetchNameValue(papszParams, "nodata");
3897 22 : poPowerOpts->dfNoDataValue = pszValue ? CPLAtofM(pszValue) : 0.0;
3898 :
3899 : pszValue =
3900 22 : CSLFetchNameValue(papszParams, "min_points_per_quadrant");
3901 22 : poPowerOpts->nMinPointsPerQuadrant =
3902 22 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3903 :
3904 : pszValue =
3905 22 : CSLFetchNameValue(papszParams, "max_points_per_quadrant");
3906 22 : poPowerOpts->nMaxPointsPerQuadrant =
3907 22 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3908 :
3909 : static const char *const apszKnownOptions[] = {
3910 : "power",
3911 : "smoothing",
3912 : "radius",
3913 : "max_points",
3914 : "min_points",
3915 : "nodata",
3916 : "min_points_per_quadrant",
3917 : "max_points_per_quadrant",
3918 : nullptr};
3919 22 : papszKnownOptions = apszKnownOptions;
3920 :
3921 22 : break;
3922 : }
3923 26 : case GGA_MovingAverage:
3924 : {
3925 26 : *ppOptions = CPLMalloc(sizeof(GDALGridMovingAverageOptions));
3926 :
3927 26 : GDALGridMovingAverageOptions *const poAverageOpts =
3928 : static_cast<GDALGridMovingAverageOptions *>(*ppOptions);
3929 :
3930 26 : poAverageOpts->nSizeOfStructure = sizeof(*poAverageOpts);
3931 :
3932 26 : const char *pszValue = CSLFetchNameValue(papszParams, "radius");
3933 26 : if (pszValue)
3934 : {
3935 14 : poAverageOpts->dfRadius1 = CPLAtofM(pszValue);
3936 14 : poAverageOpts->dfRadius2 = poAverageOpts->dfRadius1;
3937 : }
3938 : else
3939 : {
3940 12 : pszValue = CSLFetchNameValue(papszParams, "radius1");
3941 12 : poAverageOpts->dfRadius1 = pszValue ? CPLAtofM(pszValue) : 0.0;
3942 :
3943 12 : pszValue = CSLFetchNameValue(papszParams, "radius2");
3944 12 : poAverageOpts->dfRadius2 = pszValue ? CPLAtofM(pszValue) : 0.0;
3945 : }
3946 :
3947 26 : pszValue = CSLFetchNameValue(papszParams, "angle");
3948 26 : poAverageOpts->dfAngle = pszValue ? CPLAtofM(pszValue) : 0.0;
3949 :
3950 26 : pszValue = CSLFetchNameValue(papszParams, "min_points");
3951 26 : poAverageOpts->nMinPoints =
3952 26 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3953 :
3954 26 : pszValue = CSLFetchNameValue(papszParams, "max_points");
3955 26 : poAverageOpts->nMaxPoints =
3956 26 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3957 :
3958 26 : pszValue = CSLFetchNameValue(papszParams, "nodata");
3959 26 : poAverageOpts->dfNoDataValue = pszValue ? CPLAtofM(pszValue) : 0.0;
3960 :
3961 : pszValue =
3962 26 : CSLFetchNameValue(papszParams, "min_points_per_quadrant");
3963 26 : poAverageOpts->nMinPointsPerQuadrant =
3964 26 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3965 :
3966 : pszValue =
3967 26 : CSLFetchNameValue(papszParams, "max_points_per_quadrant");
3968 26 : poAverageOpts->nMaxPointsPerQuadrant =
3969 26 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3970 :
3971 26 : if (poAverageOpts->nMinPointsPerQuadrant != 0 ||
3972 18 : poAverageOpts->nMaxPointsPerQuadrant != 0)
3973 : {
3974 9 : if (!(poAverageOpts->dfRadius1 > 0) ||
3975 9 : !(poAverageOpts->dfRadius2 > 0))
3976 : {
3977 0 : CPLError(CE_Failure, CPLE_IllegalArg,
3978 : "Radius value should be strictly positive when "
3979 : "per quadrant parameters are specified");
3980 0 : CSLDestroy(papszParams);
3981 0 : return CE_Failure;
3982 : }
3983 9 : if (poAverageOpts->dfAngle != 0)
3984 : {
3985 0 : CPLError(CE_Failure, CPLE_NotSupported,
3986 : "angle != 0 not supported when "
3987 : "per quadrant parameters are specified");
3988 0 : CSLDestroy(papszParams);
3989 0 : return CE_Failure;
3990 : }
3991 : }
3992 17 : else if (poAverageOpts->nMaxPoints > 0)
3993 : {
3994 0 : CPLError(CE_Warning, CPLE_AppDefined,
3995 : "max_points is ignored unless one of "
3996 : "min_points_per_quadrant or max_points_per_quadrant "
3997 : "is >= 1");
3998 : }
3999 :
4000 : static const char *const apszKnownOptions[] = {
4001 : "radius",
4002 : "radius1",
4003 : "radius2",
4004 : "angle",
4005 : "min_points",
4006 : "max_points",
4007 : "nodata",
4008 : "min_points_per_quadrant",
4009 : "max_points_per_quadrant",
4010 : nullptr};
4011 26 : papszKnownOptions = apszKnownOptions;
4012 :
4013 26 : break;
4014 : }
4015 19 : case GGA_NearestNeighbor:
4016 : {
4017 19 : *ppOptions = CPLMalloc(sizeof(GDALGridNearestNeighborOptions));
4018 :
4019 19 : GDALGridNearestNeighborOptions *const poNeighborOpts =
4020 : static_cast<GDALGridNearestNeighborOptions *>(*ppOptions);
4021 :
4022 19 : poNeighborOpts->nSizeOfStructure = sizeof(*poNeighborOpts);
4023 :
4024 19 : const char *pszValue = CSLFetchNameValue(papszParams, "radius");
4025 19 : if (pszValue)
4026 : {
4027 2 : poNeighborOpts->dfRadius1 = CPLAtofM(pszValue);
4028 2 : poNeighborOpts->dfRadius2 = poNeighborOpts->dfRadius1;
4029 : }
4030 : else
4031 : {
4032 17 : pszValue = CSLFetchNameValue(papszParams, "radius1");
4033 17 : poNeighborOpts->dfRadius1 = pszValue ? CPLAtofM(pszValue) : 0.0;
4034 :
4035 17 : pszValue = CSLFetchNameValue(papszParams, "radius2");
4036 17 : poNeighborOpts->dfRadius2 = pszValue ? CPLAtofM(pszValue) : 0.0;
4037 : }
4038 :
4039 19 : pszValue = CSLFetchNameValue(papszParams, "angle");
4040 19 : poNeighborOpts->dfAngle = pszValue ? CPLAtofM(pszValue) : 0.0;
4041 :
4042 19 : pszValue = CSLFetchNameValue(papszParams, "nodata");
4043 19 : poNeighborOpts->dfNoDataValue = pszValue ? CPLAtofM(pszValue) : 0.0;
4044 :
4045 : static const char *const apszKnownOptions[] = {
4046 : "radius", "radius1", "radius2", "angle", "nodata", nullptr};
4047 19 : papszKnownOptions = apszKnownOptions;
4048 :
4049 19 : break;
4050 : }
4051 63 : case GGA_MetricMinimum:
4052 : case GGA_MetricMaximum:
4053 : case GGA_MetricRange:
4054 : case GGA_MetricCount:
4055 : case GGA_MetricAverageDistance:
4056 : case GGA_MetricAverageDistancePts:
4057 : {
4058 63 : *ppOptions = CPLMalloc(sizeof(GDALGridDataMetricsOptions));
4059 :
4060 63 : GDALGridDataMetricsOptions *const poMetricsOptions =
4061 : static_cast<GDALGridDataMetricsOptions *>(*ppOptions);
4062 :
4063 63 : poMetricsOptions->nSizeOfStructure = sizeof(*poMetricsOptions);
4064 :
4065 63 : const char *pszValue = CSLFetchNameValue(papszParams, "radius");
4066 63 : if (pszValue)
4067 : {
4068 31 : poMetricsOptions->dfRadius1 = CPLAtofM(pszValue);
4069 31 : poMetricsOptions->dfRadius2 = poMetricsOptions->dfRadius1;
4070 : }
4071 : else
4072 : {
4073 32 : pszValue = CSLFetchNameValue(papszParams, "radius1");
4074 32 : poMetricsOptions->dfRadius1 =
4075 32 : pszValue ? CPLAtofM(pszValue) : 0.0;
4076 :
4077 32 : pszValue = CSLFetchNameValue(papszParams, "radius2");
4078 32 : poMetricsOptions->dfRadius2 =
4079 32 : pszValue ? CPLAtofM(pszValue) : 0.0;
4080 : }
4081 :
4082 63 : pszValue = CSLFetchNameValue(papszParams, "angle");
4083 63 : poMetricsOptions->dfAngle = pszValue ? CPLAtofM(pszValue) : 0.0;
4084 :
4085 63 : pszValue = CSLFetchNameValue(papszParams, "min_points");
4086 63 : poMetricsOptions->nMinPoints = pszValue ? atoi(pszValue) : 0;
4087 :
4088 63 : pszValue = CSLFetchNameValue(papszParams, "nodata");
4089 63 : poMetricsOptions->dfNoDataValue =
4090 63 : pszValue ? CPLAtofM(pszValue) : 0.0;
4091 :
4092 : pszValue =
4093 63 : CSLFetchNameValue(papszParams, "min_points_per_quadrant");
4094 63 : poMetricsOptions->nMinPointsPerQuadrant =
4095 63 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
4096 :
4097 : pszValue =
4098 63 : CSLFetchNameValue(papszParams, "max_points_per_quadrant");
4099 63 : poMetricsOptions->nMaxPointsPerQuadrant =
4100 63 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
4101 :
4102 63 : if (poMetricsOptions->nMinPointsPerQuadrant != 0 ||
4103 37 : poMetricsOptions->nMaxPointsPerQuadrant != 0)
4104 : {
4105 27 : if (*peAlgorithm == GGA_MetricAverageDistancePts)
4106 : {
4107 0 : CPLError(CE_Failure, CPLE_NotSupported,
4108 : "Algorithm %s not supported when "
4109 : "per quadrant parameters are specified",
4110 : szAlgNameAverageDistancePts);
4111 0 : CSLDestroy(papszParams);
4112 0 : return CE_Failure;
4113 : }
4114 27 : if (!(poMetricsOptions->dfRadius1 > 0) ||
4115 27 : !(poMetricsOptions->dfRadius2 > 0))
4116 : {
4117 0 : CPLError(CE_Failure, CPLE_IllegalArg,
4118 : "Radius value should be strictly positive when "
4119 : "per quadrant parameters are specified");
4120 0 : CSLDestroy(papszParams);
4121 0 : return CE_Failure;
4122 : }
4123 27 : if (poMetricsOptions->dfAngle != 0)
4124 : {
4125 0 : CPLError(CE_Failure, CPLE_NotSupported,
4126 : "angle != 0 not supported when "
4127 : "per quadrant parameters are specified");
4128 0 : CSLDestroy(papszParams);
4129 0 : return CE_Failure;
4130 : }
4131 : }
4132 :
4133 : static const char *const apszKnownOptions[] = {
4134 : "radius",
4135 : "radius1",
4136 : "radius2",
4137 : "angle",
4138 : "min_points",
4139 : "nodata",
4140 : "min_points_per_quadrant",
4141 : "max_points_per_quadrant",
4142 : nullptr};
4143 63 : papszKnownOptions = apszKnownOptions;
4144 :
4145 63 : break;
4146 : }
4147 7 : case GGA_Linear:
4148 : {
4149 7 : *ppOptions = CPLMalloc(sizeof(GDALGridLinearOptions));
4150 :
4151 7 : GDALGridLinearOptions *const poLinearOpts =
4152 : static_cast<GDALGridLinearOptions *>(*ppOptions);
4153 :
4154 7 : poLinearOpts->nSizeOfStructure = sizeof(*poLinearOpts);
4155 :
4156 7 : const char *pszValue = CSLFetchNameValue(papszParams, "radius");
4157 7 : poLinearOpts->dfRadius = pszValue ? CPLAtofM(pszValue) : -1.0;
4158 :
4159 7 : pszValue = CSLFetchNameValue(papszParams, "nodata");
4160 7 : poLinearOpts->dfNoDataValue = pszValue ? CPLAtofM(pszValue) : 0.0;
4161 :
4162 : static const char *const apszKnownOptions[] = {"radius", "nodata",
4163 : nullptr};
4164 7 : papszKnownOptions = apszKnownOptions;
4165 :
4166 7 : break;
4167 : }
4168 : }
4169 :
4170 382 : if (papszKnownOptions)
4171 : {
4172 1084 : for (int i = 1; papszParams[i] != nullptr; ++i)
4173 : {
4174 702 : char *pszKey = nullptr;
4175 702 : CPL_IGNORE_RET_VAL(CPLParseNameValue(papszParams[i], &pszKey));
4176 702 : if (pszKey)
4177 : {
4178 702 : bool bKnownKey = false;
4179 2992 : for (const char *const *papszKnownKeyIter = papszKnownOptions;
4180 2992 : *papszKnownKeyIter; ++papszKnownKeyIter)
4181 : {
4182 2992 : if (EQUAL(*papszKnownKeyIter, pszKey))
4183 : {
4184 702 : bKnownKey = true;
4185 702 : break;
4186 : }
4187 : }
4188 702 : if (!bKnownKey)
4189 : {
4190 0 : CPLError(CE_Warning, CPLE_AppDefined, "Option %s ignored",
4191 : pszKey);
4192 : }
4193 : }
4194 702 : CPLFree(pszKey);
4195 : }
4196 : }
4197 :
4198 382 : CSLDestroy(papszParams);
4199 382 : return CE_None;
4200 : }
|
