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 760148 : static void GDALGridGetPointBounds(const void *hFeature, CPLRectObj *pBounds)
47 : {
48 760148 : const GDALGridPoint *psPoint = static_cast<const GDALGridPoint *>(hFeature);
49 760148 : GDALGridXYArrays *psXYArrays = psPoint->psXYArrays;
50 760148 : const int i = psPoint->i;
51 760148 : const double dfX = psXYArrays->padfX[i];
52 760148 : const double dfY = psXYArrays->padfY[i];
53 760148 : pBounds->minx = dfX;
54 760148 : pBounds->miny = dfY;
55 760148 : pBounds->maxx = dfX;
56 760148 : pBounds->maxy = dfY;
57 760148 : }
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 160400 : for (GUInt32 i = 0; i < nPoints; i++)
141 : {
142 160000 : double dfRX = padfX[i] - dfXPoint;
143 160000 : double dfRY = padfY[i] - dfYPoint;
144 160000 : const double dfR2 =
145 160000 : dfRX * dfRX + dfRY * dfRY + dfSmoothing * dfSmoothing;
146 :
147 160000 : 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 160000 : 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 3364 : if (dfR2 < 0.0000000000001)
163 : {
164 0 : *pdfValue = padfZ[i];
165 0 : return CE_None;
166 : }
167 :
168 3364 : const double dfW = pow(dfR2, dfPowerDiv2);
169 3364 : const double dfInvW = 1.0 / dfW;
170 3364 : dfNominator += dfInvW * padfZ[i];
171 3364 : dfDenominator += dfInvW;
172 3364 : n++;
173 3364 : 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 1202 : 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 1202 : 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 52787 : for (int k = 0; k < nFeatureCount; k++)
280 : {
281 51581 : const int i = papsPoints[k]->i;
282 51581 : const double dfRX = padfX[i] - dfXPoint;
283 51581 : const double dfRY = padfY[i] - dfYPoint;
284 :
285 51581 : const double dfR2 = dfRX * dfRX + dfRY * dfRY;
286 : // real distance + smoothing
287 51581 : const double dfRsmoothed2 = dfR2 + dfSmoothing2;
288 51581 : 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 51581 : if (dfR2 <= dfRPower2)
296 : {
297 : oMapDistanceToZValues.insert(
298 38174 : std::make_pair(dfRsmoothed2, padfZ[i]));
299 : }
300 : }
301 : }
302 1206 : CPLFree(papsPoints);
303 :
304 1202 : double dfNominator = 0.0;
305 1202 : double dfDenominator = 0.0;
306 1202 : 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 8320 : for (std::multimap<double, double>::iterator oMapDistanceToZValuesIter =
312 1202 : oMapDistanceToZValues.begin();
313 9522 : oMapDistanceToZValuesIter != oMapDistanceToZValues.end();
314 8320 : ++oMapDistanceToZValuesIter)
315 : {
316 9121 : const double dfR2 = oMapDistanceToZValuesIter->first;
317 9121 : const double dfZ = oMapDistanceToZValuesIter->second;
318 :
319 9120 : const double dfW = pow(dfR2, dfPowerDiv2);
320 9120 : const double dfInvW = 1.0 / dfW;
321 9120 : dfNominator += dfInvW * dfZ;
322 9120 : dfDenominator += dfInvW;
323 9120 : n++;
324 9120 : 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 80401 : for (i = 0; i < nPoints; i++)
531 : {
532 80301 : const double dfRX = padfX[i] - dfXPoint;
533 80301 : const double dfRY = padfY[i] - dfYPoint;
534 80301 : 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 80301 : if (dfR2 < 0.0000000000001)
539 : {
540 401 : break;
541 : }
542 :
543 79900 : const double dfInvR2 = 1.0 / dfR2;
544 79900 : dfNominator += dfInvR2 * padfZ[i];
545 79900 : 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 2415 : 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 2415 : const GDALGridMovingAverageOptions *const poOptions =
639 : static_cast<const GDALGridMovingAverageOptions *>(poOptionsIn);
640 : // Pre-compute search ellipse parameters.
641 2415 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
642 2415 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
643 : const double dfSearchRadius =
644 2415 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
645 2415 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
646 :
647 2415 : GDALGridExtraParameters *psExtraParams =
648 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
649 2415 : const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
650 :
651 : // Compute coefficients for coordinate system rotation.
652 2415 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
653 2415 : const bool bRotated = dfAngle != 0.0;
654 :
655 2415 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
656 2415 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
657 :
658 2415 : double dfAccumulator = 0.0;
659 :
660 2415 : GUInt32 n = 0; // Used after for.
661 2415 : if (phQuadTree != nullptr)
662 : {
663 : CPLRectObj sAoi;
664 1600 : sAoi.minx = dfXPoint - dfSearchRadius;
665 1600 : sAoi.miny = dfYPoint - dfSearchRadius;
666 1600 : sAoi.maxx = dfXPoint + dfSearchRadius;
667 1600 : sAoi.maxy = dfYPoint + dfSearchRadius;
668 1600 : int nFeatureCount = 0;
669 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
670 1600 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
671 1600 : if (nFeatureCount != 0)
672 : {
673 41096 : for (int k = 0; k < nFeatureCount; k++)
674 : {
675 39496 : const int i = papsPoints[k]->i;
676 39496 : const double dfRX = padfX[i] - dfXPoint;
677 39496 : const double dfRY = padfY[i] - dfYPoint;
678 :
679 39496 : if (dfRadius2Square * dfRX * dfRX +
680 39496 : dfRadius1Square * dfRY * dfRY <=
681 : dfR12Square)
682 : {
683 32288 : dfAccumulator += padfZ[i];
684 32288 : n++;
685 : }
686 : }
687 : }
688 1600 : CPLFree(papsPoints);
689 : }
690 : else
691 : {
692 320905 : for (GUInt32 i = 0; i < nPoints; i++)
693 : {
694 320090 : double dfRX = padfX[i] - dfXPoint;
695 320090 : double dfRY = padfY[i] - dfYPoint;
696 :
697 320090 : if (bRotated)
698 : {
699 160000 : const double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
700 160000 : const double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
701 :
702 160000 : dfRX = dfRXRotated;
703 160000 : dfRY = dfRYRotated;
704 : }
705 :
706 : // Is this point located inside the search ellipse?
707 320090 : if (dfRadius2Square * dfRX * dfRX + dfRadius1Square * dfRY * dfRY <=
708 : dfR12Square)
709 : {
710 169310 : dfAccumulator += padfZ[i];
711 169310 : n++;
712 : }
713 : }
714 : }
715 :
716 2415 : if (n < poOptions->nMinPoints || n == 0)
717 : {
718 152 : *pdfValue = poOptions->dfNoDataValue;
719 : }
720 : else
721 : {
722 2263 : *pdfValue = dfAccumulator / n;
723 : }
724 :
725 2415 : 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 44855 : 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 44855 : const GDALGridNearestNeighborOptions *const poOptions =
889 : static_cast<const GDALGridNearestNeighborOptions *>(poOptionsIn);
890 : // Pre-compute search ellipse parameters.
891 44855 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
892 44855 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
893 44855 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
894 44855 : GDALGridExtraParameters *psExtraParams =
895 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
896 44855 : CPLQuadTree *hQuadTree = psExtraParams->hQuadTree;
897 :
898 : // Compute coefficients for coordinate system rotation.
899 44855 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
900 44855 : const bool bRotated = dfAngle != 0.0;
901 44855 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
902 44855 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
903 :
904 : // If the nearest point will not be found, its value remains as NODATA.
905 44855 : double dfNearestValue = poOptions->dfNoDataValue;
906 44855 : GUInt32 i = 0;
907 :
908 44855 : double dfSearchRadius = psExtraParams->dfInitialSearchRadius;
909 44855 : if (hQuadTree != nullptr)
910 : {
911 14369 : if (poOptions->dfRadius1 > 0 || poOptions->dfRadius2 > 0)
912 3199 : dfSearchRadius =
913 3207 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
914 : CPLRectObj sAoi;
915 14964 : while (dfSearchRadius > 0)
916 : {
917 14964 : sAoi.minx = dfXPoint - dfSearchRadius;
918 14964 : sAoi.miny = dfYPoint - dfSearchRadius;
919 14964 : sAoi.maxx = dfXPoint + dfSearchRadius;
920 14964 : sAoi.maxy = dfYPoint + dfSearchRadius;
921 14964 : int nFeatureCount = 0;
922 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
923 14964 : CPLQuadTreeSearch(hQuadTree, &sAoi, &nFeatureCount));
924 14912 : if (nFeatureCount != 0)
925 : {
926 : // Nearest distance will be initialized with the distance to the
927 : // first point in array.
928 14316 : double dfNearestRSquare = std::numeric_limits<double>::max();
929 104874 : for (int k = 0; k < nFeatureCount; k++)
930 : {
931 90558 : const int idx = papsPoints[k]->i;
932 90558 : const double dfRX = padfX[idx] - dfXPoint;
933 90558 : const double dfRY = padfY[idx] - dfYPoint;
934 :
935 90558 : const double dfR2 = dfRX * dfRX + dfRY * dfRY;
936 90558 : if (dfR2 <= dfNearestRSquare)
937 : {
938 32282 : dfNearestRSquare = dfR2;
939 32282 : dfNearestValue = padfZ[idx];
940 : }
941 : }
942 :
943 14316 : CPLFree(papsPoints);
944 14347 : break;
945 : }
946 :
947 596 : CPLFree(papsPoints);
948 603 : if (poOptions->dfRadius1 > 0 || poOptions->dfRadius2 > 0)
949 : break;
950 603 : 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 214614000 : while (i < nPoints)
961 : {
962 214584000 : double dfRX = padfX[i] - dfXPoint;
963 214584000 : double dfRY = padfY[i] - dfYPoint;
964 :
965 214584000 : 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 214584000 : const double dfRXSquare = dfRX * dfRX;
976 214584000 : const double dfRYSquare = dfRY * dfRY;
977 214584000 : if (dfRadius2Square * dfRXSquare + dfRadius1Square * dfRYSquare <=
978 : dfR12Square)
979 : {
980 193230000 : const double dfR2 = dfRXSquare + dfRYSquare;
981 193230000 : if (dfR2 <= dfNearestRSquare)
982 : {
983 15860100 : dfNearestRSquare = dfR2;
984 15860100 : dfNearestValue = padfZ[i];
985 : }
986 : }
987 :
988 214584000 : i++;
989 : }
990 : }
991 :
992 44833 : *pdfValue = dfNearestValue;
993 :
994 44833 : 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 2400 : 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 2400 : const GDALGridDataMetricsOptions *const poOptions =
1043 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
1044 :
1045 : // Pre-compute search ellipse parameters.
1046 2400 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
1047 2400 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
1048 : const double dfSearchRadius =
1049 2400 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
1050 2400 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
1051 :
1052 2400 : GDALGridExtraParameters *psExtraParams =
1053 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
1054 2400 : CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
1055 :
1056 : // Compute coefficients for coordinate system rotation.
1057 2400 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
1058 2400 : const bool bRotated = dfAngle != 0.0;
1059 2400 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
1060 2400 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
1061 :
1062 2400 : double dfMinimumValue = std::numeric_limits<double>::max();
1063 2400 : GUInt32 n = 0;
1064 2400 : if (phQuadTree != nullptr)
1065 : {
1066 : CPLRectObj sAoi;
1067 1600 : sAoi.minx = dfXPoint - dfSearchRadius;
1068 1600 : sAoi.miny = dfYPoint - dfSearchRadius;
1069 1600 : sAoi.maxx = dfXPoint + dfSearchRadius;
1070 1600 : sAoi.maxy = dfYPoint + dfSearchRadius;
1071 1600 : int nFeatureCount = 0;
1072 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
1073 1600 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
1074 1600 : if (nFeatureCount != 0)
1075 : {
1076 35168 : for (int k = 0; k < nFeatureCount; k++)
1077 : {
1078 33568 : const int i = papsPoints[k]->i;
1079 33568 : const double dfRX = padfX[i] - dfXPoint;
1080 33568 : const double dfRY = padfY[i] - dfYPoint;
1081 :
1082 33568 : if (dfRadius2Square * dfRX * dfRX +
1083 33568 : dfRadius1Square * dfRY * dfRY <=
1084 : dfR12Square)
1085 : {
1086 21192 : if (dfMinimumValue > padfZ[i])
1087 3280 : dfMinimumValue = padfZ[i];
1088 21192 : n++;
1089 : }
1090 : }
1091 : }
1092 1600 : CPLFree(papsPoints);
1093 : }
1094 : else
1095 : {
1096 800 : GUInt32 i = 0;
1097 320800 : while (i < nPoints)
1098 : {
1099 320000 : double dfRX = padfX[i] - dfXPoint;
1100 320000 : double dfRY = padfY[i] - dfYPoint;
1101 :
1102 320000 : if (bRotated)
1103 : {
1104 160000 : const double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
1105 160000 : const double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
1106 :
1107 160000 : dfRX = dfRXRotated;
1108 160000 : dfRY = dfRYRotated;
1109 : }
1110 :
1111 : // Is this point located inside the search ellipse?
1112 320000 : if (dfRadius2Square * dfRX * dfRX + dfRadius1Square * dfRY * dfRY <=
1113 : dfR12Square)
1114 : {
1115 171394 : if (dfMinimumValue > padfZ[i])
1116 2873 : dfMinimumValue = padfZ[i];
1117 171394 : n++;
1118 : }
1119 :
1120 320000 : i++;
1121 : }
1122 : }
1123 :
1124 2400 : if (n < poOptions->nMinPoints || n == 0)
1125 : {
1126 0 : *pdfValue = poOptions->dfNoDataValue;
1127 : }
1128 : else
1129 : {
1130 2400 : *pdfValue = dfMinimumValue;
1131 : }
1132 :
1133 2400 : 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 2400 : 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 2400 : const GDALGridDataMetricsOptions *const poOptions =
1337 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
1338 :
1339 : // Pre-compute search ellipse parameters.
1340 2400 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
1341 2400 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
1342 : const double dfSearchRadius =
1343 2400 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
1344 2400 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
1345 :
1346 2400 : GDALGridExtraParameters *psExtraParams =
1347 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
1348 2400 : CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
1349 :
1350 : // Compute coefficients for coordinate system rotation.
1351 2400 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
1352 2400 : const bool bRotated = dfAngle != 0.0;
1353 2400 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
1354 2400 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
1355 :
1356 2400 : double dfMaximumValue = -std::numeric_limits<double>::max();
1357 2400 : GUInt32 n = 0;
1358 2400 : if (phQuadTree != nullptr)
1359 : {
1360 : CPLRectObj sAoi;
1361 1200 : sAoi.minx = dfXPoint - dfSearchRadius;
1362 1200 : sAoi.miny = dfYPoint - dfSearchRadius;
1363 1200 : sAoi.maxx = dfXPoint + dfSearchRadius;
1364 1200 : sAoi.maxy = dfYPoint + dfSearchRadius;
1365 1200 : int nFeatureCount = 0;
1366 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
1367 1200 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
1368 1200 : if (nFeatureCount != 0)
1369 : {
1370 37332 : for (int k = 0; k < nFeatureCount; k++)
1371 : {
1372 36132 : const int i = papsPoints[k]->i;
1373 36132 : const double dfRX = padfX[i] - dfXPoint;
1374 36132 : const double dfRY = padfY[i] - dfYPoint;
1375 :
1376 36132 : if (dfRadius2Square * dfRX * dfRX +
1377 36132 : dfRadius1Square * dfRY * dfRY <=
1378 : dfR12Square)
1379 : {
1380 23756 : if (dfMaximumValue < padfZ[i])
1381 3338 : dfMaximumValue = padfZ[i];
1382 23756 : n++;
1383 : }
1384 : }
1385 : }
1386 1200 : CPLFree(papsPoints);
1387 : }
1388 : else
1389 : {
1390 1200 : GUInt32 i = 0;
1391 481200 : while (i < nPoints)
1392 : {
1393 480000 : double dfRX = padfX[i] - dfXPoint;
1394 480000 : double dfRY = padfY[i] - dfYPoint;
1395 :
1396 480000 : if (bRotated)
1397 : {
1398 320000 : const double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
1399 320000 : const double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
1400 :
1401 320000 : dfRX = dfRXRotated;
1402 320000 : dfRY = dfRYRotated;
1403 : }
1404 :
1405 : // Is this point located inside the search ellipse?
1406 480000 : if (dfRadius2Square * dfRX * dfRX + dfRadius1Square * dfRY * dfRY <=
1407 : dfR12Square)
1408 : {
1409 160800 : if (dfMaximumValue < padfZ[i])
1410 4800 : dfMaximumValue = padfZ[i];
1411 160800 : n++;
1412 : }
1413 :
1414 480000 : i++;
1415 : }
1416 : }
1417 :
1418 2400 : if (n < poOptions->nMinPoints || n == 0)
1419 : {
1420 0 : *pdfValue = poOptions->dfNoDataValue;
1421 : }
1422 : else
1423 : {
1424 2400 : *pdfValue = dfMaximumValue;
1425 : }
1426 :
1427 2400 : 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 1200 : 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 1200 : const GDALGridDataMetricsOptions *const poOptions =
1494 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
1495 : // Pre-compute search ellipse parameters.
1496 1200 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
1497 1200 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
1498 : const double dfSearchRadius =
1499 1200 : 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 800 : if (nFeatureCount != 0)
1526 : {
1527 7528 : for (int k = 0; k < nFeatureCount; k++)
1528 : {
1529 6728 : const int i = papsPoints[k]->i;
1530 6728 : const double dfRX = padfX[i] - dfXPoint;
1531 6728 : const double dfRY = padfY[i] - dfYPoint;
1532 :
1533 6728 : if (dfRadius2Square * dfRX * dfRX +
1534 6728 : dfRadius1Square * dfRY * dfRY <=
1535 : dfR12Square)
1536 : {
1537 6728 : if (dfMinimumValue > padfZ[i])
1538 1948 : dfMinimumValue = padfZ[i];
1539 6728 : if (dfMaximumValue < padfZ[i])
1540 1836 : dfMaximumValue = padfZ[i];
1541 6728 : n++;
1542 : }
1543 : }
1544 : }
1545 800 : CPLFree(papsPoints);
1546 : }
1547 : else
1548 : {
1549 400 : GUInt32 i = 0;
1550 160400 : while (i < nPoints)
1551 : {
1552 160000 : double dfRX = padfX[i] - dfXPoint;
1553 160000 : double dfRY = padfY[i] - dfYPoint;
1554 :
1555 160000 : 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 160000 : if (dfRadius2Square * dfRX * dfRX + dfRadius1Square * dfRY * dfRY <=
1566 : dfR12Square)
1567 : {
1568 160000 : if (dfMinimumValue > padfZ[i])
1569 1600 : dfMinimumValue = padfZ[i];
1570 160000 : if (dfMaximumValue < padfZ[i])
1571 4000 : dfMaximumValue = padfZ[i];
1572 160000 : n++;
1573 : }
1574 :
1575 160000 : i++;
1576 : }
1577 : }
1578 :
1579 1200 : if (n < poOptions->nMinPoints || n == 0)
1580 : {
1581 152 : *pdfValue = poOptions->dfNoDataValue;
1582 : }
1583 : else
1584 : {
1585 1048 : *pdfValue = dfMaximumValue - dfMinimumValue;
1586 : }
1587 :
1588 1200 : 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 1998 : 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 1998 : const GDALGridDataMetricsOptions *const poOptions =
1765 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
1766 :
1767 : // Pre-compute search ellipse parameters.
1768 1998 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
1769 1998 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
1770 : const double dfSearchRadius =
1771 1998 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
1772 2000 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
1773 :
1774 2000 : GDALGridExtraParameters *psExtraParams =
1775 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
1776 2000 : CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
1777 :
1778 : // Compute coefficients for coordinate system rotation.
1779 2000 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
1780 2000 : const bool bRotated = dfAngle != 0.0;
1781 2000 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
1782 2000 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
1783 :
1784 2000 : GUInt32 n = 0;
1785 2000 : 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 2000 : if (nFeatureCount != 0)
1796 : {
1797 79983 : for (int k = 0; k < nFeatureCount; k++)
1798 : {
1799 77983 : const int i = papsPoints[k]->i;
1800 77983 : const double dfRX = padfX[i] - dfXPoint;
1801 77983 : const double dfRY = padfY[i] - dfYPoint;
1802 :
1803 77983 : if (dfRadius2Square * dfRX * dfRX +
1804 77983 : dfRadius1Square * dfRY * dfRY <=
1805 : dfR12Square)
1806 : {
1807 55149 : n++;
1808 : }
1809 : }
1810 : }
1811 2000 : CPLFree(papsPoints);
1812 : }
1813 : else
1814 : {
1815 1 : GUInt32 i = 0;
1816 1 : 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 2000 : if (n < poOptions->nMinPoints)
1842 : {
1843 0 : *pdfValue = poOptions->dfNoDataValue;
1844 : }
1845 : else
1846 : {
1847 2000 : *pdfValue = static_cast<double>(n);
1848 : }
1849 :
1850 2000 : 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 800 : sAoi.minx = dfXPoint - dfSearchRadius;
2050 800 : sAoi.miny = dfYPoint - dfSearchRadius;
2051 800 : sAoi.maxx = dfXPoint + dfSearchRadius;
2052 800 : sAoi.maxy = dfYPoint + dfSearchRadius;
2053 800 : int nFeatureCount = 0;
2054 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
2055 800 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
2056 800 : if (nFeatureCount != 0)
2057 : {
2058 18327 : for (int k = 0; k < nFeatureCount; k++)
2059 : {
2060 17527 : const int i = papsPoints[k]->i;
2061 17527 : const double dfRX = padfX[i] - dfXPoint;
2062 17527 : const double dfRY = padfY[i] - dfYPoint;
2063 :
2064 17527 : if (dfRadius2Square * dfRX * dfRX +
2065 17527 : dfRadius1Square * dfRY * dfRY <=
2066 : dfR12Square)
2067 : {
2068 14993 : dfAccumulator += sqrt(dfRX * dfRX + dfRY * dfRY);
2069 14993 : n++;
2070 : }
2071 : }
2072 : }
2073 800 : CPLFree(papsPoints);
2074 : }
2075 : else
2076 : {
2077 399 : GUInt32 i = 0;
2078 :
2079 106952 : while (i < nPoints)
2080 : {
2081 106553 : double dfRX = padfX[i] - dfXPoint;
2082 106553 : double dfRY = padfY[i] - dfYPoint;
2083 :
2084 106553 : 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 106553 : if (dfRadius2Square * dfRX * dfRX + dfRadius1Square * dfRY * dfRY <=
2095 : dfR12Square)
2096 : {
2097 109712 : dfAccumulator += sqrt(dfRX * dfRX + dfRY * dfRY);
2098 109712 : n++;
2099 : }
2100 :
2101 106553 : i++;
2102 : }
2103 : }
2104 :
2105 1199 : if (n < poOptions->nMinPoints || n == 0)
2106 : {
2107 0 : *pdfValue = poOptions->dfNoDataValue;
2108 : }
2109 : else
2110 : {
2111 1199 : *pdfValue = dfAccumulator / n;
2112 : }
2113 :
2114 1199 : 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 1200 : 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 1200 : const GDALGridDataMetricsOptions *const poOptions =
2292 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
2293 : // Pre-compute search ellipse parameters.
2294 1200 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
2295 1200 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
2296 : const double dfSearchRadius =
2297 1200 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
2298 1200 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
2299 :
2300 1200 : GDALGridExtraParameters *psExtraParams =
2301 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
2302 1200 : CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
2303 :
2304 : // Compute coefficients for coordinate system rotation.
2305 1200 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
2306 1200 : const bool bRotated = dfAngle != 0.0;
2307 1200 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
2308 1200 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
2309 :
2310 1200 : double dfAccumulator = 0.0;
2311 1200 : GUInt32 n = 0;
2312 1200 : if (phQuadTree != nullptr)
2313 : {
2314 : CPLRectObj sAoi;
2315 800 : sAoi.minx = dfXPoint - dfSearchRadius;
2316 800 : sAoi.miny = dfYPoint - dfSearchRadius;
2317 800 : sAoi.maxx = dfXPoint + dfSearchRadius;
2318 800 : sAoi.maxy = dfYPoint + dfSearchRadius;
2319 800 : int nFeatureCount = 0;
2320 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
2321 800 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
2322 800 : if (nFeatureCount != 0)
2323 : {
2324 17672 : for (int k = 0; k < nFeatureCount - 1; k++)
2325 : {
2326 16872 : const int i = papsPoints[k]->i;
2327 16872 : const double dfRX1 = padfX[i] - dfXPoint;
2328 16872 : const double dfRY1 = padfY[i] - dfYPoint;
2329 :
2330 16872 : if (dfRadius2Square * dfRX1 * dfRX1 +
2331 16872 : dfRadius1Square * dfRY1 * dfRY1 <=
2332 : dfR12Square)
2333 : {
2334 193952 : 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 179318 : const int ji = papsPoints[j]->i;
2339 179318 : double dfRX2 = padfX[ji] - dfXPoint;
2340 179318 : double dfRY2 = padfY[ji] - dfYPoint;
2341 :
2342 179318 : if (dfRadius2Square * dfRX2 * dfRX2 +
2343 179318 : dfRadius1Square * dfRY2 * dfRY2 <=
2344 : dfR12Square)
2345 : {
2346 153666 : const double dfRX = padfX[ji] - padfX[i];
2347 153666 : const double dfRY = padfY[ji] - padfY[i];
2348 :
2349 153666 : dfAccumulator += sqrt(dfRX * dfRX + dfRY * dfRY);
2350 153666 : n++;
2351 : }
2352 : }
2353 : }
2354 : }
2355 : }
2356 800 : CPLFree(papsPoints);
2357 : }
2358 : else
2359 : {
2360 400 : GUInt32 i = 0;
2361 160000 : while (i < nPoints - 1)
2362 : {
2363 159600 : double dfRX1 = padfX[i] - dfXPoint;
2364 159600 : double dfRY1 = padfY[i] - dfYPoint;
2365 :
2366 159600 : if (bRotated)
2367 : {
2368 159600 : const double dfRXRotated = dfRX1 * dfCoeff1 + dfRY1 * dfCoeff2;
2369 159600 : const double dfRYRotated = dfRY1 * dfCoeff1 - dfRX1 * dfCoeff2;
2370 :
2371 159600 : dfRX1 = dfRXRotated;
2372 159600 : dfRY1 = dfRYRotated;
2373 : }
2374 :
2375 : // Is this point located inside the search ellipse?
2376 159600 : if (dfRadius2Square * dfRX1 * dfRX1 +
2377 159600 : 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 409618 : for (GUInt32 j = i + 1; j < nPoints; j++)
2383 : {
2384 407020 : double dfRX2 = padfX[j] - dfXPoint;
2385 407020 : double dfRY2 = padfY[j] - dfYPoint;
2386 :
2387 407020 : if (bRotated)
2388 : {
2389 519099 : const double dfRXRotated =
2390 519099 : dfRX2 * dfCoeff1 + dfRY2 * dfCoeff2;
2391 519099 : const double dfRYRotated =
2392 519099 : dfRY2 * dfCoeff1 - dfRX2 * dfCoeff2;
2393 :
2394 519099 : dfRX2 = dfRXRotated;
2395 519099 : dfRY2 = dfRYRotated;
2396 : }
2397 :
2398 407020 : if (dfRadius2Square * dfRX2 * dfRX2 +
2399 407020 : dfRadius1Square * dfRY2 * dfRY2 <=
2400 : dfR12Square)
2401 : {
2402 7342 : const double dfRX = padfX[j] - padfX[i];
2403 7342 : const double dfRY = padfY[j] - padfY[i];
2404 :
2405 7342 : dfAccumulator += sqrt(dfRX * dfRX + dfRY * dfRY);
2406 7342 : n++;
2407 : }
2408 : }
2409 : }
2410 :
2411 159600 : i++;
2412 : }
2413 : }
2414 :
2415 : // Search for the first point within the search ellipse.
2416 1200 : if (n < poOptions->nMinPoints || n == 0)
2417 : {
2418 0 : *pdfValue = poOptions->dfNoDataValue;
2419 : }
2420 : else
2421 : {
2422 1200 : *pdfValue = dfAccumulator / n;
2423 : }
2424 :
2425 1200 : 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 25316 : 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 25316 : GDALGridExtraParameters *psExtraParams =
2465 : static_cast<GDALGridExtraParameters *>(hExtraParams);
2466 25316 : GDALTriangulation *psTriangulation = psExtraParams->psTriangulation;
2467 :
2468 25316 : int nOutputFacetIdx = -1;
2469 25316 : const bool bRet = CPL_TO_BOOL(GDALTriangulationFindFacetDirected(
2470 : psTriangulation, psExtraParams->nInitialFacetIdx, dfXPoint, dfYPoint,
2471 : &nOutputFacetIdx));
2472 :
2473 25323 : if (bRet)
2474 : {
2475 14161 : CPLAssert(nOutputFacetIdx >= 0);
2476 : // Reuse output facet idx as next initial index since we proceed line by
2477 : // line.
2478 14161 : psExtraParams->nInitialFacetIdx = nOutputFacetIdx;
2479 :
2480 14161 : double lambda1 = 0.0;
2481 14161 : double lambda2 = 0.0;
2482 14161 : double lambda3 = 0.0;
2483 14161 : GDALTriangulationComputeBarycentricCoordinates(
2484 : psTriangulation, nOutputFacetIdx, dfXPoint, dfYPoint, &lambda1,
2485 : &lambda2, &lambda3);
2486 14161 : const int i1 =
2487 14161 : psTriangulation->pasFacets[nOutputFacetIdx].anVertexIdx[0];
2488 14161 : const int i2 =
2489 14161 : psTriangulation->pasFacets[nOutputFacetIdx].anVertexIdx[1];
2490 14161 : const int i3 =
2491 14161 : psTriangulation->pasFacets[nOutputFacetIdx].anVertexIdx[2];
2492 14161 : *pdfValue =
2493 14161 : lambda1 * padfZ[i1] + lambda2 * padfZ[i2] + lambda3 * padfZ[i3];
2494 : }
2495 : else
2496 : {
2497 11162 : 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 11162 : const GDALGridLinearOptions *const poOptions =
2505 : static_cast<const GDALGridLinearOptions *>(poOptionsIn);
2506 11162 : const double dfRadius = poOptions->dfRadius;
2507 11162 : if (dfRadius == 0.0)
2508 : {
2509 0 : *pdfValue = poOptions->dfNoDataValue;
2510 : }
2511 : else
2512 : {
2513 : GDALGridNearestNeighborOptions sNeighbourOptions;
2514 11162 : sNeighbourOptions.nSizeOfStructure = sizeof(sNeighbourOptions);
2515 11162 : sNeighbourOptions.dfRadius1 = dfRadius < 0.0 ? 0.0 : dfRadius;
2516 11162 : sNeighbourOptions.dfRadius2 = dfRadius < 0.0 ? 0.0 : dfRadius;
2517 11162 : sNeighbourOptions.dfAngle = 0.0;
2518 11162 : sNeighbourOptions.dfNoDataValue = poOptions->dfNoDataValue;
2519 11162 : return GDALGridNearestNeighbor(&sNeighbourOptions, nPoints, padfX,
2520 : padfY, padfZ, dfXPoint, dfYPoint,
2521 10349 : pdfValue, hExtraParams);
2522 : }
2523 : }
2524 :
2525 14161 : 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 1544 : static int GDALGridProgressMultiThread(GDALGridJob *psJob)
2571 : {
2572 1544 : 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 426 : 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 426 : const GUInt32 nYStart = psJob->nYStart;
2626 426 : const GUInt32 nYStep = psJob->nYStep;
2627 426 : GByte *pabyData = psJob->pabyData;
2628 :
2629 426 : const GUInt32 nYSize = psJob->nYSize;
2630 426 : const double dfXMin = psJob->dfXMin;
2631 426 : const double dfYMin = psJob->dfYMin;
2632 426 : const double dfDeltaX = psJob->dfDeltaX;
2633 426 : const double dfDeltaY = psJob->dfDeltaY;
2634 426 : const GUInt32 nPoints = psJob->nPoints;
2635 426 : const double *padfX = psJob->padfX;
2636 426 : const double *padfY = psJob->padfY;
2637 426 : const double *padfZ = psJob->padfZ;
2638 426 : const void *poOptions = psJob->poOptions;
2639 426 : GDALGridFunction pfnGDALGridMethod = psJob->pfnGDALGridMethod;
2640 : // Have a local copy of sExtraParameters since we want to modify
2641 : // nInitialFacetIdx.
2642 426 : GDALGridExtraParameters sExtraParameters = *psJob->psExtraParameters;
2643 426 : const GDALDataType eType = psJob->eType;
2644 :
2645 426 : const int nDataTypeSize = GDALGetDataTypeSizeBytes(eType);
2646 426 : const int nLineSpace = nXSize * nDataTypeSize;
2647 :
2648 2010 : for (GUInt32 nYPoint = nYStart; nYPoint < nYSize; nYPoint += nYStep)
2649 : {
2650 1584 : const double dfYPoint = dfYMin + (nYPoint + 0.5) * dfDeltaY;
2651 :
2652 77367 : for (GUInt32 nXPoint = 0; nXPoint < nXSize; nXPoint++)
2653 : {
2654 75781 : const double dfXPoint = dfXMin + (nXPoint + 0.5) * dfDeltaX;
2655 :
2656 151564 : if ((*pfnGDALGridMethod)(poOptions, nPoints, padfX, padfY, padfZ,
2657 75781 : dfXPoint, dfYPoint, padfValues + nXPoint,
2658 75783 : &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 1586 : GDALCopyWords(padfValues, GDT_Float64, sizeof(double),
2672 1586 : pabyData + nYPoint * nLineSpace, eType, nDataTypeSize,
2673 : nXSize);
2674 :
2675 1583 : if (*psJob->pbStop || (pfnProgress != nullptr && pfnProgress(psJob)))
2676 0 : break;
2677 : }
2678 :
2679 426 : 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(CPLGetConfigOption("GDAL_USE_SSE", "YES"))
2846 : #if !defined(USE_NEON_OPTIMIZATIONS)
2847 20 : && CPLHaveRuntimeSSE()
2848 : #endif
2849 : )
2850 : {
2851 : pafXAligned = static_cast<float *>(
2852 3 : VSI_MALLOC_ALIGNED_AUTO_VERBOSE(sizeof(float) *
2853 : nPoints));
2854 : pafYAligned = static_cast<float *>(
2855 3 : VSI_MALLOC_ALIGNED_AUTO_VERBOSE(sizeof(float) *
2856 : nPoints));
2857 : pafZAligned = static_cast<float *>(
2858 3 : VSI_MALLOC_ALIGNED_AUTO_VERBOSE(sizeof(float) *
2859 : nPoints));
2860 3 : if (pafXAligned != nullptr && pafYAligned != nullptr &&
2861 : pafZAligned != nullptr)
2862 : {
2863 3 : CPLDebug("GDAL_GRID",
2864 : "Using SSE optimized version");
2865 3 : pfnGDALGridMethod =
2866 : GDALGridInverseDistanceToAPower2NoSmoothingNoSearchSSE;
2867 405 : for (GUInt32 i = 0; i < nPoints; i++)
2868 : {
2869 402 : pafXAligned[i] = static_cast<float>(padfX[i]);
2870 402 : pafYAligned[i] = static_cast<float>(padfY[i]);
2871 402 : pafZAligned[i] = static_cast<float>(padfZ[i]);
2872 3 : }
2873 : }
2874 : else
2875 : {
2876 0 : VSIFree(pafXAligned);
2877 0 : VSIFree(pafYAligned);
2878 0 : VSIFree(pafZAligned);
2879 0 : pafXAligned = nullptr;
2880 0 : pafYAligned = nullptr;
2881 0 : pafZAligned = nullptr;
2882 : }
2883 : }
2884 : #endif // HAVE_SSE_AT_COMPILE_TIME
2885 14 : }
2886 : }
2887 : else
2888 : {
2889 1 : pfnGDALGridMethod = GDALGridInverseDistanceToAPower;
2890 : }
2891 15 : break;
2892 : }
2893 11 : case GGA_InverseDistanceToAPowerNearestNeighbor:
2894 : {
2895 11 : const auto poOptionsOld = static_cast<
2896 : const GDALGridInverseDistanceToAPowerNearestNeighborOptions *>(
2897 : poOptions);
2898 11 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
2899 : {
2900 0 : CPLError(CE_Failure, CPLE_AppDefined,
2901 : "Wrong value of nSizeOfStructure member");
2902 0 : return nullptr;
2903 : }
2904 11 : poOptionsNew = CPLMalloc(
2905 : sizeof(GDALGridInverseDistanceToAPowerNearestNeighborOptions));
2906 11 : memcpy(
2907 : poOptionsNew, poOptions,
2908 : sizeof(GDALGridInverseDistanceToAPowerNearestNeighborOptions));
2909 :
2910 11 : if (poOptionsOld->nMinPointsPerQuadrant != 0 ||
2911 5 : poOptionsOld->nMaxPointsPerQuadrant != 0)
2912 : {
2913 6 : pfnGDALGridMethod =
2914 : GDALGridInverseDistanceToAPowerNearestNeighborPerQuadrant;
2915 : }
2916 : else
2917 : {
2918 5 : pfnGDALGridMethod =
2919 : GDALGridInverseDistanceToAPowerNearestNeighbor;
2920 : }
2921 11 : bCreateQuadTree = true;
2922 11 : break;
2923 : }
2924 13 : case GGA_MovingAverage:
2925 : {
2926 13 : const auto poOptionsOld =
2927 : static_cast<const GDALGridMovingAverageOptions *>(poOptions);
2928 13 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
2929 : {
2930 0 : CPLError(CE_Failure, CPLE_AppDefined,
2931 : "Wrong value of nSizeOfStructure member");
2932 0 : return nullptr;
2933 : }
2934 13 : poOptionsNew = CPLMalloc(sizeof(GDALGridMovingAverageOptions));
2935 13 : memcpy(poOptionsNew, poOptions,
2936 : sizeof(GDALGridMovingAverageOptions));
2937 :
2938 13 : if (poOptionsOld->nMinPointsPerQuadrant != 0 ||
2939 7 : poOptionsOld->nMaxPointsPerQuadrant != 0)
2940 : {
2941 6 : pfnGDALGridMethod = GDALGridMovingAveragePerQuadrant;
2942 6 : bCreateQuadTree = true;
2943 : }
2944 : else
2945 : {
2946 7 : pfnGDALGridMethod = GDALGridMovingAverage;
2947 7 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
2948 12 : poOptionsOld->dfAngle == 0.0 &&
2949 5 : (poOptionsOld->dfRadius1 > 0.0 ||
2950 1 : poOptionsOld->dfRadius2 > 0.0));
2951 : }
2952 13 : break;
2953 : }
2954 16 : case GGA_NearestNeighbor:
2955 : {
2956 16 : const auto poOptionsOld =
2957 : static_cast<const GDALGridNearestNeighborOptions *>(poOptions);
2958 16 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
2959 : {
2960 0 : CPLError(CE_Failure, CPLE_AppDefined,
2961 : "Wrong value of nSizeOfStructure member");
2962 0 : return nullptr;
2963 : }
2964 16 : poOptionsNew = CPLMalloc(sizeof(GDALGridNearestNeighborOptions));
2965 16 : memcpy(poOptionsNew, poOptions,
2966 : sizeof(GDALGridNearestNeighborOptions));
2967 :
2968 16 : pfnGDALGridMethod = GDALGridNearestNeighbor;
2969 29 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
2970 29 : poOptionsOld->dfAngle == 0.0 &&
2971 13 : (poOptionsOld->dfRadius1 > 0.0 ||
2972 5 : poOptionsOld->dfRadius2 > 0.0));
2973 16 : break;
2974 : }
2975 11 : case GGA_MetricMinimum:
2976 : {
2977 11 : const auto poOptionsOld =
2978 : static_cast<const GDALGridDataMetricsOptions *>(poOptions);
2979 11 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
2980 : {
2981 0 : CPLError(CE_Failure, CPLE_AppDefined,
2982 : "Wrong value of nSizeOfStructure member");
2983 0 : return nullptr;
2984 : }
2985 11 : poOptionsNew = CPLMalloc(sizeof(GDALGridDataMetricsOptions));
2986 11 : memcpy(poOptionsNew, poOptions, sizeof(GDALGridDataMetricsOptions));
2987 :
2988 11 : if (poOptionsOld->nMinPointsPerQuadrant != 0 ||
2989 6 : poOptionsOld->nMaxPointsPerQuadrant != 0)
2990 : {
2991 5 : pfnGDALGridMethod = GDALGridDataMetricMinimumPerQuadrant;
2992 5 : bCreateQuadTree = true;
2993 : }
2994 : else
2995 : {
2996 6 : pfnGDALGridMethod = GDALGridDataMetricMinimum;
2997 6 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
2998 11 : poOptionsOld->dfAngle == 0.0 &&
2999 5 : (poOptionsOld->dfRadius1 > 0.0 ||
3000 1 : poOptionsOld->dfRadius2 > 0.0));
3001 : }
3002 11 : break;
3003 : }
3004 11 : case GGA_MetricMaximum:
3005 : {
3006 11 : const auto poOptionsOld =
3007 : static_cast<const GDALGridDataMetricsOptions *>(poOptions);
3008 11 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
3009 : {
3010 0 : CPLError(CE_Failure, CPLE_AppDefined,
3011 : "Wrong value of nSizeOfStructure member");
3012 0 : return nullptr;
3013 : }
3014 11 : poOptionsNew = CPLMalloc(sizeof(GDALGridDataMetricsOptions));
3015 11 : memcpy(poOptionsNew, poOptions, sizeof(GDALGridDataMetricsOptions));
3016 :
3017 11 : if (poOptionsOld->nMinPointsPerQuadrant != 0 ||
3018 6 : poOptionsOld->nMaxPointsPerQuadrant != 0)
3019 : {
3020 5 : pfnGDALGridMethod = GDALGridDataMetricMaximumPerQuadrant;
3021 5 : bCreateQuadTree = true;
3022 : }
3023 : else
3024 : {
3025 6 : pfnGDALGridMethod = GDALGridDataMetricMaximum;
3026 6 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
3027 10 : poOptionsOld->dfAngle == 0.0 &&
3028 4 : (poOptionsOld->dfRadius1 > 0.0 ||
3029 1 : poOptionsOld->dfRadius2 > 0.0));
3030 : }
3031 :
3032 11 : break;
3033 : }
3034 8 : case GGA_MetricRange:
3035 : {
3036 8 : const auto poOptionsOld =
3037 : static_cast<const GDALGridDataMetricsOptions *>(poOptions);
3038 8 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
3039 : {
3040 0 : CPLError(CE_Failure, CPLE_AppDefined,
3041 : "Wrong value of nSizeOfStructure member");
3042 0 : return nullptr;
3043 : }
3044 8 : poOptionsNew = CPLMalloc(sizeof(GDALGridDataMetricsOptions));
3045 8 : memcpy(poOptionsNew, poOptions, sizeof(GDALGridDataMetricsOptions));
3046 :
3047 8 : if (poOptionsOld->nMinPointsPerQuadrant != 0 ||
3048 3 : poOptionsOld->nMaxPointsPerQuadrant != 0)
3049 : {
3050 5 : pfnGDALGridMethod = GDALGridDataMetricRangePerQuadrant;
3051 5 : bCreateQuadTree = true;
3052 : }
3053 : else
3054 : {
3055 3 : pfnGDALGridMethod = GDALGridDataMetricRange;
3056 3 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
3057 6 : poOptionsOld->dfAngle == 0.0 &&
3058 3 : (poOptionsOld->dfRadius1 > 0.0 ||
3059 1 : poOptionsOld->dfRadius2 > 0.0));
3060 : }
3061 :
3062 8 : break;
3063 : }
3064 10 : case GGA_MetricCount:
3065 : {
3066 10 : const auto poOptionsOld =
3067 : static_cast<const GDALGridDataMetricsOptions *>(poOptions);
3068 10 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
3069 : {
3070 0 : CPLError(CE_Failure, CPLE_AppDefined,
3071 : "Wrong value of nSizeOfStructure member");
3072 0 : return nullptr;
3073 : }
3074 10 : poOptionsNew = CPLMalloc(sizeof(GDALGridDataMetricsOptions));
3075 10 : memcpy(poOptionsNew, poOptions, sizeof(GDALGridDataMetricsOptions));
3076 :
3077 10 : if (poOptionsOld->nMinPointsPerQuadrant != 0 ||
3078 5 : poOptionsOld->nMaxPointsPerQuadrant != 0)
3079 : {
3080 5 : pfnGDALGridMethod = GDALGridDataMetricCountPerQuadrant;
3081 5 : bCreateQuadTree = true;
3082 : }
3083 : else
3084 : {
3085 5 : pfnGDALGridMethod = GDALGridDataMetricCount;
3086 5 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
3087 10 : poOptionsOld->dfAngle == 0.0 &&
3088 5 : (poOptionsOld->dfRadius1 > 0.0 ||
3089 0 : poOptionsOld->dfRadius2 > 0.0));
3090 : }
3091 :
3092 10 : break;
3093 : }
3094 8 : case GGA_MetricAverageDistance:
3095 : {
3096 8 : const auto poOptionsOld =
3097 : static_cast<const GDALGridDataMetricsOptions *>(poOptions);
3098 8 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
3099 : {
3100 0 : CPLError(CE_Failure, CPLE_AppDefined,
3101 : "Wrong value of nSizeOfStructure member");
3102 0 : return nullptr;
3103 : }
3104 8 : poOptionsNew = CPLMalloc(sizeof(GDALGridDataMetricsOptions));
3105 8 : memcpy(poOptionsNew, poOptions, sizeof(GDALGridDataMetricsOptions));
3106 :
3107 8 : if (poOptionsOld->nMinPointsPerQuadrant != 0 ||
3108 3 : poOptionsOld->nMaxPointsPerQuadrant != 0)
3109 : {
3110 5 : pfnGDALGridMethod =
3111 : GDALGridDataMetricAverageDistancePerQuadrant;
3112 5 : bCreateQuadTree = true;
3113 : }
3114 : else
3115 : {
3116 3 : pfnGDALGridMethod = GDALGridDataMetricAverageDistance;
3117 3 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
3118 6 : poOptionsOld->dfAngle == 0.0 &&
3119 3 : (poOptionsOld->dfRadius1 > 0.0 ||
3120 1 : poOptionsOld->dfRadius2 > 0.0));
3121 : }
3122 :
3123 8 : break;
3124 : }
3125 3 : case GGA_MetricAverageDistancePts:
3126 : {
3127 3 : const auto poOptionsOld =
3128 : static_cast<const GDALGridDataMetricsOptions *>(poOptions);
3129 3 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
3130 : {
3131 0 : CPLError(CE_Failure, CPLE_AppDefined,
3132 : "Wrong value of nSizeOfStructure member");
3133 0 : return nullptr;
3134 : }
3135 3 : poOptionsNew = CPLMalloc(sizeof(GDALGridDataMetricsOptions));
3136 3 : memcpy(poOptionsNew, poOptions, sizeof(GDALGridDataMetricsOptions));
3137 :
3138 3 : pfnGDALGridMethod = GDALGridDataMetricAverageDistancePts;
3139 6 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
3140 5 : poOptionsOld->dfAngle == 0.0 &&
3141 2 : (poOptionsOld->dfRadius1 > 0.0 ||
3142 0 : poOptionsOld->dfRadius2 > 0.0));
3143 :
3144 3 : break;
3145 : }
3146 2 : case GGA_Linear:
3147 : {
3148 2 : const auto poOptionsOld =
3149 : static_cast<const GDALGridLinearOptions *>(poOptions);
3150 2 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
3151 : {
3152 0 : CPLError(CE_Failure, CPLE_AppDefined,
3153 : "Wrong value of nSizeOfStructure member");
3154 0 : return nullptr;
3155 : }
3156 2 : poOptionsNew = CPLMalloc(sizeof(GDALGridLinearOptions));
3157 2 : memcpy(poOptionsNew, poOptions, sizeof(GDALGridLinearOptions));
3158 :
3159 2 : pfnGDALGridMethod = GDALGridLinear;
3160 2 : break;
3161 : }
3162 0 : default:
3163 0 : CPLError(CE_Failure, CPLE_IllegalArg,
3164 : "GDAL does not support gridding method %d", eAlgorithm);
3165 0 : return nullptr;
3166 : }
3167 :
3168 108 : if (pafXAligned == nullptr && !bCallerWillKeepPointArraysAlive)
3169 : {
3170 : double *padfXNew =
3171 0 : static_cast<double *>(VSI_MALLOC2_VERBOSE(nPoints, sizeof(double)));
3172 : double *padfYNew =
3173 0 : static_cast<double *>(VSI_MALLOC2_VERBOSE(nPoints, sizeof(double)));
3174 : double *padfZNew =
3175 0 : static_cast<double *>(VSI_MALLOC2_VERBOSE(nPoints, sizeof(double)));
3176 0 : if (padfXNew == nullptr || padfYNew == nullptr || padfZNew == nullptr)
3177 : {
3178 0 : VSIFree(padfXNew);
3179 0 : VSIFree(padfYNew);
3180 0 : VSIFree(padfZNew);
3181 0 : CPLFree(poOptionsNew);
3182 0 : return nullptr;
3183 : }
3184 0 : memcpy(padfXNew, padfX, nPoints * sizeof(double));
3185 0 : memcpy(padfYNew, padfY, nPoints * sizeof(double));
3186 0 : memcpy(padfZNew, padfZ, nPoints * sizeof(double));
3187 0 : padfX = padfXNew;
3188 0 : padfY = padfYNew;
3189 0 : padfZ = padfZNew;
3190 : }
3191 : GDALGridContext *psContext =
3192 108 : static_cast<GDALGridContext *>(CPLCalloc(1, sizeof(GDALGridContext)));
3193 108 : psContext->eAlgorithm = eAlgorithm;
3194 108 : psContext->poOptions = poOptionsNew;
3195 108 : psContext->pfnGDALGridMethod = pfnGDALGridMethod;
3196 108 : psContext->nPoints = nPoints;
3197 108 : psContext->pasGridPoints = nullptr;
3198 108 : psContext->sXYArrays.padfX = padfX;
3199 108 : psContext->sXYArrays.padfY = padfY;
3200 108 : psContext->sExtraParameters.hQuadTree = nullptr;
3201 108 : psContext->sExtraParameters.dfInitialSearchRadius = 0.0;
3202 108 : psContext->sExtraParameters.pafX = pafXAligned;
3203 108 : psContext->sExtraParameters.pafY = pafYAligned;
3204 108 : psContext->sExtraParameters.pafZ = pafZAligned;
3205 108 : psContext->sExtraParameters.psTriangulation = nullptr;
3206 108 : psContext->sExtraParameters.nInitialFacetIdx = 0;
3207 108 : psContext->padfX = pafXAligned ? nullptr : const_cast<double *>(padfX);
3208 108 : psContext->padfY = pafXAligned ? nullptr : const_cast<double *>(padfY);
3209 108 : psContext->padfZ = pafXAligned ? nullptr : const_cast<double *>(padfZ);
3210 108 : psContext->bFreePadfXYZArrays =
3211 108 : pafXAligned ? false : !bCallerWillKeepPointArraysAlive;
3212 :
3213 : /* -------------------------------------------------------------------- */
3214 : /* Create quadtree if requested and possible. */
3215 : /* -------------------------------------------------------------------- */
3216 108 : if (bCreateQuadTree)
3217 : {
3218 72 : GDALGridContextCreateQuadTree(psContext);
3219 72 : if (psContext->sExtraParameters.hQuadTree == nullptr &&
3220 0 : (eAlgorithm == GGA_InverseDistanceToAPowerNearestNeighbor ||
3221 : pfnGDALGridMethod == GDALGridMovingAveragePerQuadrant))
3222 : {
3223 : // shouldn't happen unless memory allocation failure occurs
3224 0 : GDALGridContextFree(psContext);
3225 0 : return nullptr;
3226 : }
3227 : }
3228 :
3229 : /* -------------------------------------------------------------------- */
3230 : /* Pre-compute extra parameters in GDALGridExtraParameters */
3231 : /* -------------------------------------------------------------------- */
3232 108 : if (eAlgorithm == GGA_InverseDistanceToAPowerNearestNeighbor)
3233 : {
3234 11 : const double dfPower =
3235 : static_cast<
3236 : const GDALGridInverseDistanceToAPowerNearestNeighborOptions *>(
3237 : poOptions)
3238 : ->dfPower;
3239 11 : psContext->sExtraParameters.dfPowerDiv2PreComp = dfPower / 2;
3240 :
3241 11 : const double dfRadius =
3242 : static_cast<
3243 : const GDALGridInverseDistanceToAPowerNearestNeighborOptions *>(
3244 : poOptions)
3245 : ->dfRadius;
3246 11 : psContext->sExtraParameters.dfRadiusPower2PreComp = pow(dfRadius, 2);
3247 : }
3248 :
3249 108 : if (eAlgorithm == GGA_Linear)
3250 : {
3251 2 : psContext->sExtraParameters.psTriangulation =
3252 2 : GDALTriangulationCreateDelaunay(nPoints, padfX, padfY);
3253 2 : if (psContext->sExtraParameters.psTriangulation == nullptr)
3254 : {
3255 0 : GDALGridContextFree(psContext);
3256 0 : return nullptr;
3257 : }
3258 2 : GDALTriangulationComputeBarycentricCoefficients(
3259 : psContext->sExtraParameters.psTriangulation, padfX, padfY);
3260 : }
3261 :
3262 : /* -------------------------------------------------------------------- */
3263 : /* Start thread pool. */
3264 : /* -------------------------------------------------------------------- */
3265 108 : const char *pszThreads = CPLGetConfigOption("GDAL_NUM_THREADS", "ALL_CPUS");
3266 108 : int nThreads = 0;
3267 108 : if (EQUAL(pszThreads, "ALL_CPUS"))
3268 106 : nThreads = CPLGetNumCPUs();
3269 : else
3270 2 : nThreads = atoi(pszThreads);
3271 108 : if (nThreads > 128)
3272 0 : nThreads = 128;
3273 108 : if (nThreads > 1)
3274 : {
3275 106 : psContext->poWorkerThreadPool = new CPLWorkerThreadPool();
3276 106 : if (!psContext->poWorkerThreadPool->Setup(nThreads, nullptr, nullptr))
3277 : {
3278 0 : delete psContext->poWorkerThreadPool;
3279 0 : psContext->poWorkerThreadPool = nullptr;
3280 : }
3281 : else
3282 : {
3283 106 : CPLDebug("GDAL_GRID", "Using %d threads", nThreads);
3284 : }
3285 : }
3286 : else
3287 2 : psContext->poWorkerThreadPool = nullptr;
3288 :
3289 108 : return psContext;
3290 : }
3291 :
3292 : /************************************************************************/
3293 : /* GDALGridContextCreateQuadTree() */
3294 : /************************************************************************/
3295 :
3296 74 : void GDALGridContextCreateQuadTree(GDALGridContext *psContext)
3297 : {
3298 74 : const GUInt32 nPoints = psContext->nPoints;
3299 74 : psContext->pasGridPoints = static_cast<GDALGridPoint *>(
3300 74 : VSI_MALLOC2_VERBOSE(nPoints, sizeof(GDALGridPoint)));
3301 74 : if (psContext->pasGridPoints != nullptr)
3302 : {
3303 74 : const double *const padfX = psContext->padfX;
3304 74 : const double *const padfY = psContext->padfY;
3305 :
3306 : // Determine point extents.
3307 : CPLRectObj sRect;
3308 74 : sRect.minx = padfX[0];
3309 74 : sRect.miny = padfY[0];
3310 74 : sRect.maxx = padfX[0];
3311 74 : sRect.maxy = padfY[0];
3312 28031 : for (GUInt32 i = 1; i < nPoints; i++)
3313 : {
3314 27957 : if (padfX[i] < sRect.minx)
3315 40 : sRect.minx = padfX[i];
3316 27957 : if (padfY[i] < sRect.miny)
3317 788 : sRect.miny = padfY[i];
3318 27957 : if (padfX[i] > sRect.maxx)
3319 774 : sRect.maxx = padfX[i];
3320 27957 : if (padfY[i] > sRect.maxy)
3321 10 : sRect.maxy = padfY[i];
3322 : }
3323 :
3324 : // Initial value for search radius is the typical dimension of a
3325 : // "pixel" of the point array (assuming rather uniform distribution).
3326 74 : psContext->sExtraParameters.dfInitialSearchRadius = sqrt(
3327 74 : (sRect.maxx - sRect.minx) * (sRect.maxy - sRect.miny) / nPoints);
3328 :
3329 74 : psContext->sExtraParameters.hQuadTree =
3330 74 : CPLQuadTreeCreate(&sRect, GDALGridGetPointBounds);
3331 :
3332 28105 : for (GUInt32 i = 0; i < nPoints; i++)
3333 : {
3334 28031 : psContext->pasGridPoints[i].psXYArrays = &(psContext->sXYArrays);
3335 28031 : psContext->pasGridPoints[i].i = i;
3336 28031 : CPLQuadTreeInsert(psContext->sExtraParameters.hQuadTree,
3337 28031 : psContext->pasGridPoints + i);
3338 : }
3339 : }
3340 74 : }
3341 :
3342 : /************************************************************************/
3343 : /* GDALGridContextFree() */
3344 : /************************************************************************/
3345 :
3346 : /**
3347 : * Free a context used created by GDALGridContextCreate()
3348 : *
3349 : * @param psContext the context.
3350 : *
3351 : * @since GDAL 2.1
3352 : */
3353 108 : void GDALGridContextFree(GDALGridContext *psContext)
3354 : {
3355 108 : if (psContext)
3356 : {
3357 108 : CPLFree(psContext->poOptions);
3358 108 : CPLFree(psContext->pasGridPoints);
3359 108 : if (psContext->sExtraParameters.hQuadTree != nullptr)
3360 74 : CPLQuadTreeDestroy(psContext->sExtraParameters.hQuadTree);
3361 108 : if (psContext->bFreePadfXYZArrays)
3362 : {
3363 0 : CPLFree(psContext->padfX);
3364 0 : CPLFree(psContext->padfY);
3365 0 : CPLFree(psContext->padfZ);
3366 : }
3367 108 : VSIFreeAligned(psContext->sExtraParameters.pafX);
3368 108 : VSIFreeAligned(psContext->sExtraParameters.pafY);
3369 108 : VSIFreeAligned(psContext->sExtraParameters.pafZ);
3370 108 : if (psContext->sExtraParameters.psTriangulation)
3371 2 : GDALTriangulationFree(psContext->sExtraParameters.psTriangulation);
3372 108 : delete psContext->poWorkerThreadPool;
3373 108 : CPLFree(psContext);
3374 : }
3375 108 : }
3376 :
3377 : /************************************************************************/
3378 : /* GDALGridContextProcess() */
3379 : /************************************************************************/
3380 :
3381 : /**
3382 : * Do the gridding of a window of a raster.
3383 : *
3384 : * This function takes the gridding context as input to preprare computation
3385 : * of regular grid (or call it a raster) from these scattered data.
3386 : * You should supply the extent of the output grid and allocate array
3387 : * sufficient to hold such a grid.
3388 : *
3389 : * @param psContext Gridding context.
3390 : * @param dfXMin Lowest X border of output grid.
3391 : * @param dfXMax Highest X border of output grid.
3392 : * @param dfYMin Lowest Y border of output grid.
3393 : * @param dfYMax Highest Y border of output grid.
3394 : * @param nXSize Number of columns in output grid.
3395 : * @param nYSize Number of rows in output grid.
3396 : * @param eType Data type of output array.
3397 : * @param pData Pointer to array where the computed grid will be stored.
3398 : * @param pfnProgress a GDALProgressFunc() compatible callback function for
3399 : * reporting progress or NULL.
3400 : * @param pProgressArg argument to be passed to pfnProgress. May be NULL.
3401 : *
3402 : * @return CE_None on success or CE_Failure if something goes wrong.
3403 : *
3404 : * @since GDAL 2.1
3405 : */
3406 :
3407 108 : CPLErr GDALGridContextProcess(GDALGridContext *psContext, double dfXMin,
3408 : double dfXMax, double dfYMin, double dfYMax,
3409 : GUInt32 nXSize, GUInt32 nYSize,
3410 : GDALDataType eType, void *pData,
3411 : GDALProgressFunc pfnProgress, void *pProgressArg)
3412 : {
3413 108 : CPLAssert(psContext);
3414 108 : CPLAssert(pData);
3415 :
3416 108 : if (nXSize == 0 || nYSize == 0)
3417 : {
3418 0 : CPLError(CE_Failure, CPLE_IllegalArg,
3419 : "Output raster dimensions should have non-zero size.");
3420 0 : return CE_Failure;
3421 : }
3422 :
3423 108 : const double dfDeltaX = (dfXMax - dfXMin) / nXSize;
3424 108 : const double dfDeltaY = (dfYMax - dfYMin) / nYSize;
3425 :
3426 : // For linear, check if we will need to fallback to nearest neighbour
3427 : // by sampling along the edges. If all points on edges are within
3428 : // triangles, then interior points will also be.
3429 108 : if (psContext->eAlgorithm == GGA_Linear &&
3430 2 : psContext->sExtraParameters.hQuadTree == nullptr)
3431 : {
3432 2 : bool bNeedNearest = false;
3433 2 : int nStartLeft = 0;
3434 2 : int nStartRight = 0;
3435 2 : const double dfXPointMin = dfXMin + (0 + 0.5) * dfDeltaX;
3436 2 : const double dfXPointMax = dfXMin + (nXSize - 1 + 0.5) * dfDeltaX;
3437 4 : for (GUInt32 nYPoint = 0; !bNeedNearest && nYPoint < nYSize; nYPoint++)
3438 : {
3439 2 : const double dfYPoint = dfYMin + (nYPoint + 0.5) * dfDeltaY;
3440 :
3441 2 : if (!GDALTriangulationFindFacetDirected(
3442 2 : psContext->sExtraParameters.psTriangulation, nStartLeft,
3443 : dfXPointMin, dfYPoint, &nStartLeft))
3444 : {
3445 2 : bNeedNearest = true;
3446 : }
3447 2 : if (!GDALTriangulationFindFacetDirected(
3448 2 : psContext->sExtraParameters.psTriangulation, nStartRight,
3449 : dfXPointMax, dfYPoint, &nStartRight))
3450 : {
3451 2 : bNeedNearest = true;
3452 : }
3453 : }
3454 2 : int nStartTop = 0;
3455 2 : int nStartBottom = 0;
3456 2 : const double dfYPointMin = dfYMin + (0 + 0.5) * dfDeltaY;
3457 2 : const double dfYPointMax = dfYMin + (nYSize - 1 + 0.5) * dfDeltaY;
3458 2 : for (GUInt32 nXPoint = 1; !bNeedNearest && nXPoint + 1 < nXSize;
3459 : nXPoint++)
3460 : {
3461 0 : const double dfXPoint = dfXMin + (nXPoint + 0.5) * dfDeltaX;
3462 :
3463 0 : if (!GDALTriangulationFindFacetDirected(
3464 0 : psContext->sExtraParameters.psTriangulation, nStartTop,
3465 : dfXPoint, dfYPointMin, &nStartTop))
3466 : {
3467 0 : bNeedNearest = true;
3468 : }
3469 0 : if (!GDALTriangulationFindFacetDirected(
3470 0 : psContext->sExtraParameters.psTriangulation, nStartBottom,
3471 : dfXPoint, dfYPointMax, &nStartBottom))
3472 : {
3473 0 : bNeedNearest = true;
3474 : }
3475 : }
3476 2 : if (bNeedNearest)
3477 : {
3478 2 : CPLDebug("GDAL_GRID", "Will need nearest neighbour");
3479 2 : GDALGridContextCreateQuadTree(psContext);
3480 : }
3481 : }
3482 :
3483 108 : int nCounter = 0;
3484 108 : volatile int bStop = FALSE;
3485 : GDALGridJob sJob;
3486 108 : sJob.nYStart = 0;
3487 108 : sJob.pabyData = static_cast<GByte *>(pData);
3488 108 : sJob.nYStep = 1;
3489 108 : sJob.nXSize = nXSize;
3490 108 : sJob.nYSize = nYSize;
3491 108 : sJob.dfXMin = dfXMin;
3492 108 : sJob.dfYMin = dfYMin;
3493 108 : sJob.dfDeltaX = dfDeltaX;
3494 108 : sJob.dfDeltaY = dfDeltaY;
3495 108 : sJob.nPoints = psContext->nPoints;
3496 108 : sJob.padfX = psContext->padfX;
3497 108 : sJob.padfY = psContext->padfY;
3498 108 : sJob.padfZ = psContext->padfZ;
3499 108 : sJob.poOptions = psContext->poOptions;
3500 108 : sJob.pfnGDALGridMethod = psContext->pfnGDALGridMethod;
3501 108 : sJob.psExtraParameters = &psContext->sExtraParameters;
3502 108 : sJob.pfnProgress = nullptr;
3503 108 : sJob.eType = eType;
3504 108 : sJob.pfnRealProgress = pfnProgress;
3505 108 : sJob.pRealProgressArg = pProgressArg;
3506 108 : sJob.pnCounter = &nCounter;
3507 108 : sJob.pbStop = &bStop;
3508 108 : sJob.hCond = nullptr;
3509 108 : sJob.hCondMutex = nullptr;
3510 :
3511 108 : if (psContext->poWorkerThreadPool == nullptr)
3512 : {
3513 2 : if (sJob.pfnRealProgress != nullptr &&
3514 2 : sJob.pfnRealProgress != GDALDummyProgress)
3515 : {
3516 2 : sJob.pfnProgress = GDALGridProgressMonoThread;
3517 : }
3518 :
3519 2 : GDALGridJobProcess(&sJob);
3520 : }
3521 : else
3522 : {
3523 106 : int nThreads = psContext->poWorkerThreadPool->GetThreadCount();
3524 : GDALGridJob *pasJobs = static_cast<GDALGridJob *>(
3525 106 : CPLMalloc(sizeof(GDALGridJob) * nThreads));
3526 :
3527 106 : sJob.nYStep = nThreads;
3528 106 : sJob.hCondMutex = CPLCreateMutex(); /* and implicitly take the mutex */
3529 106 : sJob.hCond = CPLCreateCond();
3530 106 : sJob.pfnProgress = GDALGridProgressMultiThread;
3531 :
3532 : /* --------------------------------------------------------------------
3533 : */
3534 : /* Start threads. */
3535 : /* --------------------------------------------------------------------
3536 : */
3537 530 : for (int i = 0; i < nThreads && !bStop; i++)
3538 : {
3539 424 : memcpy(&pasJobs[i], &sJob, sizeof(GDALGridJob));
3540 424 : pasJobs[i].nYStart = i;
3541 424 : psContext->poWorkerThreadPool->SubmitJob(GDALGridJobProcess,
3542 424 : &pasJobs[i]);
3543 : }
3544 :
3545 : /* --------------------------------------------------------------------
3546 : */
3547 : /* Report progress. */
3548 : /* --------------------------------------------------------------------
3549 : */
3550 527 : while (*(sJob.pnCounter) < static_cast<int>(nYSize) && !bStop)
3551 : {
3552 421 : CPLCondWait(sJob.hCond, sJob.hCondMutex);
3553 :
3554 421 : int nLocalCounter = *(sJob.pnCounter);
3555 421 : CPLReleaseMutex(sJob.hCondMutex);
3556 :
3557 835 : if (pfnProgress != nullptr &&
3558 414 : !pfnProgress(nLocalCounter / static_cast<double>(nYSize), "",
3559 : pProgressArg))
3560 : {
3561 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
3562 0 : bStop = TRUE;
3563 : }
3564 :
3565 421 : CPLAcquireMutex(sJob.hCondMutex, 1.0);
3566 : }
3567 :
3568 : // Release mutex before joining threads, otherwise they will dead-lock
3569 : // forever in GDALGridProgressMultiThread().
3570 106 : CPLReleaseMutex(sJob.hCondMutex);
3571 :
3572 : /* --------------------------------------------------------------------
3573 : */
3574 : /* Wait for all threads to complete and finish. */
3575 : /* --------------------------------------------------------------------
3576 : */
3577 106 : psContext->poWorkerThreadPool->WaitCompletion();
3578 :
3579 106 : CPLFree(pasJobs);
3580 106 : CPLDestroyCond(sJob.hCond);
3581 106 : CPLDestroyMutex(sJob.hCondMutex);
3582 : }
3583 :
3584 108 : return bStop ? CE_Failure : CE_None;
3585 : }
3586 :
3587 : /************************************************************************/
3588 : /* GDALGridCreate() */
3589 : /************************************************************************/
3590 :
3591 : /**
3592 : * Create regular grid from the scattered data.
3593 : *
3594 : * This function takes the arrays of X and Y coordinates and corresponding Z
3595 : * values as input and computes regular grid (or call it a raster) from these
3596 : * scattered data. You should supply geometry and extent of the output grid
3597 : * and allocate array sufficient to hold such a grid.
3598 : *
3599 : * Starting with GDAL 1.10, it is possible to set the GDAL_NUM_THREADS
3600 : * configuration option to parallelize the processing. The value to set is
3601 : * the number of worker threads, or ALL_CPUS to use all the cores/CPUs of the
3602 : * computer (default value).
3603 : *
3604 : * Starting with GDAL 1.10, on Intel/AMD i386/x86_64 architectures, some
3605 : * gridding methods will be optimized with SSE instructions (provided GDAL
3606 : * has been compiled with such support, and it is available at runtime).
3607 : * Currently, only 'invdist' algorithm with default parameters has an optimized
3608 : * implementation.
3609 : * This can provide substantial speed-up, but sometimes at the expense of
3610 : * reduced floating point precision. This can be disabled by setting the
3611 : * GDAL_USE_SSE configuration option to NO.
3612 : * Starting with GDAL 1.11, a further optimized version can use the AVX
3613 : * instruction set. This can be disabled by setting the GDAL_USE_AVX
3614 : * configuration option to NO.
3615 : *
3616 : * Note: it will be more efficient to use GDALGridContextCreate(),
3617 : * GDALGridContextProcess() and GDALGridContextFree() when doing repeated
3618 : * gridding operations with the same algorithm, parameters and points, and
3619 : * moving the window in the output grid.
3620 : *
3621 : * @param eAlgorithm Gridding method.
3622 : * @param poOptions Options to control chosen gridding method.
3623 : * @param nPoints Number of elements in input arrays.
3624 : * @param padfX Input array of X coordinates.
3625 : * @param padfY Input array of Y coordinates.
3626 : * @param padfZ Input array of Z values.
3627 : * @param dfXMin Lowest X border of output grid.
3628 : * @param dfXMax Highest X border of output grid.
3629 : * @param dfYMin Lowest Y border of output grid.
3630 : * @param dfYMax Highest Y border of output grid.
3631 : * @param nXSize Number of columns in output grid.
3632 : * @param nYSize Number of rows in output grid.
3633 : * @param eType Data type of output array.
3634 : * @param pData Pointer to array where the computed grid will be stored.
3635 : * @param pfnProgress a GDALProgressFunc() compatible callback function for
3636 : * reporting progress or NULL.
3637 : * @param pProgressArg argument to be passed to pfnProgress. May be NULL.
3638 : *
3639 : * @return CE_None on success or CE_Failure if something goes wrong.
3640 : */
3641 :
3642 5 : CPLErr GDALGridCreate(GDALGridAlgorithm eAlgorithm, const void *poOptions,
3643 : GUInt32 nPoints, const double *padfX, const double *padfY,
3644 : const double *padfZ, double dfXMin, double dfXMax,
3645 : double dfYMin, double dfYMax, GUInt32 nXSize,
3646 : GUInt32 nYSize, GDALDataType eType, void *pData,
3647 : GDALProgressFunc pfnProgress, void *pProgressArg)
3648 : {
3649 5 : GDALGridContext *psContext = GDALGridContextCreate(
3650 : eAlgorithm, poOptions, nPoints, padfX, padfY, padfZ, TRUE);
3651 5 : CPLErr eErr = CE_Failure;
3652 5 : if (psContext)
3653 : {
3654 5 : eErr = GDALGridContextProcess(psContext, dfXMin, dfXMax, dfYMin, dfYMax,
3655 : nXSize, nYSize, eType, pData, pfnProgress,
3656 : pProgressArg);
3657 : }
3658 :
3659 5 : GDALGridContextFree(psContext);
3660 5 : return eErr;
3661 : }
3662 :
3663 : /************************************************************************/
3664 : /* GDALGridParseAlgorithmAndOptions() */
3665 : /************************************************************************/
3666 :
3667 : /** Translates mnemonic gridding algorithm names into GDALGridAlgorithm code,
3668 : * parse control parameters and assign defaults.
3669 : */
3670 229 : CPLErr GDALGridParseAlgorithmAndOptions(const char *pszAlgorithm,
3671 : GDALGridAlgorithm *peAlgorithm,
3672 : void **ppOptions)
3673 : {
3674 229 : CPLAssert(pszAlgorithm);
3675 229 : CPLAssert(peAlgorithm);
3676 229 : CPLAssert(ppOptions);
3677 :
3678 229 : *ppOptions = nullptr;
3679 :
3680 229 : char **papszParams = CSLTokenizeString2(pszAlgorithm, ":", FALSE);
3681 :
3682 229 : if (CSLCount(papszParams) < 1)
3683 : {
3684 0 : CSLDestroy(papszParams);
3685 0 : return CE_Failure;
3686 : }
3687 :
3688 229 : if (EQUAL(papszParams[0], szAlgNameInvDist))
3689 : {
3690 276 : if (CSLFetchNameValue(papszParams, "min_points_per_quadrant") ||
3691 137 : CSLFetchNameValue(papszParams, "max_points_per_quadrant"))
3692 : {
3693 : // Remap invdist to invdistnn if per quadrant is required
3694 2 : if (CSLFetchNameValue(papszParams, "radius") == nullptr)
3695 : {
3696 : const double dfRadius1 =
3697 0 : CPLAtofM(CSLFetchNameValueDef(papszParams, "radius1", "1"));
3698 : const double dfRadius2 =
3699 0 : CPLAtofM(CSLFetchNameValueDef(papszParams, "radius2", "1"));
3700 0 : if (dfRadius1 != dfRadius2)
3701 : {
3702 0 : CPLError(CE_Failure, CPLE_NotSupported,
3703 : "radius1 != radius2 not supported when "
3704 : "min_points_per_quadrant and/or "
3705 : "max_points_per_quadrant is specified");
3706 0 : CSLDestroy(papszParams);
3707 0 : return CE_Failure;
3708 : }
3709 : }
3710 :
3711 2 : if (CPLAtofM(CSLFetchNameValueDef(papszParams, "angle", "0")) != 0)
3712 : {
3713 0 : CPLError(CE_Failure, CPLE_NotSupported,
3714 : "angle != 0 not supported when "
3715 : "min_points_per_quadrant and/or "
3716 : "max_points_per_quadrant is specified");
3717 0 : CSLDestroy(papszParams);
3718 0 : return CE_Failure;
3719 : }
3720 :
3721 2 : char **papszNewParams = CSLAddString(nullptr, "invdistnn");
3722 2 : if (CSLFetchNameValue(papszParams, "radius") == nullptr)
3723 : {
3724 0 : papszNewParams = CSLSetNameValue(
3725 : papszNewParams, "radius",
3726 : CSLFetchNameValueDef(papszParams, "radius1", "1"));
3727 : }
3728 12 : for (const char *pszItem : {"radius", "power", "smoothing",
3729 14 : "max_points", "min_points", "nodata"})
3730 : {
3731 12 : const char *pszValue = CSLFetchNameValue(papszParams, pszItem);
3732 12 : if (pszValue)
3733 : papszNewParams =
3734 6 : CSLSetNameValue(papszNewParams, pszItem, pszValue);
3735 : }
3736 2 : CSLDestroy(papszParams);
3737 2 : papszParams = papszNewParams;
3738 :
3739 2 : *peAlgorithm = GGA_InverseDistanceToAPowerNearestNeighbor;
3740 : }
3741 : else
3742 : {
3743 137 : *peAlgorithm = GGA_InverseDistanceToAPower;
3744 : }
3745 : }
3746 90 : else if (EQUAL(papszParams[0], szAlgNameInvDistNearestNeighbor))
3747 : {
3748 9 : *peAlgorithm = GGA_InverseDistanceToAPowerNearestNeighbor;
3749 : }
3750 81 : else if (EQUAL(papszParams[0], szAlgNameAverage))
3751 : {
3752 13 : *peAlgorithm = GGA_MovingAverage;
3753 : }
3754 68 : else if (EQUAL(papszParams[0], szAlgNameNearest))
3755 : {
3756 14 : *peAlgorithm = GGA_NearestNeighbor;
3757 : }
3758 54 : else if (EQUAL(papszParams[0], szAlgNameMinimum))
3759 : {
3760 11 : *peAlgorithm = GGA_MetricMinimum;
3761 : }
3762 43 : else if (EQUAL(papszParams[0], szAlgNameMaximum))
3763 : {
3764 11 : *peAlgorithm = GGA_MetricMaximum;
3765 : }
3766 32 : else if (EQUAL(papszParams[0], szAlgNameRange))
3767 : {
3768 8 : *peAlgorithm = GGA_MetricRange;
3769 : }
3770 24 : else if (EQUAL(papszParams[0], szAlgNameCount))
3771 : {
3772 10 : *peAlgorithm = GGA_MetricCount;
3773 : }
3774 14 : else if (EQUAL(papszParams[0], szAlgNameAverageDistance))
3775 : {
3776 8 : *peAlgorithm = GGA_MetricAverageDistance;
3777 : }
3778 6 : else if (EQUAL(papszParams[0], szAlgNameAverageDistancePts))
3779 : {
3780 3 : *peAlgorithm = GGA_MetricAverageDistancePts;
3781 : }
3782 3 : else if (EQUAL(papszParams[0], szAlgNameLinear))
3783 : {
3784 2 : *peAlgorithm = GGA_Linear;
3785 : }
3786 : else
3787 : {
3788 1 : CPLError(CE_Failure, CPLE_IllegalArg,
3789 : "Unsupported gridding method \"%s\"", papszParams[0]);
3790 1 : CSLDestroy(papszParams);
3791 1 : return CE_Failure;
3792 : }
3793 :
3794 : /* -------------------------------------------------------------------- */
3795 : /* Parse algorithm parameters and assign defaults. */
3796 : /* -------------------------------------------------------------------- */
3797 228 : const char *const *papszKnownOptions = nullptr;
3798 :
3799 228 : switch (*peAlgorithm)
3800 : {
3801 137 : case GGA_InverseDistanceToAPower:
3802 : default:
3803 : {
3804 137 : *ppOptions =
3805 137 : CPLMalloc(sizeof(GDALGridInverseDistanceToAPowerOptions));
3806 :
3807 137 : GDALGridInverseDistanceToAPowerOptions *const poPowerOpts =
3808 : static_cast<GDALGridInverseDistanceToAPowerOptions *>(
3809 : *ppOptions);
3810 :
3811 137 : poPowerOpts->nSizeOfStructure = sizeof(*poPowerOpts);
3812 :
3813 137 : const char *pszValue = CSLFetchNameValue(papszParams, "power");
3814 137 : poPowerOpts->dfPower = pszValue ? CPLAtofM(pszValue) : 2.0;
3815 :
3816 137 : pszValue = CSLFetchNameValue(papszParams, "smoothing");
3817 137 : poPowerOpts->dfSmoothing = pszValue ? CPLAtofM(pszValue) : 0.0;
3818 :
3819 137 : pszValue = CSLFetchNameValue(papszParams, "radius");
3820 137 : if (pszValue)
3821 : {
3822 0 : poPowerOpts->dfRadius1 = CPLAtofM(pszValue);
3823 0 : poPowerOpts->dfRadius2 = poPowerOpts->dfRadius1;
3824 : }
3825 : else
3826 : {
3827 137 : pszValue = CSLFetchNameValue(papszParams, "radius1");
3828 137 : poPowerOpts->dfRadius1 = pszValue ? CPLAtofM(pszValue) : 0.0;
3829 :
3830 137 : pszValue = CSLFetchNameValue(papszParams, "radius2");
3831 137 : poPowerOpts->dfRadius2 = pszValue ? CPLAtofM(pszValue) : 0.0;
3832 : }
3833 :
3834 137 : pszValue = CSLFetchNameValue(papszParams, "angle");
3835 137 : poPowerOpts->dfAngle = pszValue ? CPLAtofM(pszValue) : 0.0;
3836 :
3837 137 : pszValue = CSLFetchNameValue(papszParams, "max_points");
3838 137 : poPowerOpts->nMaxPoints =
3839 137 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3840 :
3841 137 : pszValue = CSLFetchNameValue(papszParams, "min_points");
3842 137 : poPowerOpts->nMinPoints =
3843 137 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3844 :
3845 137 : pszValue = CSLFetchNameValue(papszParams, "nodata");
3846 137 : poPowerOpts->dfNoDataValue = pszValue ? CPLAtofM(pszValue) : 0.0;
3847 :
3848 : static const char *const apszKnownOptions[] = {
3849 : "power", "smoothing", "radius1", "radius2", "angle",
3850 : "max_points", "min_points", "nodata", nullptr};
3851 137 : papszKnownOptions = apszKnownOptions;
3852 :
3853 137 : break;
3854 : }
3855 11 : case GGA_InverseDistanceToAPowerNearestNeighbor:
3856 : {
3857 11 : *ppOptions = CPLMalloc(
3858 : sizeof(GDALGridInverseDistanceToAPowerNearestNeighborOptions));
3859 :
3860 : GDALGridInverseDistanceToAPowerNearestNeighborOptions
3861 11 : *const poPowerOpts = static_cast<
3862 : GDALGridInverseDistanceToAPowerNearestNeighborOptions *>(
3863 : *ppOptions);
3864 :
3865 11 : poPowerOpts->nSizeOfStructure = sizeof(*poPowerOpts);
3866 :
3867 11 : const char *pszValue = CSLFetchNameValue(papszParams, "power");
3868 11 : poPowerOpts->dfPower = pszValue ? CPLAtofM(pszValue) : 2.0;
3869 :
3870 11 : pszValue = CSLFetchNameValue(papszParams, "smoothing");
3871 11 : poPowerOpts->dfSmoothing = pszValue ? CPLAtofM(pszValue) : 0.0;
3872 :
3873 11 : pszValue = CSLFetchNameValue(papszParams, "radius");
3874 11 : poPowerOpts->dfRadius = pszValue ? CPLAtofM(pszValue) : 1.0;
3875 11 : if (!(poPowerOpts->dfRadius > 0))
3876 : {
3877 0 : CPLError(CE_Failure, CPLE_IllegalArg,
3878 : "Radius value should be strictly positive");
3879 0 : CSLDestroy(papszParams);
3880 0 : return CE_Failure;
3881 : }
3882 :
3883 11 : pszValue = CSLFetchNameValue(papszParams, "max_points");
3884 11 : poPowerOpts->nMaxPoints =
3885 11 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 12);
3886 :
3887 11 : pszValue = CSLFetchNameValue(papszParams, "min_points");
3888 11 : poPowerOpts->nMinPoints =
3889 11 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3890 :
3891 11 : pszValue = CSLFetchNameValue(papszParams, "nodata");
3892 11 : poPowerOpts->dfNoDataValue = pszValue ? CPLAtofM(pszValue) : 0.0;
3893 :
3894 : pszValue =
3895 11 : CSLFetchNameValue(papszParams, "min_points_per_quadrant");
3896 11 : poPowerOpts->nMinPointsPerQuadrant =
3897 11 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3898 :
3899 : pszValue =
3900 11 : CSLFetchNameValue(papszParams, "max_points_per_quadrant");
3901 11 : poPowerOpts->nMaxPointsPerQuadrant =
3902 11 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3903 :
3904 : static const char *const apszKnownOptions[] = {
3905 : "power",
3906 : "smoothing",
3907 : "radius",
3908 : "max_points",
3909 : "min_points",
3910 : "nodata",
3911 : "min_points_per_quadrant",
3912 : "max_points_per_quadrant",
3913 : nullptr};
3914 11 : papszKnownOptions = apszKnownOptions;
3915 :
3916 11 : break;
3917 : }
3918 13 : case GGA_MovingAverage:
3919 : {
3920 13 : *ppOptions = CPLMalloc(sizeof(GDALGridMovingAverageOptions));
3921 :
3922 13 : GDALGridMovingAverageOptions *const poAverageOpts =
3923 : static_cast<GDALGridMovingAverageOptions *>(*ppOptions);
3924 :
3925 13 : poAverageOpts->nSizeOfStructure = sizeof(*poAverageOpts);
3926 :
3927 13 : const char *pszValue = CSLFetchNameValue(papszParams, "radius");
3928 13 : if (pszValue)
3929 : {
3930 6 : poAverageOpts->dfRadius1 = CPLAtofM(pszValue);
3931 6 : poAverageOpts->dfRadius2 = poAverageOpts->dfRadius1;
3932 : }
3933 : else
3934 : {
3935 7 : pszValue = CSLFetchNameValue(papszParams, "radius1");
3936 7 : poAverageOpts->dfRadius1 = pszValue ? CPLAtofM(pszValue) : 0.0;
3937 :
3938 7 : pszValue = CSLFetchNameValue(papszParams, "radius2");
3939 7 : poAverageOpts->dfRadius2 = pszValue ? CPLAtofM(pszValue) : 0.0;
3940 : }
3941 :
3942 13 : pszValue = CSLFetchNameValue(papszParams, "angle");
3943 13 : poAverageOpts->dfAngle = pszValue ? CPLAtofM(pszValue) : 0.0;
3944 :
3945 13 : pszValue = CSLFetchNameValue(papszParams, "min_points");
3946 13 : poAverageOpts->nMinPoints =
3947 13 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3948 :
3949 13 : pszValue = CSLFetchNameValue(papszParams, "max_points");
3950 13 : poAverageOpts->nMaxPoints =
3951 13 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3952 :
3953 13 : pszValue = CSLFetchNameValue(papszParams, "nodata");
3954 13 : poAverageOpts->dfNoDataValue = pszValue ? CPLAtofM(pszValue) : 0.0;
3955 :
3956 : pszValue =
3957 13 : CSLFetchNameValue(papszParams, "min_points_per_quadrant");
3958 13 : poAverageOpts->nMinPointsPerQuadrant =
3959 13 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3960 :
3961 : pszValue =
3962 13 : CSLFetchNameValue(papszParams, "max_points_per_quadrant");
3963 13 : poAverageOpts->nMaxPointsPerQuadrant =
3964 13 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3965 :
3966 13 : if (poAverageOpts->nMinPointsPerQuadrant != 0 ||
3967 7 : poAverageOpts->nMaxPointsPerQuadrant != 0)
3968 : {
3969 6 : if (!(poAverageOpts->dfRadius1 > 0) ||
3970 6 : !(poAverageOpts->dfRadius2 > 0))
3971 : {
3972 0 : CPLError(CE_Failure, CPLE_IllegalArg,
3973 : "Radius value should be strictly positive when "
3974 : "per quadrant parameters are specified");
3975 0 : CSLDestroy(papszParams);
3976 0 : return CE_Failure;
3977 : }
3978 6 : if (poAverageOpts->dfAngle != 0)
3979 : {
3980 0 : CPLError(CE_Failure, CPLE_NotSupported,
3981 : "angle != 0 not supported when "
3982 : "per quadrant parameters are specified");
3983 0 : CSLDestroy(papszParams);
3984 0 : return CE_Failure;
3985 : }
3986 : }
3987 7 : else if (poAverageOpts->nMaxPoints > 0)
3988 : {
3989 0 : CPLError(CE_Warning, CPLE_AppDefined,
3990 : "max_points is ignored unless one of "
3991 : "min_points_per_quadrant or max_points_per_quadrant "
3992 : "is >= 1");
3993 : }
3994 :
3995 : static const char *const apszKnownOptions[] = {
3996 : "radius",
3997 : "radius1",
3998 : "radius2",
3999 : "angle",
4000 : "min_points",
4001 : "max_points",
4002 : "nodata",
4003 : "min_points_per_quadrant",
4004 : "max_points_per_quadrant",
4005 : nullptr};
4006 13 : papszKnownOptions = apszKnownOptions;
4007 :
4008 13 : break;
4009 : }
4010 14 : case GGA_NearestNeighbor:
4011 : {
4012 14 : *ppOptions = CPLMalloc(sizeof(GDALGridNearestNeighborOptions));
4013 :
4014 14 : GDALGridNearestNeighborOptions *const poNeighborOpts =
4015 : static_cast<GDALGridNearestNeighborOptions *>(*ppOptions);
4016 :
4017 14 : poNeighborOpts->nSizeOfStructure = sizeof(*poNeighborOpts);
4018 :
4019 14 : const char *pszValue = CSLFetchNameValue(papszParams, "radius");
4020 14 : if (pszValue)
4021 : {
4022 0 : poNeighborOpts->dfRadius1 = CPLAtofM(pszValue);
4023 0 : poNeighborOpts->dfRadius2 = poNeighborOpts->dfRadius1;
4024 : }
4025 : else
4026 : {
4027 14 : pszValue = CSLFetchNameValue(papszParams, "radius1");
4028 14 : poNeighborOpts->dfRadius1 = pszValue ? CPLAtofM(pszValue) : 0.0;
4029 :
4030 14 : pszValue = CSLFetchNameValue(papszParams, "radius2");
4031 14 : poNeighborOpts->dfRadius2 = pszValue ? CPLAtofM(pszValue) : 0.0;
4032 : }
4033 :
4034 14 : pszValue = CSLFetchNameValue(papszParams, "angle");
4035 14 : poNeighborOpts->dfAngle = pszValue ? CPLAtofM(pszValue) : 0.0;
4036 :
4037 14 : pszValue = CSLFetchNameValue(papszParams, "nodata");
4038 14 : poNeighborOpts->dfNoDataValue = pszValue ? CPLAtofM(pszValue) : 0.0;
4039 :
4040 : static const char *const apszKnownOptions[] = {
4041 : "radius", "radius1", "radius2", "angle", "nodata", nullptr};
4042 14 : papszKnownOptions = apszKnownOptions;
4043 :
4044 14 : break;
4045 : }
4046 51 : case GGA_MetricMinimum:
4047 : case GGA_MetricMaximum:
4048 : case GGA_MetricRange:
4049 : case GGA_MetricCount:
4050 : case GGA_MetricAverageDistance:
4051 : case GGA_MetricAverageDistancePts:
4052 : {
4053 51 : *ppOptions = CPLMalloc(sizeof(GDALGridDataMetricsOptions));
4054 :
4055 51 : GDALGridDataMetricsOptions *const poMetricsOptions =
4056 : static_cast<GDALGridDataMetricsOptions *>(*ppOptions);
4057 :
4058 51 : poMetricsOptions->nSizeOfStructure = sizeof(*poMetricsOptions);
4059 :
4060 51 : const char *pszValue = CSLFetchNameValue(papszParams, "radius");
4061 51 : if (pszValue)
4062 : {
4063 25 : poMetricsOptions->dfRadius1 = CPLAtofM(pszValue);
4064 25 : poMetricsOptions->dfRadius2 = poMetricsOptions->dfRadius1;
4065 : }
4066 : else
4067 : {
4068 26 : pszValue = CSLFetchNameValue(papszParams, "radius1");
4069 26 : poMetricsOptions->dfRadius1 =
4070 26 : pszValue ? CPLAtofM(pszValue) : 0.0;
4071 :
4072 26 : pszValue = CSLFetchNameValue(papszParams, "radius2");
4073 26 : poMetricsOptions->dfRadius2 =
4074 26 : pszValue ? CPLAtofM(pszValue) : 0.0;
4075 : }
4076 :
4077 51 : pszValue = CSLFetchNameValue(papszParams, "angle");
4078 51 : poMetricsOptions->dfAngle = pszValue ? CPLAtofM(pszValue) : 0.0;
4079 :
4080 51 : pszValue = CSLFetchNameValue(papszParams, "min_points");
4081 51 : poMetricsOptions->nMinPoints = pszValue ? atoi(pszValue) : 0;
4082 :
4083 51 : pszValue = CSLFetchNameValue(papszParams, "nodata");
4084 51 : poMetricsOptions->dfNoDataValue =
4085 51 : pszValue ? CPLAtofM(pszValue) : 0.0;
4086 :
4087 : pszValue =
4088 51 : CSLFetchNameValue(papszParams, "min_points_per_quadrant");
4089 51 : poMetricsOptions->nMinPointsPerQuadrant =
4090 51 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
4091 :
4092 : pszValue =
4093 51 : CSLFetchNameValue(papszParams, "max_points_per_quadrant");
4094 51 : poMetricsOptions->nMaxPointsPerQuadrant =
4095 51 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
4096 :
4097 51 : if (poMetricsOptions->nMinPointsPerQuadrant != 0 ||
4098 26 : poMetricsOptions->nMaxPointsPerQuadrant != 0)
4099 : {
4100 25 : if (*peAlgorithm == GGA_MetricAverageDistancePts)
4101 : {
4102 0 : CPLError(CE_Failure, CPLE_NotSupported,
4103 : "Algorithm %s not supported when "
4104 : "per quadrant parameters are specified",
4105 : szAlgNameAverageDistancePts);
4106 0 : CSLDestroy(papszParams);
4107 0 : return CE_Failure;
4108 : }
4109 25 : if (!(poMetricsOptions->dfRadius1 > 0) ||
4110 25 : !(poMetricsOptions->dfRadius2 > 0))
4111 : {
4112 0 : CPLError(CE_Failure, CPLE_IllegalArg,
4113 : "Radius value should be strictly positive when "
4114 : "per quadrant parameters are specified");
4115 0 : CSLDestroy(papszParams);
4116 0 : return CE_Failure;
4117 : }
4118 25 : if (poMetricsOptions->dfAngle != 0)
4119 : {
4120 0 : CPLError(CE_Failure, CPLE_NotSupported,
4121 : "angle != 0 not supported when "
4122 : "per quadrant parameters are specified");
4123 0 : CSLDestroy(papszParams);
4124 0 : return CE_Failure;
4125 : }
4126 : }
4127 :
4128 : static const char *const apszKnownOptions[] = {
4129 : "radius",
4130 : "radius1",
4131 : "radius2",
4132 : "angle",
4133 : "min_points",
4134 : "nodata",
4135 : "min_points_per_quadrant",
4136 : "max_points_per_quadrant",
4137 : nullptr};
4138 51 : papszKnownOptions = apszKnownOptions;
4139 :
4140 51 : break;
4141 : }
4142 2 : case GGA_Linear:
4143 : {
4144 2 : *ppOptions = CPLMalloc(sizeof(GDALGridLinearOptions));
4145 :
4146 2 : GDALGridLinearOptions *const poLinearOpts =
4147 : static_cast<GDALGridLinearOptions *>(*ppOptions);
4148 :
4149 2 : poLinearOpts->nSizeOfStructure = sizeof(*poLinearOpts);
4150 :
4151 2 : const char *pszValue = CSLFetchNameValue(papszParams, "radius");
4152 2 : poLinearOpts->dfRadius = pszValue ? CPLAtofM(pszValue) : -1.0;
4153 :
4154 2 : pszValue = CSLFetchNameValue(papszParams, "nodata");
4155 2 : poLinearOpts->dfNoDataValue = pszValue ? CPLAtofM(pszValue) : 0.0;
4156 :
4157 : static const char *const apszKnownOptions[] = {"radius", "nodata",
4158 : nullptr};
4159 2 : papszKnownOptions = apszKnownOptions;
4160 :
4161 2 : break;
4162 : }
4163 : }
4164 :
4165 228 : if (papszKnownOptions)
4166 : {
4167 624 : for (int i = 1; papszParams[i] != nullptr; ++i)
4168 : {
4169 396 : char *pszKey = nullptr;
4170 396 : CPL_IGNORE_RET_VAL(CPLParseNameValue(papszParams[i], &pszKey));
4171 396 : if (pszKey)
4172 : {
4173 396 : bool bKnownKey = false;
4174 1677 : for (const char *const *papszKnownKeyIter = papszKnownOptions;
4175 1677 : *papszKnownKeyIter; ++papszKnownKeyIter)
4176 : {
4177 1677 : if (EQUAL(*papszKnownKeyIter, pszKey))
4178 : {
4179 396 : bKnownKey = true;
4180 396 : break;
4181 : }
4182 : }
4183 396 : if (!bKnownKey)
4184 : {
4185 0 : CPLError(CE_Warning, CPLE_AppDefined, "Option %s ignored",
4186 : pszKey);
4187 : }
4188 : }
4189 396 : CPLFree(pszKey);
4190 : }
4191 : }
4192 :
4193 228 : CSLDestroy(papszParams);
4194 228 : return CE_None;
4195 : }
|
