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