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