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 5895980 : static void GDALGridGetPointBounds(const void *hFeature, CPLRectObj *pBounds)
47 : {
48 5895980 : const GDALGridPoint *psPoint = static_cast<const GDALGridPoint *>(hFeature);
49 5895980 : GDALGridXYArrays *psXYArrays = psPoint->psXYArrays;
50 5895980 : const int i = psPoint->i;
51 5895980 : const double dfX = psXYArrays->padfX[i];
52 5895980 : const double dfY = psXYArrays->padfY[i];
53 5895980 : pBounds->minx = dfX;
54 5895980 : pBounds->miny = dfY;
55 5895980 : pBounds->maxx = dfX;
56 5895980 : pBounds->maxy = dfY;
57 5895980 : }
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 433756 : 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 433756 : const GDALGridInverseDistanceToAPowerOptions *const poOptions =
119 : static_cast<const GDALGridInverseDistanceToAPowerOptions *>(
120 : poOptionsIn);
121 :
122 : // Pre-compute search ellipse parameters.
123 433756 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
124 433756 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
125 433756 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
126 :
127 : // Compute coefficients for coordinate system rotation.
128 433756 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
129 433756 : const bool bRotated = dfAngle != 0.0;
130 433756 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
131 433756 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
132 :
133 433756 : const double dfPowerDiv2 = poOptions->dfPower / 2;
134 433756 : const double dfSmoothing = poOptions->dfSmoothing;
135 433756 : const GUInt32 nMaxPoints = poOptions->nMaxPoints;
136 433756 : double dfNominator = 0.0;
137 433756 : double dfDenominator = 0.0;
138 433756 : GUInt32 n = 0;
139 :
140 1840480 : for (GUInt32 i = 0; i < nPoints; i++)
141 : {
142 1458870 : double dfRX = padfX[i] - dfXPoint;
143 1458870 : double dfRY = padfY[i] - dfYPoint;
144 1458870 : const double dfR2 =
145 1458870 : dfRX * dfRX + dfRY * dfRY + dfSmoothing * dfSmoothing;
146 :
147 1458870 : if (bRotated)
148 : {
149 168479 : const double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
150 168479 : const double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
151 :
152 168479 : dfRX = dfRXRotated;
153 168479 : dfRY = dfRYRotated;
154 : }
155 :
156 : // Is this point located inside the search ellipse?
157 1458870 : 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 960582 : if (dfR2 < 0.0000000000001)
163 : {
164 0 : *pdfValue = padfZ[i];
165 0 : return CE_None;
166 : }
167 :
168 960582 : const double dfW = pow(dfR2, dfPowerDiv2);
169 960582 : const double dfInvW = 1.0 / dfW;
170 960582 : dfNominator += dfInvW * padfZ[i];
171 960582 : dfDenominator += dfInvW;
172 960582 : n++;
173 960582 : if (nMaxPoints > 0 && n > nMaxPoints)
174 52152 : break;
175 : }
176 : }
177 :
178 433756 : if (n < poOptions->nMinPoints || dfDenominator == 0.0)
179 : {
180 9916 : *pdfValue = poOptions->dfNoDataValue;
181 : }
182 : else
183 : {
184 423840 : *pdfValue = dfNominator / dfDenominator;
185 : }
186 :
187 433756 : 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 453663 : 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 453663 : CPL_IGNORE_RET_VAL(nPoints);
247 :
248 : const GDALGridInverseDistanceToAPowerNearestNeighborOptions
249 449235 : *const poOptions = static_cast<
250 : const GDALGridInverseDistanceToAPowerNearestNeighborOptions *>(
251 : poOptionsIn);
252 449235 : const double dfRadius = poOptions->dfRadius;
253 449235 : const double dfSmoothing = poOptions->dfSmoothing;
254 449235 : const double dfSmoothing2 = dfSmoothing * dfSmoothing;
255 :
256 449235 : const GUInt32 nMaxPoints = poOptions->nMaxPoints;
257 :
258 449235 : GDALGridExtraParameters *psExtraParams =
259 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
260 449235 : const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
261 449235 : CPLAssert(phQuadTree);
262 :
263 449235 : const double dfRPower2 = psExtraParams->dfRadiusPower2PreComp;
264 449235 : const double dfPowerDiv2 = psExtraParams->dfPowerDiv2PreComp;
265 :
266 899889 : std::multimap<double, double> oMapDistanceToZValues;
267 :
268 445996 : const double dfSearchRadius = dfRadius;
269 : CPLRectObj sAoi;
270 445996 : sAoi.minx = dfXPoint - dfSearchRadius;
271 445996 : sAoi.miny = dfYPoint - dfSearchRadius;
272 445996 : sAoi.maxx = dfXPoint + dfSearchRadius;
273 445996 : sAoi.maxy = dfYPoint + dfSearchRadius;
274 445996 : int nFeatureCount = 0;
275 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
276 445996 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
277 450835 : if (nFeatureCount != 0)
278 : {
279 2600840 : for (int k = 0; k < nFeatureCount; k++)
280 : {
281 2202470 : const int i = papsPoints[k]->i;
282 2202470 : const double dfRX = padfX[i] - dfXPoint;
283 2202470 : const double dfRY = padfY[i] - dfYPoint;
284 :
285 2202470 : const double dfR2 = dfRX * dfRX + dfRY * dfRY;
286 : // real distance + smoothing
287 2202470 : const double dfRsmoothed2 = dfR2 + dfSmoothing2;
288 2202470 : 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 2202470 : if (dfR2 <= dfRPower2)
296 : {
297 : oMapDistanceToZValues.insert(
298 2180340 : std::make_pair(dfRsmoothed2, padfZ[i]));
299 : }
300 : }
301 : }
302 399248 : CPLFree(papsPoints);
303 :
304 450991 : double dfNominator = 0.0;
305 450991 : double dfDenominator = 0.0;
306 450991 : 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 1818010 : for (std::multimap<double, double>::iterator oMapDistanceToZValuesIter =
312 450991 : oMapDistanceToZValues.begin();
313 2259420 : oMapDistanceToZValuesIter != oMapDistanceToZValues.end();
314 1804590 : ++oMapDistanceToZValuesIter)
315 : {
316 1864750 : const double dfR2 = oMapDistanceToZValuesIter->first;
317 1856760 : const double dfZ = oMapDistanceToZValuesIter->second;
318 :
319 1869490 : const double dfW = pow(dfR2, dfPowerDiv2);
320 1869490 : const double dfInvW = 1.0 / dfW;
321 1869490 : dfNominator += dfInvW * dfZ;
322 1869490 : dfDenominator += dfInvW;
323 1869490 : n++;
324 1869490 : if (nMaxPoints > 0 && n >= nMaxPoints)
325 : {
326 64897 : break;
327 : }
328 : }
329 :
330 450654 : if (n < poOptions->nMinPoints || dfDenominator == 0.0)
331 : {
332 64192 : *pdfValue = poOptions->dfNoDataValue;
333 : }
334 : else
335 : {
336 386462 : *pdfValue = dfNominator / dfDenominator;
337 : }
338 :
339 450654 : 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 258131 : 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 258131 : *const poOptions = static_cast<
357 : const GDALGridInverseDistanceToAPowerNearestNeighborOptions *>(
358 : poOptionsIn);
359 258131 : const double dfRadius = poOptions->dfRadius;
360 258131 : const double dfSmoothing = poOptions->dfSmoothing;
361 258131 : const double dfSmoothing2 = dfSmoothing * dfSmoothing;
362 :
363 258131 : const GUInt32 nMaxPoints = poOptions->nMaxPoints;
364 258131 : const GUInt32 nMinPointsPerQuadrant = poOptions->nMinPointsPerQuadrant;
365 258131 : const GUInt32 nMaxPointsPerQuadrant = poOptions->nMaxPointsPerQuadrant;
366 :
367 258131 : GDALGridExtraParameters *psExtraParams =
368 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
369 258131 : const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
370 258131 : CPLAssert(phQuadTree);
371 :
372 258131 : const double dfRPower2 = psExtraParams->dfRadiusPower2PreComp;
373 258131 : const double dfPowerDiv2 = psExtraParams->dfPowerDiv2PreComp;
374 2519130 : std::multimap<double, double> oMapDistanceToZValuesPerQuadrant[4];
375 :
376 257352 : const double dfSearchRadius = dfRadius;
377 : CPLRectObj sAoi;
378 257352 : sAoi.minx = dfXPoint - dfSearchRadius;
379 257352 : sAoi.miny = dfYPoint - dfSearchRadius;
380 257352 : sAoi.maxx = dfXPoint + dfSearchRadius;
381 257352 : sAoi.maxy = dfYPoint + dfSearchRadius;
382 257352 : int nFeatureCount = 0;
383 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
384 257352 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
385 256865 : if (nFeatureCount != 0)
386 : {
387 1454630 : for (int k = 0; k < nFeatureCount; k++)
388 : {
389 1220190 : const int i = papsPoints[k]->i;
390 1220190 : const double dfRX = padfX[i] - dfXPoint;
391 1220190 : const double dfRY = padfY[i] - dfYPoint;
392 :
393 1220190 : const double dfR2 = dfRX * dfRX + dfRY * dfRY;
394 : // real distance + smoothing
395 1220190 : const double dfRsmoothed2 = dfR2 + dfSmoothing2;
396 1220190 : 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 1220190 : if (dfR2 <= dfRPower2)
404 : {
405 1125550 : const int iQuadrant =
406 1125550 : ((dfRX >= 0) ? 1 : 0) | (((dfRY >= 0) ? 1 : 0) << 1);
407 : oMapDistanceToZValuesPerQuadrant[iQuadrant].insert(
408 1125550 : std::make_pair(dfRsmoothed2, padfZ[i]));
409 : }
410 : }
411 : }
412 235028 : CPLFree(papsPoints);
413 :
414 : std::multimap<double, double>::iterator aoIter[] = {
415 257146 : oMapDistanceToZValuesPerQuadrant[0].begin(),
416 251465 : oMapDistanceToZValuesPerQuadrant[1].begin(),
417 250475 : oMapDistanceToZValuesPerQuadrant[2].begin(),
418 251591 : oMapDistanceToZValuesPerQuadrant[3].begin(),
419 257146 : };
420 251751 : constexpr int ALL_QUADRANT_FLAGS = 1 + 2 + 4 + 8;
421 :
422 : // Examine all "neighbors" within the radius (sorted by distance via the
423 : // multimap), and use the closest n points based on distance until the max
424 : // is reached.
425 : // Do that by fetching the nearest point in quadrant 0, then the nearest
426 : // point in quadrant 1, 2 and 3, and starting again with the next nearest
427 : // point in quarant 0, etc.
428 251751 : int nQuadrantIterFinishedFlag = 0;
429 251751 : GUInt32 anPerQuadrant[4] = {0};
430 251751 : double dfNominator = 0.0;
431 251751 : double dfDenominator = 0.0;
432 251751 : GUInt32 n = 0;
433 251751 : for (int iQuadrant = 0; /* true */; iQuadrant = (iQuadrant + 1) % 4)
434 : {
435 2238350 : if (aoIter[iQuadrant] ==
436 5180140 : oMapDistanceToZValuesPerQuadrant[iQuadrant].end() ||
437 559455 : (nMaxPointsPerQuadrant > 0 &&
438 559455 : anPerQuadrant[iQuadrant] >= nMaxPointsPerQuadrant))
439 : {
440 1297860 : nQuadrantIterFinishedFlag |= 1 << iQuadrant;
441 1297860 : if (nQuadrantIterFinishedFlag == ALL_QUADRANT_FLAGS)
442 257496 : break;
443 1040370 : continue;
444 : }
445 :
446 1029330 : const double dfR2 = aoIter[iQuadrant]->first;
447 1005000 : const double dfZ = aoIter[iQuadrant]->second;
448 1008480 : ++aoIter[iQuadrant];
449 :
450 1001380 : const double dfW = pow(dfR2, dfPowerDiv2);
451 1001380 : const double dfInvW = 1.0 / dfW;
452 1001380 : dfNominator += dfInvW * dfZ;
453 1001380 : dfDenominator += dfInvW;
454 1001380 : n++;
455 1001380 : anPerQuadrant[iQuadrant]++;
456 1001380 : if (nMaxPoints > 0 && n >= nMaxPoints)
457 : {
458 1 : break;
459 : }
460 2041740 : }
461 :
462 257497 : if (nMinPointsPerQuadrant > 0 &&
463 128847 : (anPerQuadrant[0] < nMinPointsPerQuadrant ||
464 52152 : anPerQuadrant[1] < nMinPointsPerQuadrant ||
465 42366 : anPerQuadrant[2] < nMinPointsPerQuadrant ||
466 33292 : anPerQuadrant[3] < nMinPointsPerQuadrant))
467 : {
468 100094 : *pdfValue = poOptions->dfNoDataValue;
469 : }
470 157403 : else if (n < poOptions->nMinPoints || dfDenominator == 0.0)
471 : {
472 0 : *pdfValue = poOptions->dfNoDataValue;
473 : }
474 : else
475 : {
476 158097 : *pdfValue = dfNominator / dfDenominator;
477 : }
478 :
479 257497 : 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 103121 : 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 103121 : const GDALGridInverseDistanceToAPowerOptions *const poOptions =
503 : static_cast<const GDALGridInverseDistanceToAPowerOptions *>(
504 : poOptionsIn);
505 103121 : const double dfPowerDiv2 = poOptions->dfPower / 2.0;
506 103121 : const double dfSmoothing = poOptions->dfSmoothing;
507 103121 : const double dfSmoothing2 = dfSmoothing * dfSmoothing;
508 103121 : double dfNominator = 0.0;
509 103121 : double dfDenominator = 0.0;
510 103121 : const bool bPower2 = dfPowerDiv2 == 1.0;
511 :
512 103121 : GUInt32 i = 0; // Used after if.
513 103121 : if (bPower2)
514 : {
515 55639 : if (dfSmoothing2 > 0)
516 : {
517 234244 : for (i = 0; i < nPoints; i++)
518 : {
519 179790 : const double dfRX = padfX[i] - dfXPoint;
520 179790 : const double dfRY = padfY[i] - dfYPoint;
521 179790 : const double dfR2 = dfRX * dfRX + dfRY * dfRY + dfSmoothing2;
522 :
523 179790 : const double dfInvR2 = 1.0 / dfR2;
524 179790 : dfNominator += dfInvR2 * padfZ[i];
525 179790 : dfDenominator += dfInvR2;
526 : }
527 : }
528 : else
529 : {
530 60017 : for (i = 0; i < nPoints; i++)
531 : {
532 59233 : const double dfRX = padfX[i] - dfXPoint;
533 59233 : const double dfRY = padfY[i] - dfYPoint;
534 59233 : 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 59233 : if (dfR2 < 0.0000000000001)
539 : {
540 401 : break;
541 : }
542 :
543 58832 : const double dfInvR2 = 1.0 / dfR2;
544 58832 : dfNominator += dfInvR2 * padfZ[i];
545 58832 : dfDenominator += dfInvR2;
546 : }
547 : }
548 : }
549 : else
550 : {
551 172927 : for (i = 0; i < nPoints; i++)
552 : {
553 125445 : const double dfRX = padfX[i] - dfXPoint;
554 125445 : const double dfRY = padfY[i] - dfYPoint;
555 125445 : 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 125445 : if (dfR2 < 0.0000000000001)
560 : {
561 0 : break;
562 : }
563 :
564 125445 : const double dfW = pow(dfR2, dfPowerDiv2);
565 125445 : const double dfInvW = 1.0 / dfW;
566 125445 : dfNominator += dfInvW * padfZ[i];
567 125445 : dfDenominator += dfInvW;
568 : }
569 : }
570 :
571 103121 : if (i != nPoints)
572 : {
573 401 : *pdfValue = padfZ[i];
574 : }
575 102720 : else if (dfDenominator == 0.0)
576 : {
577 0 : *pdfValue = poOptions->dfNoDataValue;
578 : }
579 : else
580 : {
581 102720 : *pdfValue = dfNominator / dfDenominator;
582 : }
583 :
584 103121 : 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 489556 : 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 489556 : const GDALGridMovingAverageOptions *const poOptions =
639 : static_cast<const GDALGridMovingAverageOptions *>(poOptionsIn);
640 : // Pre-compute search ellipse parameters.
641 489556 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
642 489556 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
643 : const double dfSearchRadius =
644 489556 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
645 525482 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
646 :
647 525482 : GDALGridExtraParameters *psExtraParams =
648 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
649 525482 : const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
650 :
651 : // Compute coefficients for coordinate system rotation.
652 525482 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
653 525482 : const bool bRotated = dfAngle != 0.0;
654 :
655 525482 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
656 525482 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
657 :
658 525482 : double dfAccumulator = 0.0;
659 :
660 525482 : GUInt32 n = 0; // Used after for.
661 525482 : if (phQuadTree != nullptr)
662 : {
663 : CPLRectObj sAoi;
664 1597 : sAoi.minx = dfXPoint - dfSearchRadius;
665 1597 : sAoi.miny = dfYPoint - dfSearchRadius;
666 1597 : sAoi.maxx = dfXPoint + dfSearchRadius;
667 1597 : sAoi.maxy = dfYPoint + dfSearchRadius;
668 1597 : int nFeatureCount = 0;
669 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
670 1597 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
671 1596 : if (nFeatureCount != 0)
672 : {
673 39283 : for (int k = 0; k < nFeatureCount; k++)
674 : {
675 37690 : const int i = papsPoints[k]->i;
676 37690 : const double dfRX = padfX[i] - dfXPoint;
677 37690 : const double dfRY = padfY[i] - dfYPoint;
678 :
679 37690 : if (dfRadius2Square * dfRX * dfRX +
680 37690 : dfRadius1Square * dfRY * dfRY <=
681 : dfR12Square)
682 : {
683 31010 : dfAccumulator += padfZ[i];
684 31010 : n++;
685 : }
686 : }
687 : }
688 1596 : CPLFree(papsPoints);
689 : }
690 : else
691 : {
692 2785640 : for (GUInt32 i = 0; i < nPoints; i++)
693 : {
694 2261750 : double dfRX = padfX[i] - dfXPoint;
695 2261750 : double dfRY = padfY[i] - dfYPoint;
696 :
697 2261750 : if (bRotated)
698 : {
699 337889 : const double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
700 337889 : const double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
701 :
702 337889 : dfRX = dfRXRotated;
703 337889 : dfRY = dfRYRotated;
704 : }
705 :
706 : // Is this point located inside the search ellipse?
707 2261750 : if (dfRadius2Square * dfRX * dfRX + dfRadius1Square * dfRY * dfRY <=
708 : dfR12Square)
709 : {
710 1570320 : dfAccumulator += padfZ[i];
711 1570320 : n++;
712 : }
713 : }
714 : }
715 :
716 525478 : if (n < poOptions->nMinPoints || n == 0)
717 : {
718 80781 : *pdfValue = poOptions->dfNoDataValue;
719 : }
720 : else
721 : {
722 444697 : *pdfValue = dfAccumulator / n;
723 : }
724 :
725 525478 : return CE_None;
726 : }
727 :
728 : /************************************************************************/
729 : /* GDALGridMovingAveragePerQuadrant() */
730 : /************************************************************************/
731 :
732 : /**
733 : * Moving average, with a per-quadrant search logic.
734 : */
735 193875 : 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 193875 : const GDALGridMovingAverageOptions *const poOptions =
741 : static_cast<const GDALGridMovingAverageOptions *>(poOptionsIn);
742 193875 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
743 193875 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
744 193875 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
745 :
746 193875 : const GUInt32 nMaxPoints = poOptions->nMaxPoints;
747 193875 : const GUInt32 nMinPointsPerQuadrant = poOptions->nMinPointsPerQuadrant;
748 193875 : const GUInt32 nMaxPointsPerQuadrant = poOptions->nMaxPointsPerQuadrant;
749 :
750 193875 : GDALGridExtraParameters *psExtraParams =
751 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
752 193875 : const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
753 193875 : CPLAssert(phQuadTree);
754 :
755 1883600 : std::multimap<double, double> oMapDistanceToZValuesPerQuadrant[4];
756 :
757 : const double dfSearchRadius =
758 191796 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
759 : CPLRectObj sAoi;
760 189404 : sAoi.minx = dfXPoint - dfSearchRadius;
761 189404 : sAoi.miny = dfYPoint - dfSearchRadius;
762 189404 : sAoi.maxx = dfXPoint + dfSearchRadius;
763 189404 : sAoi.maxy = dfYPoint + dfSearchRadius;
764 189404 : int nFeatureCount = 0;
765 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
766 189404 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
767 192760 : if (nFeatureCount != 0)
768 : {
769 1104720 : for (int k = 0; k < nFeatureCount; k++)
770 : {
771 936692 : const int i = papsPoints[k]->i;
772 936692 : const double dfRX = padfX[i] - dfXPoint;
773 936692 : const double dfRY = padfY[i] - dfYPoint;
774 936692 : const double dfRXSquare = dfRX * dfRX;
775 936692 : const double dfRYSquare = dfRY * dfRY;
776 :
777 936692 : if (dfRadius2Square * dfRXSquare + dfRadius1Square * dfRYSquare <=
778 : dfR12Square)
779 : {
780 945169 : const int iQuadrant =
781 945169 : ((dfRX >= 0) ? 1 : 0) | (((dfRY >= 0) ? 1 : 0) << 1);
782 : oMapDistanceToZValuesPerQuadrant[iQuadrant].insert(
783 945169 : std::make_pair(dfRXSquare + dfRYSquare, padfZ[i]));
784 : }
785 : }
786 : }
787 168887 : CPLFree(papsPoints);
788 :
789 : std::multimap<double, double>::iterator aoIter[] = {
790 192252 : oMapDistanceToZValuesPerQuadrant[0].begin(),
791 187338 : oMapDistanceToZValuesPerQuadrant[1].begin(),
792 188803 : oMapDistanceToZValuesPerQuadrant[2].begin(),
793 189366 : oMapDistanceToZValuesPerQuadrant[3].begin(),
794 192252 : };
795 187547 : constexpr int ALL_QUADRANT_FLAGS = 1 + 2 + 4 + 8;
796 :
797 : // Examine all "neighbors" within the radius (sorted by distance via the
798 : // multimap), and use the closest n points based on distance until the max
799 : // is reached.
800 : // Do that by fetching the nearest point in quadrant 0, then the nearest
801 : // point in quadrant 1, 2 and 3, and starting again with the next nearest
802 : // point in quarant 0, etc.
803 187547 : int nQuadrantIterFinishedFlag = 0;
804 187547 : GUInt32 anPerQuadrant[4] = {0};
805 187547 : double dfNominator = 0.0;
806 187547 : GUInt32 n = 0;
807 187547 : for (int iQuadrant = 0; /* true */; iQuadrant = (iQuadrant + 1) % 4)
808 : {
809 921754 : if (aoIter[iQuadrant] ==
810 2160290 : oMapDistanceToZValuesPerQuadrant[iQuadrant].end() ||
811 307934 : (nMaxPointsPerQuadrant > 0 &&
812 307934 : anPerQuadrant[iQuadrant] >= nMaxPointsPerQuadrant))
813 : {
814 244798 : nQuadrantIterFinishedFlag |= 1 << iQuadrant;
815 244798 : if (nQuadrantIterFinishedFlag == ALL_QUADRANT_FLAGS)
816 64236 : break;
817 180562 : continue;
818 : }
819 :
820 682531 : const double dfZ = aoIter[iQuadrant]->second;
821 676657 : ++aoIter[iQuadrant];
822 :
823 684449 : dfNominator += dfZ;
824 684449 : n++;
825 684449 : anPerQuadrant[iQuadrant]++;
826 684449 : if (nMaxPoints > 0 && n >= nMaxPoints)
827 : {
828 127533 : break;
829 : }
830 737478 : }
831 :
832 191769 : if (nMinPointsPerQuadrant > 0 &&
833 127639 : (anPerQuadrant[0] < nMinPointsPerQuadrant ||
834 127728 : anPerQuadrant[1] < nMinPointsPerQuadrant ||
835 127745 : anPerQuadrant[2] < nMinPointsPerQuadrant ||
836 127821 : anPerQuadrant[3] < nMinPointsPerQuadrant))
837 : {
838 63175 : *pdfValue = poOptions->dfNoDataValue;
839 : }
840 128594 : else if (n < poOptions->nMinPoints || n == 0)
841 : {
842 0 : *pdfValue = poOptions->dfNoDataValue;
843 : }
844 : else
845 : {
846 128796 : *pdfValue = dfNominator / n;
847 : }
848 :
849 384116 : 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 448690 : 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 448690 : const GDALGridNearestNeighborOptions *const poOptions =
889 : static_cast<const GDALGridNearestNeighborOptions *>(poOptionsIn);
890 : // Pre-compute search ellipse parameters.
891 448690 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
892 448690 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
893 448690 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
894 448690 : GDALGridExtraParameters *psExtraParams =
895 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
896 448690 : CPLQuadTree *hQuadTree = psExtraParams->hQuadTree;
897 :
898 : // Compute coefficients for coordinate system rotation.
899 448690 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
900 448690 : const bool bRotated = dfAngle != 0.0;
901 448690 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
902 448690 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
903 :
904 : // If the nearest point will not be found, its value remains as NODATA.
905 448690 : double dfNearestValue = poOptions->dfNoDataValue;
906 448690 : GUInt32 i = 0;
907 :
908 448690 : double dfSearchRadius = psExtraParams->dfInitialSearchRadius;
909 448690 : if (hQuadTree != nullptr)
910 : {
911 129203 : if (poOptions->dfRadius1 > 0 || poOptions->dfRadius2 > 0)
912 62243 : dfSearchRadius =
913 56215 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
914 : CPLRectObj sAoi;
915 323724 : while (dfSearchRadius > 0)
916 : {
917 313537 : sAoi.minx = dfXPoint - dfSearchRadius;
918 313537 : sAoi.miny = dfYPoint - dfSearchRadius;
919 313537 : sAoi.maxx = dfXPoint + dfSearchRadius;
920 313537 : sAoi.maxy = dfYPoint + dfSearchRadius;
921 313537 : int nFeatureCount = 0;
922 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
923 313537 : CPLQuadTreeSearch(hQuadTree, &sAoi, &nFeatureCount));
924 312236 : if (nFeatureCount != 0)
925 : {
926 : // Nearest distance will be initialized with the distance to the
927 : // first point in array.
928 75362 : double dfNearestRSquare = std::numeric_limits<double>::max();
929 378910 : for (int k = 0; k < nFeatureCount; k++)
930 : {
931 303548 : const int idx = papsPoints[k]->i;
932 303548 : const double dfRX = padfX[idx] - dfXPoint;
933 303548 : const double dfRY = padfY[idx] - dfYPoint;
934 :
935 303548 : const double dfR2 = dfRX * dfRX + dfRY * dfRY;
936 303548 : if (dfR2 <= dfNearestRSquare)
937 : {
938 165467 : dfNearestRSquare = dfR2;
939 165467 : dfNearestValue = padfZ[idx];
940 : }
941 : }
942 :
943 75362 : CPLFree(papsPoints);
944 129996 : break;
945 : }
946 :
947 236874 : CPLFree(papsPoints);
948 242713 : if (poOptions->dfRadius1 > 0 || poOptions->dfRadius2 > 0)
949 : break;
950 188493 : 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 319487 : double dfNearestRSquare = std::numeric_limits<double>::max();
960 184343000 : while (i < nPoints)
961 : {
962 184023000 : double dfRX = padfX[i] - dfXPoint;
963 184023000 : double dfRY = padfY[i] - dfYPoint;
964 :
965 184023000 : 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 184023000 : const double dfRXSquare = dfRX * dfRX;
976 184023000 : const double dfRYSquare = dfRY * dfRY;
977 184023000 : if (dfRadius2Square * dfRXSquare + dfRadius1Square * dfRYSquare <=
978 : dfR12Square)
979 : {
980 182670000 : const double dfR2 = dfRXSquare + dfRYSquare;
981 182670000 : if (dfR2 <= dfNearestRSquare)
982 : {
983 16273700 : dfNearestRSquare = dfR2;
984 16273700 : dfNearestValue = padfZ[i];
985 : }
986 : }
987 :
988 184023000 : i++;
989 : }
990 : }
991 :
992 459670 : *pdfValue = dfNearestValue;
993 :
994 459670 : 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 323049 : 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 323049 : const GDALGridDataMetricsOptions *const poOptions =
1043 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
1044 :
1045 : // Pre-compute search ellipse parameters.
1046 323049 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
1047 323049 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
1048 : const double dfSearchRadius =
1049 323049 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
1050 323776 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
1051 :
1052 323776 : GDALGridExtraParameters *psExtraParams =
1053 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
1054 323776 : CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
1055 :
1056 : // Compute coefficients for coordinate system rotation.
1057 323776 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
1058 323776 : const bool bRotated = dfAngle != 0.0;
1059 323776 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
1060 323776 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
1061 :
1062 323776 : double dfMinimumValue = std::numeric_limits<double>::max();
1063 323776 : GUInt32 n = 0;
1064 323776 : if (phQuadTree != nullptr)
1065 : {
1066 : CPLRectObj sAoi;
1067 1589 : sAoi.minx = dfXPoint - dfSearchRadius;
1068 1589 : sAoi.miny = dfYPoint - dfSearchRadius;
1069 1589 : sAoi.maxx = dfXPoint + dfSearchRadius;
1070 1589 : sAoi.maxy = dfYPoint + dfSearchRadius;
1071 1589 : int nFeatureCount = 0;
1072 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
1073 1589 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
1074 1598 : if (nFeatureCount != 0)
1075 : {
1076 33905 : for (int k = 0; k < nFeatureCount; k++)
1077 : {
1078 32309 : const int i = papsPoints[k]->i;
1079 32309 : const double dfRX = padfX[i] - dfXPoint;
1080 32309 : const double dfRY = padfY[i] - dfYPoint;
1081 :
1082 32309 : if (dfRadius2Square * dfRX * dfRX +
1083 32309 : dfRadius1Square * dfRY * dfRY <=
1084 : dfR12Square)
1085 : {
1086 20575 : if (dfMinimumValue > padfZ[i])
1087 3256 : dfMinimumValue = padfZ[i];
1088 20575 : n++;
1089 : }
1090 : }
1091 : }
1092 1598 : CPLFree(papsPoints);
1093 : }
1094 : else
1095 : {
1096 322187 : GUInt32 i = 0;
1097 2024060 : while (i < nPoints)
1098 : {
1099 1701870 : double dfRX = padfX[i] - dfXPoint;
1100 1701870 : double dfRY = padfY[i] - dfYPoint;
1101 :
1102 1701870 : if (bRotated)
1103 : {
1104 71792 : const double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
1105 71792 : const double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
1106 :
1107 71792 : dfRX = dfRXRotated;
1108 71792 : dfRY = dfRYRotated;
1109 : }
1110 :
1111 : // Is this point located inside the search ellipse?
1112 1701870 : if (dfRadius2Square * dfRX * dfRX + dfRadius1Square * dfRY * dfRY <=
1113 : dfR12Square)
1114 : {
1115 682049 : if (dfMinimumValue > padfZ[i])
1116 411357 : dfMinimumValue = padfZ[i];
1117 682049 : n++;
1118 : }
1119 :
1120 1701870 : i++;
1121 : }
1122 : }
1123 :
1124 323781 : if (n < poOptions->nMinPoints || n == 0)
1125 : {
1126 77508 : *pdfValue = poOptions->dfNoDataValue;
1127 : }
1128 : else
1129 : {
1130 246273 : *pdfValue = dfMinimumValue;
1131 : }
1132 :
1133 323781 : 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 129493 : 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 129493 : const GDALGridDataMetricsOptions *const poOptions =
1151 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
1152 :
1153 : // Pre-compute search ellipse parameters.
1154 129493 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
1155 129493 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
1156 127385 : const double dfSearchRadius =
1157 129493 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
1158 127385 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
1159 :
1160 : // const GUInt32 nMaxPoints = poOptions->nMaxPoints;
1161 127385 : const GUInt32 nMinPointsPerQuadrant = poOptions->nMinPointsPerQuadrant;
1162 127385 : const GUInt32 nMaxPointsPerQuadrant = poOptions->nMaxPointsPerQuadrant;
1163 :
1164 127385 : GDALGridExtraParameters *psExtraParams =
1165 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
1166 127385 : const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
1167 127385 : CPLAssert(phQuadTree);
1168 :
1169 : CPLRectObj sAoi;
1170 127385 : sAoi.minx = dfXPoint - dfSearchRadius;
1171 127385 : sAoi.miny = dfYPoint - dfSearchRadius;
1172 127385 : sAoi.maxx = dfXPoint + dfSearchRadius;
1173 127385 : sAoi.maxy = dfYPoint + dfSearchRadius;
1174 127385 : int nFeatureCount = 0;
1175 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
1176 127385 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
1177 1260000 : std::multimap<double, double> oMapDistanceToZValuesPerQuadrant[4];
1178 :
1179 127867 : if (nFeatureCount != 0)
1180 : {
1181 520048 : for (int k = 0; k < nFeatureCount; k++)
1182 : {
1183 398500 : const int i = papsPoints[k]->i;
1184 398500 : const double dfRX = padfX[i] - dfXPoint;
1185 398500 : const double dfRY = padfY[i] - dfYPoint;
1186 398500 : const double dfRXSquare = dfRX * dfRX;
1187 398500 : const double dfRYSquare = dfRY * dfRY;
1188 :
1189 398500 : if (dfRadius2Square * dfRXSquare + dfRadius1Square * dfRYSquare <=
1190 : dfR12Square)
1191 : {
1192 384557 : const int iQuadrant =
1193 384557 : ((dfRX >= 0) ? 1 : 0) | (((dfRY >= 0) ? 1 : 0) << 1);
1194 378890 : oMapDistanceToZValuesPerQuadrant[iQuadrant].insert(
1195 762066 : std::make_pair(dfRXSquare + dfRYSquare, padfZ[i]));
1196 : }
1197 : }
1198 : }
1199 120819 : CPLFree(papsPoints);
1200 :
1201 629243 : std::multimap<double, double>::iterator aoIter[] = {
1202 126983 : oMapDistanceToZValuesPerQuadrant[0].begin(),
1203 124806 : oMapDistanceToZValuesPerQuadrant[1].begin(),
1204 125309 : oMapDistanceToZValuesPerQuadrant[2].begin(),
1205 125945 : oMapDistanceToZValuesPerQuadrant[3].begin(),
1206 : };
1207 126200 : constexpr int ALL_QUADRANT_FLAGS = 1 + 2 + 4 + 8;
1208 :
1209 : // Examine all "neighbors" within the radius (sorted by distance via the
1210 : // multimap), and use the closest n points based on distance until the max
1211 : // is reached.
1212 : // Do that by fetching the nearest point in quadrant 0, then the nearest
1213 : // point in quadrant 1, 2 and 3, and starting again with the next nearest
1214 : // point in quarant 0, etc.
1215 126200 : int nQuadrantIterFinishedFlag = 0;
1216 126200 : GUInt32 anPerQuadrant[4] = {0};
1217 126200 : double dfExtremum = IS_MIN ? std::numeric_limits<double>::max()
1218 : : -std::numeric_limits<double>::max();
1219 126200 : GUInt32 n = 0;
1220 907781 : for (int iQuadrant = 0; /* true */; iQuadrant = (iQuadrant + 1) % 4)
1221 : {
1222 882902 : if (aoIter[iQuadrant] ==
1223 2126910 : oMapDistanceToZValuesPerQuadrant[iQuadrant].end() ||
1224 309127 : (nMaxPointsPerQuadrant > 0 &&
1225 309127 : anPerQuadrant[iQuadrant] >= nMaxPointsPerQuadrant))
1226 : {
1227 593205 : nQuadrantIterFinishedFlag |= 1 << iQuadrant;
1228 593205 : if (nQuadrantIterFinishedFlag == ALL_QUADRANT_FLAGS)
1229 128444 : break;
1230 464761 : continue;
1231 : }
1232 :
1233 316794 : const double dfZ = aoIter[iQuadrant]->second;
1234 315079 : ++aoIter[iQuadrant];
1235 :
1236 : if (IS_MIN)
1237 : {
1238 316804 : if (dfExtremum > dfZ)
1239 243231 : dfExtremum = dfZ;
1240 : }
1241 : else
1242 : {
1243 16 : if (dfExtremum < dfZ)
1244 6 : dfExtremum = dfZ;
1245 : }
1246 316820 : n++;
1247 316820 : anPerQuadrant[iQuadrant]++;
1248 : /*if( nMaxPoints > 0 && n >= nMaxPoints )
1249 : {
1250 : break;
1251 : }*/
1252 : }
1253 :
1254 128444 : if (nMinPointsPerQuadrant > 0 &&
1255 64190 : (anPerQuadrant[0] < nMinPointsPerQuadrant ||
1256 8 : anPerQuadrant[1] < nMinPointsPerQuadrant ||
1257 6 : anPerQuadrant[2] < nMinPointsPerQuadrant ||
1258 6 : anPerQuadrant[3] < nMinPointsPerQuadrant))
1259 : {
1260 64449 : *pdfValue = poOptions->dfNoDataValue;
1261 : }
1262 63995 : else if (n < poOptions->nMinPoints || n == 0)
1263 : {
1264 1 : *pdfValue = poOptions->dfNoDataValue;
1265 : }
1266 : else
1267 : {
1268 64388 : *pdfValue = dfExtremum;
1269 : }
1270 :
1271 256557 : 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 128636 : 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 128636 : return GDALGridDataMetricMinimumOrMaximumPerQuadrant</*IS_MIN=*/true>(
1287 : poOptionsIn, padfX, padfY, padfZ, dfXPoint, dfYPoint, pdfValue,
1288 129319 : 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 59320 : 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 59320 : const GDALGridDataMetricsOptions *const poOptions =
1337 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
1338 :
1339 : // Pre-compute search ellipse parameters.
1340 59320 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
1341 59320 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
1342 : const double dfSearchRadius =
1343 59320 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
1344 60430 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
1345 :
1346 60430 : GDALGridExtraParameters *psExtraParams =
1347 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
1348 60430 : CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
1349 :
1350 : // Compute coefficients for coordinate system rotation.
1351 60430 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
1352 60430 : const bool bRotated = dfAngle != 0.0;
1353 60430 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
1354 60430 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
1355 :
1356 60430 : double dfMaximumValue = -std::numeric_limits<double>::max();
1357 60430 : GUInt32 n = 0;
1358 60430 : if (phQuadTree != nullptr)
1359 : {
1360 : CPLRectObj sAoi;
1361 1195 : sAoi.minx = dfXPoint - dfSearchRadius;
1362 1195 : sAoi.miny = dfYPoint - dfSearchRadius;
1363 1195 : sAoi.maxx = dfXPoint + dfSearchRadius;
1364 1195 : sAoi.maxy = dfYPoint + dfSearchRadius;
1365 1195 : int nFeatureCount = 0;
1366 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
1367 1195 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
1368 1197 : if (nFeatureCount != 0)
1369 : {
1370 35220 : for (int k = 0; k < nFeatureCount; k++)
1371 : {
1372 34023 : const int i = papsPoints[k]->i;
1373 34023 : const double dfRX = padfX[i] - dfXPoint;
1374 34023 : const double dfRY = padfY[i] - dfYPoint;
1375 :
1376 34023 : if (dfRadius2Square * dfRX * dfRX +
1377 34023 : dfRadius1Square * dfRY * dfRY <=
1378 : dfR12Square)
1379 : {
1380 22870 : if (dfMaximumValue < padfZ[i])
1381 3315 : dfMaximumValue = padfZ[i];
1382 22870 : n++;
1383 : }
1384 : }
1385 : }
1386 1197 : CPLFree(papsPoints);
1387 : }
1388 : else
1389 : {
1390 59235 : GUInt32 i = 0;
1391 486824 : while (i < nPoints)
1392 : {
1393 427589 : double dfRX = padfX[i] - dfXPoint;
1394 427589 : double dfRY = padfY[i] - dfYPoint;
1395 :
1396 427589 : if (bRotated)
1397 : {
1398 146779 : const double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
1399 146779 : const double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
1400 :
1401 146779 : dfRX = dfRXRotated;
1402 146779 : dfRY = dfRYRotated;
1403 : }
1404 :
1405 : // Is this point located inside the search ellipse?
1406 427589 : if (dfRadius2Square * dfRX * dfRX + dfRadius1Square * dfRY * dfRY <=
1407 : dfR12Square)
1408 : {
1409 282619 : if (dfMaximumValue < padfZ[i])
1410 137067 : dfMaximumValue = padfZ[i];
1411 282619 : n++;
1412 : }
1413 :
1414 427589 : i++;
1415 : }
1416 : }
1417 :
1418 60433 : if (n < poOptions->nMinPoints || n == 0)
1419 : {
1420 3535 : *pdfValue = poOptions->dfNoDataValue;
1421 : }
1422 : else
1423 : {
1424 56898 : *pdfValue = dfMaximumValue;
1425 : }
1426 :
1427 60433 : 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 59207 : 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 59207 : const GDALGridDataMetricsOptions *const poOptions =
1494 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
1495 : // Pre-compute search ellipse parameters.
1496 59207 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
1497 59207 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
1498 : const double dfSearchRadius =
1499 59207 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
1500 60279 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
1501 :
1502 60279 : GDALGridExtraParameters *psExtraParams =
1503 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
1504 60279 : CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
1505 :
1506 : // Compute coefficients for coordinate system rotation.
1507 60279 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
1508 60279 : const bool bRotated = dfAngle != 0.0;
1509 60279 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
1510 60279 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
1511 :
1512 60279 : double dfMaximumValue = -std::numeric_limits<double>::max();
1513 60279 : double dfMinimumValue = std::numeric_limits<double>::max();
1514 60279 : GUInt32 n = 0;
1515 60279 : 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 7477 : for (int k = 0; k < nFeatureCount; k++)
1528 : {
1529 6677 : const int i = papsPoints[k]->i;
1530 6677 : const double dfRX = padfX[i] - dfXPoint;
1531 6677 : const double dfRY = padfY[i] - dfYPoint;
1532 :
1533 6677 : if (dfRadius2Square * dfRX * dfRX +
1534 6677 : dfRadius1Square * dfRY * dfRY <=
1535 : dfR12Square)
1536 : {
1537 6666 : if (dfMinimumValue > padfZ[i])
1538 1943 : dfMinimumValue = padfZ[i];
1539 6666 : if (dfMaximumValue < padfZ[i])
1540 1832 : dfMaximumValue = padfZ[i];
1541 6666 : n++;
1542 : }
1543 : }
1544 : }
1545 800 : CPLFree(papsPoints);
1546 : }
1547 : else
1548 : {
1549 59479 : GUInt32 i = 0;
1550 341552 : while (i < nPoints)
1551 : {
1552 282073 : double dfRX = padfX[i] - dfXPoint;
1553 282073 : double dfRY = padfY[i] - dfYPoint;
1554 :
1555 282073 : 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 282073 : if (dfRadius2Square * dfRX * dfRX + dfRadius1Square * dfRY * dfRY <=
1566 : dfR12Square)
1567 : {
1568 278431 : if (dfMinimumValue > padfZ[i])
1569 148714 : dfMinimumValue = padfZ[i];
1570 278431 : if (dfMaximumValue < padfZ[i])
1571 143570 : dfMaximumValue = padfZ[i];
1572 278431 : n++;
1573 : }
1574 :
1575 282073 : i++;
1576 : }
1577 : }
1578 :
1579 60278 : if (n < poOptions->nMinPoints || n == 0)
1580 : {
1581 2172 : *pdfValue = poOptions->dfNoDataValue;
1582 : }
1583 : else
1584 : {
1585 58106 : *pdfValue = dfMaximumValue - dfMinimumValue;
1586 : }
1587 :
1588 60278 : return CE_None;
1589 : }
1590 :
1591 : /************************************************************************/
1592 : /* GDALGridDataMetricRangePerQuadrant() */
1593 : /************************************************************************/
1594 :
1595 : /**
1596 : * Data range (data metric), with a per-quadrant search logic.
1597 : */
1598 5 : static CPLErr GDALGridDataMetricRangePerQuadrant(
1599 : const void *poOptionsIn, GUInt32 /* nPoints */, const double *padfX,
1600 : const double *padfY, const double *padfZ, double dfXPoint, double dfYPoint,
1601 : double *pdfValue, void *hExtraParamsIn)
1602 : {
1603 5 : const GDALGridDataMetricsOptions *const poOptions =
1604 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
1605 :
1606 : // Pre-compute search ellipse parameters.
1607 5 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
1608 5 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
1609 : const double dfSearchRadius =
1610 5 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
1611 5 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
1612 :
1613 : // const GUInt32 nMaxPoints = poOptions->nMaxPoints;
1614 5 : const GUInt32 nMinPointsPerQuadrant = poOptions->nMinPointsPerQuadrant;
1615 5 : const GUInt32 nMaxPointsPerQuadrant = poOptions->nMaxPointsPerQuadrant;
1616 :
1617 5 : GDALGridExtraParameters *psExtraParams =
1618 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
1619 5 : const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
1620 5 : CPLAssert(phQuadTree);
1621 :
1622 : CPLRectObj sAoi;
1623 5 : sAoi.minx = dfXPoint - dfSearchRadius;
1624 5 : sAoi.miny = dfYPoint - dfSearchRadius;
1625 5 : sAoi.maxx = dfXPoint + dfSearchRadius;
1626 5 : sAoi.maxy = dfYPoint + dfSearchRadius;
1627 5 : int nFeatureCount = 0;
1628 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
1629 5 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
1630 50 : std::multimap<double, double> oMapDistanceToZValuesPerQuadrant[4];
1631 :
1632 5 : if (nFeatureCount != 0)
1633 : {
1634 26 : for (int k = 0; k < nFeatureCount; k++)
1635 : {
1636 21 : const int i = papsPoints[k]->i;
1637 21 : const double dfRX = padfX[i] - dfXPoint;
1638 21 : const double dfRY = padfY[i] - dfYPoint;
1639 21 : const double dfRXSquare = dfRX * dfRX;
1640 21 : const double dfRYSquare = dfRY * dfRY;
1641 :
1642 21 : if (dfRadius2Square * dfRXSquare + dfRadius1Square * dfRYSquare <=
1643 : dfR12Square)
1644 : {
1645 17 : const int iQuadrant =
1646 17 : ((dfRX >= 0) ? 1 : 0) | (((dfRY >= 0) ? 1 : 0) << 1);
1647 : oMapDistanceToZValuesPerQuadrant[iQuadrant].insert(
1648 17 : std::make_pair(dfRXSquare + dfRYSquare, padfZ[i]));
1649 : }
1650 : }
1651 : }
1652 5 : CPLFree(papsPoints);
1653 :
1654 : std::multimap<double, double>::iterator aoIter[] = {
1655 5 : oMapDistanceToZValuesPerQuadrant[0].begin(),
1656 5 : oMapDistanceToZValuesPerQuadrant[1].begin(),
1657 5 : oMapDistanceToZValuesPerQuadrant[2].begin(),
1658 5 : oMapDistanceToZValuesPerQuadrant[3].begin(),
1659 5 : };
1660 5 : constexpr int ALL_QUADRANT_FLAGS = 1 + 2 + 4 + 8;
1661 :
1662 : // Examine all "neighbors" within the radius (sorted by distance via the
1663 : // multimap), and use the closest n points based on distance until the max
1664 : // is reached.
1665 : // Do that by fetching the nearest point in quadrant 0, then the nearest
1666 : // point in quadrant 1, 2 and 3, and starting again with the next nearest
1667 : // point in quarant 0, etc.
1668 5 : int nQuadrantIterFinishedFlag = 0;
1669 5 : GUInt32 anPerQuadrant[4] = {0};
1670 5 : double dfMaximumValue = -std::numeric_limits<double>::max();
1671 5 : double dfMinimumValue = std::numeric_limits<double>::max();
1672 5 : GUInt32 n = 0;
1673 5 : for (int iQuadrant = 0; /* true */; iQuadrant = (iQuadrant + 1) % 4)
1674 : {
1675 40 : if (aoIter[iQuadrant] ==
1676 90 : oMapDistanceToZValuesPerQuadrant[iQuadrant].end() ||
1677 10 : (nMaxPointsPerQuadrant > 0 &&
1678 10 : anPerQuadrant[iQuadrant] >= nMaxPointsPerQuadrant))
1679 : {
1680 24 : nQuadrantIterFinishedFlag |= 1 << iQuadrant;
1681 24 : if (nQuadrantIterFinishedFlag == ALL_QUADRANT_FLAGS)
1682 5 : break;
1683 19 : continue;
1684 : }
1685 :
1686 16 : const double dfZ = aoIter[iQuadrant]->second;
1687 16 : ++aoIter[iQuadrant];
1688 :
1689 16 : if (dfMinimumValue > dfZ)
1690 7 : dfMinimumValue = dfZ;
1691 16 : if (dfMaximumValue < dfZ)
1692 6 : dfMaximumValue = dfZ;
1693 16 : n++;
1694 16 : anPerQuadrant[iQuadrant]++;
1695 : /*if( nMaxPoints > 0 && n >= nMaxPoints )
1696 : {
1697 : break;
1698 : }*/
1699 35 : }
1700 :
1701 5 : if (nMinPointsPerQuadrant > 0 &&
1702 5 : (anPerQuadrant[0] < nMinPointsPerQuadrant ||
1703 4 : anPerQuadrant[1] < nMinPointsPerQuadrant ||
1704 3 : anPerQuadrant[2] < nMinPointsPerQuadrant ||
1705 3 : anPerQuadrant[3] < nMinPointsPerQuadrant))
1706 : {
1707 2 : *pdfValue = poOptions->dfNoDataValue;
1708 : }
1709 3 : else if (n < poOptions->nMinPoints || n == 0)
1710 : {
1711 1 : *pdfValue = poOptions->dfNoDataValue;
1712 : }
1713 : else
1714 : {
1715 2 : *pdfValue = dfMaximumValue - dfMinimumValue;
1716 : }
1717 :
1718 10 : return CE_None;
1719 : }
1720 :
1721 : /************************************************************************/
1722 : /* GDALGridDataMetricCount() */
1723 : /************************************************************************/
1724 :
1725 : /**
1726 : * Number of data points (data metric).
1727 : *
1728 : * A number of data points found in grid node search ellipse.
1729 : *
1730 : * \f[
1731 : * Z=n
1732 : * \f]
1733 : *
1734 : * where
1735 : * <ul>
1736 : * <li> \f$Z\f$ is a resulting value at the grid node,
1737 : * <li> \f$n\f$ is a total number of points in search ellipse.
1738 : * </ul>
1739 : *
1740 : * @param poOptionsIn Algorithm parameters. This should point to
1741 : * GDALGridDataMetricsOptions object.
1742 : * @param nPoints Number of elements in input arrays.
1743 : * @param padfX Input array of X coordinates.
1744 : * @param padfY Input array of Y coordinates.
1745 : * @param padfZ Input array of Z values.
1746 : * @param dfXPoint X coordinate of the point to compute.
1747 : * @param dfYPoint Y coordinate of the point to compute.
1748 : * @param pdfValue Pointer to variable where the computed grid node value
1749 : * will be returned.
1750 : * @param hExtraParamsIn extra parameters (unused)
1751 : *
1752 : * @return CE_None on success or CE_Failure if something goes wrong.
1753 : */
1754 :
1755 57622 : 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 57622 : const GDALGridDataMetricsOptions *const poOptions =
1765 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
1766 :
1767 : // Pre-compute search ellipse parameters.
1768 57622 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
1769 57622 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
1770 : const double dfSearchRadius =
1771 57622 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
1772 56980 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
1773 :
1774 56980 : GDALGridExtraParameters *psExtraParams =
1775 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
1776 56980 : CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
1777 :
1778 : // Compute coefficients for coordinate system rotation.
1779 56980 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
1780 56980 : const bool bRotated = dfAngle != 0.0;
1781 56980 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
1782 56980 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
1783 :
1784 56980 : GUInt32 n = 0;
1785 56980 : if (phQuadTree != nullptr)
1786 : {
1787 : CPLRectObj sAoi;
1788 1998 : sAoi.minx = dfXPoint - dfSearchRadius;
1789 1998 : sAoi.miny = dfYPoint - dfSearchRadius;
1790 1998 : sAoi.maxx = dfXPoint + dfSearchRadius;
1791 1998 : sAoi.maxy = dfYPoint + dfSearchRadius;
1792 1998 : int nFeatureCount = 0;
1793 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
1794 1998 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
1795 1995 : if (nFeatureCount != 0)
1796 : {
1797 79255 : for (int k = 0; k < nFeatureCount; k++)
1798 : {
1799 77260 : const int i = papsPoints[k]->i;
1800 77260 : const double dfRX = padfX[i] - dfXPoint;
1801 77260 : const double dfRY = padfY[i] - dfYPoint;
1802 :
1803 77260 : if (dfRadius2Square * dfRX * dfRX +
1804 77260 : dfRadius1Square * dfRY * dfRY <=
1805 : dfR12Square)
1806 : {
1807 54693 : n++;
1808 : }
1809 : }
1810 : }
1811 1995 : CPLFree(papsPoints);
1812 : }
1813 : else
1814 : {
1815 54982 : GUInt32 i = 0;
1816 264726 : while (i < nPoints)
1817 : {
1818 209744 : double dfRX = padfX[i] - dfXPoint;
1819 209744 : double dfRY = padfY[i] - dfYPoint;
1820 :
1821 209744 : 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 209744 : if (dfRadius2Square * dfRX * dfRX + dfRadius1Square * dfRY * dfRY <=
1832 : dfR12Square)
1833 : {
1834 57382 : n++;
1835 : }
1836 :
1837 209744 : i++;
1838 : }
1839 : }
1840 :
1841 56979 : if (n < poOptions->nMinPoints)
1842 : {
1843 0 : *pdfValue = poOptions->dfNoDataValue;
1844 : }
1845 : else
1846 : {
1847 56979 : *pdfValue = static_cast<double>(n);
1848 : }
1849 :
1850 56979 : return CE_None;
1851 : }
1852 :
1853 : /************************************************************************/
1854 : /* GDALGridDataMetricCountPerQuadrant() */
1855 : /************************************************************************/
1856 :
1857 : /**
1858 : * Number of data points (data metric), with a per-quadrant search logic.
1859 : */
1860 5 : static CPLErr GDALGridDataMetricCountPerQuadrant(
1861 : const void *poOptionsIn, GUInt32 /* nPoints */, const double *padfX,
1862 : const double *padfY, const double *padfZ, double dfXPoint, double dfYPoint,
1863 : double *pdfValue, void *hExtraParamsIn)
1864 : {
1865 5 : const GDALGridDataMetricsOptions *const poOptions =
1866 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
1867 :
1868 : // Pre-compute search ellipse parameters.
1869 5 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
1870 5 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
1871 : const double dfSearchRadius =
1872 5 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
1873 5 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
1874 :
1875 : // const GUInt32 nMaxPoints = poOptions->nMaxPoints;
1876 5 : const GUInt32 nMinPointsPerQuadrant = poOptions->nMinPointsPerQuadrant;
1877 5 : const GUInt32 nMaxPointsPerQuadrant = poOptions->nMaxPointsPerQuadrant;
1878 :
1879 5 : GDALGridExtraParameters *psExtraParams =
1880 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
1881 5 : const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
1882 5 : CPLAssert(phQuadTree);
1883 :
1884 : CPLRectObj sAoi;
1885 5 : sAoi.minx = dfXPoint - dfSearchRadius;
1886 5 : sAoi.miny = dfYPoint - dfSearchRadius;
1887 5 : sAoi.maxx = dfXPoint + dfSearchRadius;
1888 5 : sAoi.maxy = dfYPoint + dfSearchRadius;
1889 5 : int nFeatureCount = 0;
1890 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
1891 5 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
1892 50 : std::multimap<double, double> oMapDistanceToZValuesPerQuadrant[4];
1893 :
1894 5 : if (nFeatureCount != 0)
1895 : {
1896 26 : for (int k = 0; k < nFeatureCount; k++)
1897 : {
1898 21 : const int i = papsPoints[k]->i;
1899 21 : const double dfRX = padfX[i] - dfXPoint;
1900 21 : const double dfRY = padfY[i] - dfYPoint;
1901 21 : const double dfRXSquare = dfRX * dfRX;
1902 21 : const double dfRYSquare = dfRY * dfRY;
1903 :
1904 21 : if (dfRadius2Square * dfRXSquare + dfRadius1Square * dfRYSquare <=
1905 : dfR12Square)
1906 : {
1907 17 : const int iQuadrant =
1908 17 : ((dfRX >= 0) ? 1 : 0) | (((dfRY >= 0) ? 1 : 0) << 1);
1909 : oMapDistanceToZValuesPerQuadrant[iQuadrant].insert(
1910 17 : std::make_pair(dfRXSquare + dfRYSquare, padfZ[i]));
1911 : }
1912 : }
1913 : }
1914 5 : CPLFree(papsPoints);
1915 :
1916 : std::multimap<double, double>::iterator aoIter[] = {
1917 5 : oMapDistanceToZValuesPerQuadrant[0].begin(),
1918 5 : oMapDistanceToZValuesPerQuadrant[1].begin(),
1919 5 : oMapDistanceToZValuesPerQuadrant[2].begin(),
1920 5 : oMapDistanceToZValuesPerQuadrant[3].begin(),
1921 5 : };
1922 5 : constexpr int ALL_QUADRANT_FLAGS = 1 + 2 + 4 + 8;
1923 :
1924 : // Examine all "neighbors" within the radius (sorted by distance via the
1925 : // multimap), and use the closest n points based on distance until the max
1926 : // is reached.
1927 : // Do that by fetching the nearest point in quadrant 0, then the nearest
1928 : // point in quadrant 1, 2 and 3, and starting again with the next nearest
1929 : // point in quarant 0, etc.
1930 5 : int nQuadrantIterFinishedFlag = 0;
1931 5 : GUInt32 anPerQuadrant[4] = {0};
1932 5 : GUInt32 n = 0;
1933 40 : for (int iQuadrant = 0; /* true */; iQuadrant = (iQuadrant + 1) % 4)
1934 : {
1935 40 : if (aoIter[iQuadrant] ==
1936 90 : oMapDistanceToZValuesPerQuadrant[iQuadrant].end() ||
1937 10 : (nMaxPointsPerQuadrant > 0 &&
1938 10 : anPerQuadrant[iQuadrant] >= nMaxPointsPerQuadrant))
1939 : {
1940 24 : nQuadrantIterFinishedFlag |= 1 << iQuadrant;
1941 24 : if (nQuadrantIterFinishedFlag == ALL_QUADRANT_FLAGS)
1942 5 : break;
1943 19 : continue;
1944 : }
1945 :
1946 16 : ++aoIter[iQuadrant];
1947 :
1948 16 : n++;
1949 16 : anPerQuadrant[iQuadrant]++;
1950 : /*if( nMaxPoints > 0 && n >= nMaxPoints )
1951 : {
1952 : break;
1953 : }*/
1954 : }
1955 :
1956 5 : if (nMinPointsPerQuadrant > 0 &&
1957 5 : (anPerQuadrant[0] < nMinPointsPerQuadrant ||
1958 4 : anPerQuadrant[1] < nMinPointsPerQuadrant ||
1959 3 : anPerQuadrant[2] < nMinPointsPerQuadrant ||
1960 3 : anPerQuadrant[3] < nMinPointsPerQuadrant))
1961 : {
1962 2 : *pdfValue = poOptions->dfNoDataValue;
1963 : }
1964 3 : else if (n < poOptions->nMinPoints)
1965 : {
1966 1 : *pdfValue = poOptions->dfNoDataValue;
1967 : }
1968 : else
1969 : {
1970 2 : *pdfValue = static_cast<double>(n);
1971 : }
1972 :
1973 10 : return CE_None;
1974 : }
1975 :
1976 : /************************************************************************/
1977 : /* GDALGridDataMetricAverageDistance() */
1978 : /************************************************************************/
1979 :
1980 : /**
1981 : * Average distance (data metric).
1982 : *
1983 : * An average distance between the grid node (center of the search ellipse)
1984 : * and all of the data points found in grid node search ellipse. If there are
1985 : * no points found, the specified NODATA value will be returned.
1986 : *
1987 : * \f[
1988 : * Z=\frac{\sum_{i = 1}^n r_i}{n}
1989 : * \f]
1990 : *
1991 : * where
1992 : * <ul>
1993 : * <li> \f$Z\f$ is a resulting value at the grid node,
1994 : * <li> \f$r_i\f$ is an Euclidean distance from the grid node
1995 : * to point \f$i\f$,
1996 : * <li> \f$n\f$ is a total number of points in search ellipse.
1997 : * </ul>
1998 : *
1999 : * @param poOptionsIn Algorithm parameters. This should point to
2000 : * GDALGridDataMetricsOptions object.
2001 : * @param nPoints Number of elements in input arrays.
2002 : * @param padfX Input array of X coordinates.
2003 : * @param padfY Input array of Y coordinates.
2004 : * @param padfZ Input array of Z values (unused)
2005 : * @param dfXPoint X coordinate of the point to compute.
2006 : * @param dfYPoint Y coordinate of the point to compute.
2007 : * @param pdfValue Pointer to variable where the computed grid node value
2008 : * will be returned.
2009 : * @param hExtraParamsIn extra parameters (unused)
2010 : *
2011 : * @return CE_None on success or CE_Failure if something goes wrong.
2012 : */
2013 :
2014 57913 : 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 57913 : const GDALGridDataMetricsOptions *const poOptions =
2025 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
2026 :
2027 : // Pre-compute search ellipse parameters.
2028 57913 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
2029 57913 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
2030 : const double dfSearchRadius =
2031 57913 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
2032 59253 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
2033 :
2034 59253 : GDALGridExtraParameters *psExtraParams =
2035 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
2036 59253 : CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
2037 :
2038 : // Compute coefficients for coordinate system rotation.
2039 59253 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
2040 59253 : const bool bRotated = dfAngle != 0.0;
2041 59253 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
2042 59253 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
2043 :
2044 59253 : double dfAccumulator = 0.0;
2045 59253 : GUInt32 n = 0;
2046 59253 : if (phQuadTree != nullptr)
2047 : {
2048 : CPLRectObj sAoi;
2049 798 : sAoi.minx = dfXPoint - dfSearchRadius;
2050 798 : sAoi.miny = dfYPoint - dfSearchRadius;
2051 798 : sAoi.maxx = dfXPoint + dfSearchRadius;
2052 798 : sAoi.maxy = dfYPoint + dfSearchRadius;
2053 798 : int nFeatureCount = 0;
2054 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
2055 798 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
2056 795 : if (nFeatureCount != 0)
2057 : {
2058 17723 : for (int k = 0; k < nFeatureCount; k++)
2059 : {
2060 16928 : const int i = papsPoints[k]->i;
2061 16928 : const double dfRX = padfX[i] - dfXPoint;
2062 16928 : const double dfRY = padfY[i] - dfYPoint;
2063 :
2064 16928 : if (dfRadius2Square * dfRX * dfRX +
2065 16928 : dfRadius1Square * dfRY * dfRY <=
2066 : dfR12Square)
2067 : {
2068 14542 : dfAccumulator += sqrt(dfRX * dfRX + dfRY * dfRY);
2069 14542 : n++;
2070 : }
2071 : }
2072 : }
2073 795 : CPLFree(papsPoints);
2074 : }
2075 : else
2076 : {
2077 58455 : GUInt32 i = 0;
2078 :
2079 373440 : while (i < nPoints)
2080 : {
2081 314985 : double dfRX = padfX[i] - dfXPoint;
2082 314985 : double dfRY = padfY[i] - dfYPoint;
2083 :
2084 314985 : 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 314985 : if (dfRadius2Square * dfRX * dfRX + dfRadius1Square * dfRY * dfRY <=
2095 : dfR12Square)
2096 : {
2097 303024 : dfAccumulator += sqrt(dfRX * dfRX + dfRY * dfRY);
2098 303024 : n++;
2099 : }
2100 :
2101 314985 : i++;
2102 : }
2103 : }
2104 :
2105 59252 : if (n < poOptions->nMinPoints || n == 0)
2106 : {
2107 2091 : *pdfValue = poOptions->dfNoDataValue;
2108 : }
2109 : else
2110 : {
2111 57161 : *pdfValue = dfAccumulator / n;
2112 : }
2113 :
2114 59252 : return CE_None;
2115 : }
2116 :
2117 : /************************************************************************/
2118 : /* GDALGridDataMetricAverageDistancePerQuadrant() */
2119 : /************************************************************************/
2120 :
2121 : /**
2122 : * Average distance (data metric), with a per-quadrant search logic.
2123 : */
2124 5 : static CPLErr GDALGridDataMetricAverageDistancePerQuadrant(
2125 : const void *poOptionsIn, GUInt32 /* nPoints */, const double *padfX,
2126 : const double *padfY, const double *padfZ, double dfXPoint, double dfYPoint,
2127 : double *pdfValue, void *hExtraParamsIn)
2128 : {
2129 5 : const GDALGridDataMetricsOptions *const poOptions =
2130 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
2131 :
2132 : // Pre-compute search ellipse parameters.
2133 5 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
2134 5 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
2135 : const double dfSearchRadius =
2136 5 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
2137 5 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
2138 :
2139 : // const GUInt32 nMaxPoints = poOptions->nMaxPoints;
2140 5 : const GUInt32 nMinPointsPerQuadrant = poOptions->nMinPointsPerQuadrant;
2141 5 : const GUInt32 nMaxPointsPerQuadrant = poOptions->nMaxPointsPerQuadrant;
2142 :
2143 5 : GDALGridExtraParameters *psExtraParams =
2144 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
2145 5 : const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
2146 5 : CPLAssert(phQuadTree);
2147 :
2148 : CPLRectObj sAoi;
2149 5 : sAoi.minx = dfXPoint - dfSearchRadius;
2150 5 : sAoi.miny = dfYPoint - dfSearchRadius;
2151 5 : sAoi.maxx = dfXPoint + dfSearchRadius;
2152 5 : sAoi.maxy = dfYPoint + dfSearchRadius;
2153 5 : int nFeatureCount = 0;
2154 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
2155 5 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
2156 50 : std::multimap<double, double> oMapDistanceToZValuesPerQuadrant[4];
2157 :
2158 5 : if (nFeatureCount != 0)
2159 : {
2160 26 : for (int k = 0; k < nFeatureCount; k++)
2161 : {
2162 21 : const int i = papsPoints[k]->i;
2163 21 : const double dfRX = padfX[i] - dfXPoint;
2164 21 : const double dfRY = padfY[i] - dfYPoint;
2165 21 : const double dfRXSquare = dfRX * dfRX;
2166 21 : const double dfRYSquare = dfRY * dfRY;
2167 :
2168 21 : if (dfRadius2Square * dfRXSquare + dfRadius1Square * dfRYSquare <=
2169 : dfR12Square)
2170 : {
2171 17 : const int iQuadrant =
2172 17 : ((dfRX >= 0) ? 1 : 0) | (((dfRY >= 0) ? 1 : 0) << 1);
2173 : oMapDistanceToZValuesPerQuadrant[iQuadrant].insert(
2174 17 : std::make_pair(dfRXSquare + dfRYSquare, padfZ[i]));
2175 : }
2176 : }
2177 : }
2178 5 : CPLFree(papsPoints);
2179 :
2180 : std::multimap<double, double>::iterator aoIter[] = {
2181 5 : oMapDistanceToZValuesPerQuadrant[0].begin(),
2182 5 : oMapDistanceToZValuesPerQuadrant[1].begin(),
2183 5 : oMapDistanceToZValuesPerQuadrant[2].begin(),
2184 5 : oMapDistanceToZValuesPerQuadrant[3].begin(),
2185 5 : };
2186 5 : constexpr int ALL_QUADRANT_FLAGS = 1 + 2 + 4 + 8;
2187 :
2188 : // Examine all "neighbors" within the radius (sorted by distance via the
2189 : // multimap), and use the closest n points based on distance until the max
2190 : // is reached.
2191 : // Do that by fetching the nearest point in quadrant 0, then the nearest
2192 : // point in quadrant 1, 2 and 3, and starting again with the next nearest
2193 : // point in quarant 0, etc.
2194 5 : int nQuadrantIterFinishedFlag = 0;
2195 5 : GUInt32 anPerQuadrant[4] = {0};
2196 5 : GUInt32 n = 0;
2197 5 : double dfAccumulator = 0;
2198 40 : for (int iQuadrant = 0; /* true */; iQuadrant = (iQuadrant + 1) % 4)
2199 : {
2200 40 : if (aoIter[iQuadrant] ==
2201 90 : oMapDistanceToZValuesPerQuadrant[iQuadrant].end() ||
2202 10 : (nMaxPointsPerQuadrant > 0 &&
2203 10 : anPerQuadrant[iQuadrant] >= nMaxPointsPerQuadrant))
2204 : {
2205 24 : nQuadrantIterFinishedFlag |= 1 << iQuadrant;
2206 24 : if (nQuadrantIterFinishedFlag == ALL_QUADRANT_FLAGS)
2207 5 : break;
2208 19 : continue;
2209 : }
2210 :
2211 16 : dfAccumulator += sqrt(aoIter[iQuadrant]->first);
2212 16 : ++aoIter[iQuadrant];
2213 :
2214 16 : n++;
2215 16 : anPerQuadrant[iQuadrant]++;
2216 : /*if( nMaxPoints > 0 && n >= nMaxPoints )
2217 : {
2218 : break;
2219 : }*/
2220 : }
2221 :
2222 5 : if (nMinPointsPerQuadrant > 0 &&
2223 5 : (anPerQuadrant[0] < nMinPointsPerQuadrant ||
2224 4 : anPerQuadrant[1] < nMinPointsPerQuadrant ||
2225 3 : anPerQuadrant[2] < nMinPointsPerQuadrant ||
2226 3 : anPerQuadrant[3] < nMinPointsPerQuadrant))
2227 : {
2228 2 : *pdfValue = poOptions->dfNoDataValue;
2229 : }
2230 3 : else if (n < poOptions->nMinPoints || n == 0)
2231 : {
2232 1 : *pdfValue = poOptions->dfNoDataValue;
2233 : }
2234 : else
2235 : {
2236 2 : *pdfValue = dfAccumulator / n;
2237 : }
2238 :
2239 10 : return CE_None;
2240 : }
2241 :
2242 : /************************************************************************/
2243 : /* GDALGridDataMetricAverageDistancePts() */
2244 : /************************************************************************/
2245 :
2246 : /**
2247 : * Average distance between points (data metric).
2248 : *
2249 : * An average distance between the data points found in grid node search
2250 : * ellipse. The distance between each pair of points within ellipse is
2251 : * calculated and average of all distances is set as a grid node value. If
2252 : * there are no points found, the specified NODATA value will be returned.
2253 :
2254 : *
2255 : * \f[
2256 : * Z=\frac{\sum_{i = 1}^{n-1}\sum_{j=i+1}^{n}
2257 : r_{ij}}{\left(n-1\right)\,n-\frac{n+{\left(n-1\right)}^{2}-1}{2}}
2258 : * \f]
2259 : *
2260 : * where
2261 : * <ul>
2262 : * <li> \f$Z\f$ is a resulting value at the grid node,
2263 : * <li> \f$r_{ij}\f$ is an Euclidean distance between points
2264 : * \f$i\f$ and \f$j\f$,
2265 : * <li> \f$n\f$ is a total number of points in search ellipse.
2266 : * </ul>
2267 : *
2268 : * @param poOptionsIn Algorithm parameters. This should point to
2269 : * GDALGridDataMetricsOptions object.
2270 : * @param nPoints Number of elements in input arrays.
2271 : * @param padfX Input array of X coordinates.
2272 : * @param padfY Input array of Y coordinates.
2273 : * @param padfZ Input array of Z values (unused)
2274 : * @param dfXPoint X coordinate of the point to compute.
2275 : * @param dfYPoint Y coordinate of the point to compute.
2276 : * @param pdfValue Pointer to variable where the computed grid node value
2277 : * will be returned.
2278 : * @param hExtraParamsIn extra parameters (unused)
2279 : *
2280 : * @return CE_None on success or CE_Failure if something goes wrong.
2281 : */
2282 :
2283 62422 : 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 62422 : const GDALGridDataMetricsOptions *const poOptions =
2292 : static_cast<const GDALGridDataMetricsOptions *>(poOptionsIn);
2293 : // Pre-compute search ellipse parameters.
2294 62422 : const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;
2295 62422 : const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;
2296 : const double dfSearchRadius =
2297 62422 : std::max(poOptions->dfRadius1, poOptions->dfRadius2);
2298 63559 : const double dfR12Square = dfRadius1Square * dfRadius2Square;
2299 :
2300 63559 : GDALGridExtraParameters *psExtraParams =
2301 : static_cast<GDALGridExtraParameters *>(hExtraParamsIn);
2302 63559 : CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;
2303 :
2304 : // Compute coefficients for coordinate system rotation.
2305 63559 : const double dfAngle = TO_RADIANS * poOptions->dfAngle;
2306 63559 : const bool bRotated = dfAngle != 0.0;
2307 63559 : const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;
2308 63559 : const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;
2309 :
2310 63559 : double dfAccumulator = 0.0;
2311 63559 : GUInt32 n = 0;
2312 63559 : if (phQuadTree != nullptr)
2313 : {
2314 : CPLRectObj sAoi;
2315 796 : sAoi.minx = dfXPoint - dfSearchRadius;
2316 796 : sAoi.miny = dfYPoint - dfSearchRadius;
2317 796 : sAoi.maxx = dfXPoint + dfSearchRadius;
2318 796 : sAoi.maxy = dfYPoint + dfSearchRadius;
2319 796 : int nFeatureCount = 0;
2320 : GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(
2321 796 : CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));
2322 800 : if (nFeatureCount != 0)
2323 : {
2324 16591 : for (int k = 0; k < nFeatureCount - 1; k++)
2325 : {
2326 15791 : const int i = papsPoints[k]->i;
2327 15791 : const double dfRX1 = padfX[i] - dfXPoint;
2328 15791 : const double dfRY1 = padfY[i] - dfYPoint;
2329 :
2330 15791 : if (dfRadius2Square * dfRX1 * dfRX1 +
2331 15791 : dfRadius1Square * dfRY1 * dfRY1 <=
2332 : dfR12Square)
2333 : {
2334 124727 : 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 110785 : const int ji = papsPoints[j]->i;
2339 110785 : double dfRX2 = padfX[ji] - dfXPoint;
2340 110785 : double dfRY2 = padfY[ji] - dfYPoint;
2341 :
2342 110785 : if (dfRadius2Square * dfRX2 * dfRX2 +
2343 110785 : dfRadius1Square * dfRY2 * dfRY2 <=
2344 : dfR12Square)
2345 : {
2346 104691 : const double dfRX = padfX[ji] - padfX[i];
2347 104691 : const double dfRY = padfY[ji] - padfY[i];
2348 :
2349 104691 : dfAccumulator += sqrt(dfRX * dfRX + dfRY * dfRY);
2350 104691 : n++;
2351 : }
2352 : }
2353 : }
2354 : }
2355 : }
2356 800 : CPLFree(papsPoints);
2357 : }
2358 : else
2359 : {
2360 62763 : GUInt32 i = 0;
2361 372600 : while (i < nPoints - 1)
2362 : {
2363 309837 : double dfRX1 = padfX[i] - dfXPoint;
2364 309837 : double dfRY1 = padfY[i] - dfYPoint;
2365 :
2366 309837 : if (bRotated)
2367 : {
2368 123998 : const double dfRXRotated = dfRX1 * dfCoeff1 + dfRY1 * dfCoeff2;
2369 123998 : const double dfRYRotated = dfRY1 * dfCoeff1 - dfRX1 * dfCoeff2;
2370 :
2371 123998 : dfRX1 = dfRXRotated;
2372 123998 : dfRY1 = dfRYRotated;
2373 : }
2374 :
2375 : // Is this point located inside the search ellipse?
2376 309837 : if (dfRadius2Square * dfRX1 * dfRX1 +
2377 309837 : 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 750643 : for (GUInt32 j = i + 1; j < nPoints; j++)
2383 : {
2384 549733 : double dfRX2 = padfX[j] - dfXPoint;
2385 549733 : double dfRY2 = padfY[j] - dfYPoint;
2386 :
2387 549733 : if (bRotated)
2388 : {
2389 227352 : const double dfRXRotated =
2390 227352 : dfRX2 * dfCoeff1 + dfRY2 * dfCoeff2;
2391 227352 : const double dfRYRotated =
2392 227352 : dfRY2 * dfCoeff1 - dfRX2 * dfCoeff2;
2393 :
2394 227352 : dfRX2 = dfRXRotated;
2395 227352 : dfRY2 = dfRYRotated;
2396 : }
2397 :
2398 549733 : if (dfRadius2Square * dfRX2 * dfRX2 +
2399 549733 : dfRadius1Square * dfRY2 * dfRY2 <=
2400 : dfR12Square)
2401 : {
2402 327684 : const double dfRX = padfX[j] - padfX[i];
2403 327684 : const double dfRY = padfY[j] - padfY[i];
2404 :
2405 327684 : dfAccumulator += sqrt(dfRX * dfRX + dfRY * dfRY);
2406 327684 : n++;
2407 : }
2408 : }
2409 : }
2410 :
2411 309837 : i++;
2412 : }
2413 : }
2414 :
2415 : // Search for the first point within the search ellipse.
2416 63560 : if (n < poOptions->nMinPoints || n == 0)
2417 : {
2418 2572 : *pdfValue = poOptions->dfNoDataValue;
2419 : }
2420 : else
2421 : {
2422 60988 : *pdfValue = dfAccumulator / n;
2423 : }
2424 :
2425 63560 : 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 311125 : 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 311125 : GDALGridExtraParameters *psExtraParams =
2465 : static_cast<GDALGridExtraParameters *>(hExtraParams);
2466 311125 : GDALTriangulation *psTriangulation = psExtraParams->psTriangulation;
2467 :
2468 311125 : int nOutputFacetIdx = -1;
2469 311125 : const bool bRet = CPL_TO_BOOL(GDALTriangulationFindFacetDirected(
2470 : psTriangulation, psExtraParams->nInitialFacetIdx, dfXPoint, dfYPoint,
2471 : &nOutputFacetIdx));
2472 :
2473 308737 : if (bRet)
2474 : {
2475 74078 : CPLAssert(nOutputFacetIdx >= 0);
2476 : // Reuse output facet idx as next initial index since we proceed line by
2477 : // line.
2478 74078 : psExtraParams->nInitialFacetIdx = nOutputFacetIdx;
2479 :
2480 74078 : double lambda1 = 0.0;
2481 74078 : double lambda2 = 0.0;
2482 74078 : double lambda3 = 0.0;
2483 74078 : GDALTriangulationComputeBarycentricCoordinates(
2484 : psTriangulation, nOutputFacetIdx, dfXPoint, dfYPoint, &lambda1,
2485 : &lambda2, &lambda3);
2486 74525 : const int i1 =
2487 74525 : psTriangulation->pasFacets[nOutputFacetIdx].anVertexIdx[0];
2488 74525 : const int i2 =
2489 74525 : psTriangulation->pasFacets[nOutputFacetIdx].anVertexIdx[1];
2490 74525 : const int i3 =
2491 74525 : psTriangulation->pasFacets[nOutputFacetIdx].anVertexIdx[2];
2492 74525 : *pdfValue =
2493 74525 : lambda1 * padfZ[i1] + lambda2 * padfZ[i2] + lambda3 * padfZ[i3];
2494 : }
2495 : else
2496 : {
2497 234659 : if (nOutputFacetIdx >= 0)
2498 : {
2499 : // Also reuse this failed output facet, when valid, as seed for
2500 : // next search.
2501 226308 : psExtraParams->nInitialFacetIdx = nOutputFacetIdx;
2502 : }
2503 :
2504 234659 : const GDALGridLinearOptions *const poOptions =
2505 : static_cast<const GDALGridLinearOptions *>(poOptionsIn);
2506 234659 : const double dfRadius = poOptions->dfRadius;
2507 234659 : if (dfRadius == 0.0)
2508 : {
2509 113062 : *pdfValue = poOptions->dfNoDataValue;
2510 : }
2511 : else
2512 : {
2513 : GDALGridNearestNeighborOptions sNeighbourOptions;
2514 121597 : sNeighbourOptions.nSizeOfStructure = sizeof(sNeighbourOptions);
2515 115792 : sNeighbourOptions.dfRadius1 =
2516 114870 : dfRadius < 0.0 || dfRadius >= std::numeric_limits<double>::max()
2517 230662 : ? 0.0
2518 : : dfRadius;
2519 115660 : sNeighbourOptions.dfRadius2 =
2520 111991 : dfRadius < 0.0 || dfRadius >= std::numeric_limits<double>::max()
2521 227651 : ? 0.0
2522 : : dfRadius;
2523 115660 : sNeighbourOptions.dfAngle = 0.0;
2524 115660 : sNeighbourOptions.dfNoDataValue = poOptions->dfNoDataValue;
2525 115660 : return GDALGridNearestNeighbor(&sNeighbourOptions, nPoints, padfX,
2526 : padfY, padfZ, dfXPoint, dfYPoint,
2527 128987 : pdfValue, hExtraParams);
2528 : }
2529 : }
2530 :
2531 187587 : return CE_None;
2532 : }
2533 :
2534 : /************************************************************************/
2535 : /* GDALGridJob */
2536 : /************************************************************************/
2537 :
2538 : typedef struct _GDALGridJob GDALGridJob;
2539 :
2540 : struct _GDALGridJob
2541 : {
2542 : GUInt32 nYStart;
2543 :
2544 : GByte *pabyData;
2545 : GUInt32 nYStep;
2546 : GUInt32 nXSize;
2547 : GUInt32 nYSize;
2548 : double dfXMin;
2549 : double dfYMin;
2550 : double dfDeltaX;
2551 : double dfDeltaY;
2552 : GUInt32 nPoints;
2553 : const double *padfX;
2554 : const double *padfY;
2555 : const double *padfZ;
2556 : const void *poOptions;
2557 : GDALGridFunction pfnGDALGridMethod;
2558 : GDALGridExtraParameters *psExtraParameters;
2559 : int (*pfnProgress)(GDALGridJob *psJob);
2560 : GDALDataType eType;
2561 :
2562 : int *pnCounter;
2563 : volatile int *pbStop;
2564 : CPLCond *hCond;
2565 : CPLMutex *hCondMutex;
2566 :
2567 : GDALProgressFunc pfnRealProgress;
2568 : void *pRealProgressArg;
2569 : };
2570 :
2571 : /************************************************************************/
2572 : /* GDALGridProgressMultiThread() */
2573 : /************************************************************************/
2574 :
2575 : // Return TRUE if the computation must be interrupted.
2576 20554 : static int GDALGridProgressMultiThread(GDALGridJob *psJob)
2577 : {
2578 20554 : CPLAcquireMutex(psJob->hCondMutex, 1.0);
2579 20596 : ++(*psJob->pnCounter);
2580 20596 : CPLCondSignal(psJob->hCond);
2581 20596 : const int bStop = *psJob->pbStop;
2582 20596 : CPLReleaseMutex(psJob->hCondMutex);
2583 :
2584 20596 : return bStop;
2585 : }
2586 :
2587 : /************************************************************************/
2588 : /* GDALGridProgressMonoThread() */
2589 : /************************************************************************/
2590 :
2591 : // Return TRUE if the computation must be interrupted.
2592 40 : static int GDALGridProgressMonoThread(GDALGridJob *psJob)
2593 : {
2594 : // coverity[missing_lock]
2595 40 : const int nCounter = ++(*psJob->pnCounter);
2596 40 : if (!psJob->pfnRealProgress(nCounter / static_cast<double>(psJob->nYSize),
2597 : "", psJob->pRealProgressArg))
2598 : {
2599 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
2600 0 : *psJob->pbStop = TRUE;
2601 0 : return TRUE;
2602 : }
2603 40 : return FALSE;
2604 : }
2605 :
2606 : /************************************************************************/
2607 : /* GDALGridJobProcess() */
2608 : /************************************************************************/
2609 :
2610 720 : static void GDALGridJobProcess(void *user_data)
2611 : {
2612 720 : GDALGridJob *const psJob = static_cast<GDALGridJob *>(user_data);
2613 720 : int (*pfnProgress)(GDALGridJob * psJob) = psJob->pfnProgress;
2614 720 : const GUInt32 nXSize = psJob->nXSize;
2615 :
2616 : /* -------------------------------------------------------------------- */
2617 : /* Allocate a buffer of scanline size, fill it with gridded values */
2618 : /* and use GDALCopyWords() to copy values into output data array with */
2619 : /* appropriate data type conversion. */
2620 : /* -------------------------------------------------------------------- */
2621 : double *padfValues =
2622 720 : static_cast<double *>(VSI_MALLOC2_VERBOSE(sizeof(double), nXSize));
2623 721 : if (padfValues == nullptr)
2624 : {
2625 0 : *(psJob->pbStop) = TRUE;
2626 0 : if (pfnProgress != nullptr)
2627 0 : pfnProgress(psJob); // To notify the main thread.
2628 0 : return;
2629 : }
2630 :
2631 721 : const GUInt32 nYStart = psJob->nYStart;
2632 721 : const GUInt32 nYStep = psJob->nYStep;
2633 721 : GByte *pabyData = psJob->pabyData;
2634 :
2635 721 : const GUInt32 nYSize = psJob->nYSize;
2636 721 : const double dfXMin = psJob->dfXMin;
2637 721 : const double dfYMin = psJob->dfYMin;
2638 721 : const double dfDeltaX = psJob->dfDeltaX;
2639 721 : const double dfDeltaY = psJob->dfDeltaY;
2640 721 : const GUInt32 nPoints = psJob->nPoints;
2641 721 : const double *padfX = psJob->padfX;
2642 721 : const double *padfY = psJob->padfY;
2643 721 : const double *padfZ = psJob->padfZ;
2644 721 : const void *poOptions = psJob->poOptions;
2645 721 : GDALGridFunction pfnGDALGridMethod = psJob->pfnGDALGridMethod;
2646 : // Have a local copy of sExtraParameters since we want to modify
2647 : // nInitialFacetIdx.
2648 721 : GDALGridExtraParameters sExtraParameters = *psJob->psExtraParameters;
2649 721 : const GDALDataType eType = psJob->eType;
2650 :
2651 721 : const int nDataTypeSize = GDALGetDataTypeSizeBytes(eType);
2652 717 : const int nLineSpace = nXSize * nDataTypeSize;
2653 :
2654 21344 : for (GUInt32 nYPoint = nYStart; nYPoint < nYSize; nYPoint += nYStep)
2655 : {
2656 20630 : const double dfYPoint = dfYMin + (nYPoint + 0.5) * dfDeltaY;
2657 :
2658 4429540 : for (GUInt32 nXPoint = 0; nXPoint < nXSize; nXPoint++)
2659 : {
2660 4446500 : const double dfXPoint = dfXMin + (nXPoint + 0.5) * dfDeltaX;
2661 :
2662 8855410 : if ((*pfnGDALGridMethod)(poOptions, nPoints, padfX, padfY, padfZ,
2663 4446500 : dfXPoint, dfYPoint, padfValues + nXPoint,
2664 4408900 : &sExtraParameters) != CE_None)
2665 : {
2666 0 : CPLError(CE_Failure, CPLE_AppDefined,
2667 : "Gridding failed at X position %lu, Y position %lu",
2668 : static_cast<long unsigned int>(nXPoint),
2669 : static_cast<long unsigned int>(nYPoint));
2670 0 : *psJob->pbStop = TRUE;
2671 0 : if (pfnProgress != nullptr)
2672 0 : pfnProgress(psJob); // To notify the main thread.
2673 0 : break;
2674 : }
2675 : }
2676 :
2677 0 : GDALCopyWords(padfValues, GDT_Float64, sizeof(double),
2678 0 : pabyData + nYPoint * nLineSpace, eType, nDataTypeSize,
2679 : nXSize);
2680 :
2681 20584 : if (*psJob->pbStop || (pfnProgress != nullptr && pfnProgress(psJob)))
2682 0 : break;
2683 : }
2684 :
2685 714 : CPLFree(padfValues);
2686 : }
2687 :
2688 : /************************************************************************/
2689 : /* GDALGridContextCreate() */
2690 : /************************************************************************/
2691 :
2692 : struct GDALGridContext
2693 : {
2694 : GDALGridAlgorithm eAlgorithm;
2695 : void *poOptions;
2696 : GDALGridFunction pfnGDALGridMethod;
2697 :
2698 : GUInt32 nPoints;
2699 : GDALGridPoint *pasGridPoints;
2700 : GDALGridXYArrays sXYArrays;
2701 :
2702 : GDALGridExtraParameters sExtraParameters;
2703 : double *padfX;
2704 : double *padfY;
2705 : double *padfZ;
2706 : bool bFreePadfXYZArrays;
2707 :
2708 : CPLWorkerThreadPool *poWorkerThreadPool;
2709 : };
2710 :
2711 : static void GDALGridContextCreateQuadTree(GDALGridContext *psContext);
2712 :
2713 : /**
2714 : * Creates a context to do regular gridding from the scattered data.
2715 : *
2716 : * This function takes the arrays of X and Y coordinates and corresponding Z
2717 : * values as input to prepare computation of regular grid (or call it a raster)
2718 : * from these scattered data.
2719 : *
2720 : * On Intel/AMD i386/x86_64 architectures, some
2721 : * gridding methods will be optimized with SSE instructions (provided GDAL
2722 : * has been compiled with such support, and it is available at runtime).
2723 : * Currently, only 'invdist' algorithm with default parameters has an optimized
2724 : * implementation.
2725 : * This can provide substantial speed-up, but sometimes at the expense of
2726 : * reduced floating point precision. This can be disabled by setting the
2727 : * GDAL_USE_SSE configuration option to NO.
2728 : * A further optimized version can use the AVX
2729 : * instruction set. This can be disabled by setting the GDAL_USE_AVX
2730 : * configuration option to NO.
2731 : *
2732 : * It is possible to set the GDAL_NUM_THREADS
2733 : * configuration option to parallelize the processing. The value to set is
2734 : * the number of worker threads, or ALL_CPUS to use all the cores/CPUs of the
2735 : * computer (default value).
2736 : *
2737 : * @param eAlgorithm Gridding method.
2738 : * @param poOptions Options to control chosen gridding method.
2739 : * @param nPoints Number of elements in input arrays.
2740 : * @param padfX Input array of X coordinates.
2741 : * @param padfY Input array of Y coordinates.
2742 : * @param padfZ Input array of Z values.
2743 : * @param bCallerWillKeepPointArraysAlive Whether the provided padfX, padfY,
2744 : * padfZ arrays will still be "alive" during the calls to
2745 : * GDALGridContextProcess(). Setting to TRUE prevent them from being
2746 : * duplicated in the context. If unsure, set to FALSE.
2747 : *
2748 : * @return the context (to be freed with GDALGridContextFree()) or NULL in case
2749 : * or error.
2750 : *
2751 : * @since GDAL 2.1
2752 : */
2753 :
2754 182 : GDALGridContext *GDALGridContextCreate(GDALGridAlgorithm eAlgorithm,
2755 : const void *poOptions, GUInt32 nPoints,
2756 : const double *padfX, const double *padfY,
2757 : const double *padfZ,
2758 : int bCallerWillKeepPointArraysAlive)
2759 : {
2760 182 : CPLAssert(poOptions);
2761 182 : CPLAssert(padfX);
2762 182 : CPLAssert(padfY);
2763 182 : CPLAssert(padfZ);
2764 182 : bool bCreateQuadTree = false;
2765 :
2766 : const unsigned int nPointCountThreshold =
2767 182 : atoi(CPLGetConfigOption("GDAL_GRID_POINT_COUNT_THRESHOLD", "100"));
2768 :
2769 : // Starting address aligned on 32-byte boundary for AVX.
2770 182 : float *pafXAligned = nullptr;
2771 182 : float *pafYAligned = nullptr;
2772 182 : float *pafZAligned = nullptr;
2773 :
2774 182 : void *poOptionsNew = nullptr;
2775 :
2776 182 : GDALGridFunction pfnGDALGridMethod = nullptr;
2777 :
2778 182 : switch (eAlgorithm)
2779 : {
2780 44 : case GGA_InverseDistanceToAPower:
2781 : {
2782 44 : const auto poOptionsOld =
2783 : static_cast<const GDALGridInverseDistanceToAPowerOptions *>(
2784 : poOptions);
2785 44 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
2786 : {
2787 0 : CPLError(CE_Failure, CPLE_AppDefined,
2788 : "Wrong value of nSizeOfStructure member");
2789 0 : return nullptr;
2790 : }
2791 : poOptionsNew =
2792 44 : CPLMalloc(sizeof(GDALGridInverseDistanceToAPowerOptions));
2793 44 : memcpy(poOptionsNew, poOptions,
2794 : sizeof(GDALGridInverseDistanceToAPowerOptions));
2795 :
2796 44 : const GDALGridInverseDistanceToAPowerOptions *const poPower =
2797 : static_cast<const GDALGridInverseDistanceToAPowerOptions *>(
2798 : poOptions);
2799 44 : if (poPower->dfRadius1 == 0.0 && poPower->dfRadius2 == 0.0)
2800 : {
2801 35 : const double dfPower = poPower->dfPower;
2802 35 : const double dfSmoothing = poPower->dfSmoothing;
2803 :
2804 35 : pfnGDALGridMethod = GDALGridInverseDistanceToAPowerNoSearch;
2805 35 : if (dfPower == 2.0 && dfSmoothing == 0.0)
2806 : {
2807 : #ifdef HAVE_AVX_AT_COMPILE_TIME
2808 :
2809 33 : if (CPLTestBool(
2810 60 : CPLGetConfigOption("GDAL_USE_AVX", "YES")) &&
2811 27 : CPLHaveRuntimeAVX())
2812 : {
2813 : pafXAligned = static_cast<float *>(
2814 27 : VSI_MALLOC_ALIGNED_AUTO_VERBOSE(sizeof(float) *
2815 : nPoints));
2816 : pafYAligned = static_cast<float *>(
2817 27 : VSI_MALLOC_ALIGNED_AUTO_VERBOSE(sizeof(float) *
2818 : nPoints));
2819 : pafZAligned = static_cast<float *>(
2820 27 : VSI_MALLOC_ALIGNED_AUTO_VERBOSE(sizeof(float) *
2821 : nPoints));
2822 27 : if (pafXAligned != nullptr && pafYAligned != nullptr &&
2823 : pafZAligned != nullptr)
2824 : {
2825 27 : CPLDebug("GDAL_GRID",
2826 : "Using AVX optimized version");
2827 27 : pfnGDALGridMethod =
2828 : GDALGridInverseDistanceToAPower2NoSmoothingNoSearchAVX;
2829 1335 : for (GUInt32 i = 0; i < nPoints; i++)
2830 : {
2831 1308 : pafXAligned[i] = static_cast<float>(padfX[i]);
2832 1308 : pafYAligned[i] = static_cast<float>(padfY[i]);
2833 1308 : pafZAligned[i] = static_cast<float>(padfZ[i]);
2834 27 : }
2835 : }
2836 : else
2837 : {
2838 0 : VSIFree(pafXAligned);
2839 0 : VSIFree(pafYAligned);
2840 0 : VSIFree(pafZAligned);
2841 0 : pafXAligned = nullptr;
2842 0 : pafYAligned = nullptr;
2843 0 : pafZAligned = nullptr;
2844 : }
2845 : }
2846 : #endif
2847 :
2848 : #ifdef HAVE_SSE_AT_COMPILE_TIME
2849 :
2850 39 : if (pafXAligned == nullptr &&
2851 6 : CPLTestBool(CPLGetConfigOption("GDAL_USE_SSE", "YES"))
2852 : #if !defined(USE_NEON_OPTIMIZATIONS)
2853 39 : && CPLHaveRuntimeSSE()
2854 : #endif
2855 : )
2856 : {
2857 : pafXAligned = static_cast<float *>(
2858 3 : VSI_MALLOC_ALIGNED_AUTO_VERBOSE(sizeof(float) *
2859 : nPoints));
2860 : pafYAligned = static_cast<float *>(
2861 3 : VSI_MALLOC_ALIGNED_AUTO_VERBOSE(sizeof(float) *
2862 : nPoints));
2863 : pafZAligned = static_cast<float *>(
2864 3 : VSI_MALLOC_ALIGNED_AUTO_VERBOSE(sizeof(float) *
2865 : nPoints));
2866 3 : if (pafXAligned != nullptr && pafYAligned != nullptr &&
2867 : pafZAligned != nullptr)
2868 : {
2869 3 : CPLDebug("GDAL_GRID",
2870 : "Using SSE optimized version");
2871 3 : pfnGDALGridMethod =
2872 : GDALGridInverseDistanceToAPower2NoSmoothingNoSearchSSE;
2873 405 : for (GUInt32 i = 0; i < nPoints; i++)
2874 : {
2875 402 : pafXAligned[i] = static_cast<float>(padfX[i]);
2876 402 : pafYAligned[i] = static_cast<float>(padfY[i]);
2877 402 : pafZAligned[i] = static_cast<float>(padfZ[i]);
2878 3 : }
2879 : }
2880 : else
2881 : {
2882 0 : VSIFree(pafXAligned);
2883 0 : VSIFree(pafYAligned);
2884 0 : VSIFree(pafZAligned);
2885 0 : pafXAligned = nullptr;
2886 0 : pafYAligned = nullptr;
2887 0 : pafZAligned = nullptr;
2888 : }
2889 : }
2890 : #endif // HAVE_SSE_AT_COMPILE_TIME
2891 35 : }
2892 : }
2893 : else
2894 : {
2895 9 : pfnGDALGridMethod = GDALGridInverseDistanceToAPower;
2896 : }
2897 44 : break;
2898 : }
2899 22 : case GGA_InverseDistanceToAPowerNearestNeighbor:
2900 : {
2901 22 : const auto poOptionsOld = static_cast<
2902 : const GDALGridInverseDistanceToAPowerNearestNeighborOptions *>(
2903 : poOptions);
2904 22 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
2905 : {
2906 0 : CPLError(CE_Failure, CPLE_AppDefined,
2907 : "Wrong value of nSizeOfStructure member");
2908 0 : return nullptr;
2909 : }
2910 22 : poOptionsNew = CPLMalloc(
2911 : sizeof(GDALGridInverseDistanceToAPowerNearestNeighborOptions));
2912 22 : memcpy(
2913 : poOptionsNew, poOptions,
2914 : sizeof(GDALGridInverseDistanceToAPowerNearestNeighborOptions));
2915 :
2916 22 : if (poOptionsOld->nMinPointsPerQuadrant != 0 ||
2917 12 : poOptionsOld->nMaxPointsPerQuadrant != 0)
2918 : {
2919 12 : pfnGDALGridMethod =
2920 : GDALGridInverseDistanceToAPowerNearestNeighborPerQuadrant;
2921 : }
2922 : else
2923 : {
2924 10 : pfnGDALGridMethod =
2925 : GDALGridInverseDistanceToAPowerNearestNeighbor;
2926 : }
2927 22 : bCreateQuadTree = true;
2928 22 : break;
2929 : }
2930 25 : case GGA_MovingAverage:
2931 : {
2932 25 : const auto poOptionsOld =
2933 : static_cast<const GDALGridMovingAverageOptions *>(poOptions);
2934 25 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
2935 : {
2936 0 : CPLError(CE_Failure, CPLE_AppDefined,
2937 : "Wrong value of nSizeOfStructure member");
2938 0 : return nullptr;
2939 : }
2940 25 : poOptionsNew = CPLMalloc(sizeof(GDALGridMovingAverageOptions));
2941 25 : memcpy(poOptionsNew, poOptions,
2942 : sizeof(GDALGridMovingAverageOptions));
2943 :
2944 25 : if (poOptionsOld->nMinPointsPerQuadrant != 0 ||
2945 17 : poOptionsOld->nMaxPointsPerQuadrant != 0)
2946 : {
2947 9 : pfnGDALGridMethod = GDALGridMovingAveragePerQuadrant;
2948 9 : bCreateQuadTree = true;
2949 : }
2950 : else
2951 : {
2952 16 : pfnGDALGridMethod = GDALGridMovingAverage;
2953 16 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
2954 21 : poOptionsOld->dfAngle == 0.0 &&
2955 5 : (poOptionsOld->dfRadius1 > 0.0 ||
2956 1 : poOptionsOld->dfRadius2 > 0.0));
2957 : }
2958 25 : break;
2959 : }
2960 21 : case GGA_NearestNeighbor:
2961 : {
2962 21 : const auto poOptionsOld =
2963 : static_cast<const GDALGridNearestNeighborOptions *>(poOptions);
2964 21 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
2965 : {
2966 0 : CPLError(CE_Failure, CPLE_AppDefined,
2967 : "Wrong value of nSizeOfStructure member");
2968 0 : return nullptr;
2969 : }
2970 21 : poOptionsNew = CPLMalloc(sizeof(GDALGridNearestNeighborOptions));
2971 21 : memcpy(poOptionsNew, poOptions,
2972 : sizeof(GDALGridNearestNeighborOptions));
2973 :
2974 21 : pfnGDALGridMethod = GDALGridNearestNeighbor;
2975 34 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
2976 34 : poOptionsOld->dfAngle == 0.0 &&
2977 13 : (poOptionsOld->dfRadius1 > 0.0 ||
2978 5 : poOptionsOld->dfRadius2 > 0.0));
2979 21 : break;
2980 : }
2981 18 : case GGA_MetricMinimum:
2982 : {
2983 18 : const auto poOptionsOld =
2984 : static_cast<const GDALGridDataMetricsOptions *>(poOptions);
2985 18 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
2986 : {
2987 0 : CPLError(CE_Failure, CPLE_AppDefined,
2988 : "Wrong value of nSizeOfStructure member");
2989 0 : return nullptr;
2990 : }
2991 18 : poOptionsNew = CPLMalloc(sizeof(GDALGridDataMetricsOptions));
2992 18 : memcpy(poOptionsNew, poOptions, sizeof(GDALGridDataMetricsOptions));
2993 :
2994 18 : if (poOptionsOld->nMinPointsPerQuadrant != 0 ||
2995 12 : poOptionsOld->nMaxPointsPerQuadrant != 0)
2996 : {
2997 7 : pfnGDALGridMethod = GDALGridDataMetricMinimumPerQuadrant;
2998 7 : bCreateQuadTree = true;
2999 : }
3000 : else
3001 : {
3002 11 : pfnGDALGridMethod = GDALGridDataMetricMinimum;
3003 11 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
3004 16 : poOptionsOld->dfAngle == 0.0 &&
3005 5 : (poOptionsOld->dfRadius1 > 0.0 ||
3006 1 : poOptionsOld->dfRadius2 > 0.0));
3007 : }
3008 18 : break;
3009 : }
3010 12 : case GGA_MetricMaximum:
3011 : {
3012 12 : const auto poOptionsOld =
3013 : static_cast<const GDALGridDataMetricsOptions *>(poOptions);
3014 12 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
3015 : {
3016 0 : CPLError(CE_Failure, CPLE_AppDefined,
3017 : "Wrong value of nSizeOfStructure member");
3018 0 : return nullptr;
3019 : }
3020 12 : poOptionsNew = CPLMalloc(sizeof(GDALGridDataMetricsOptions));
3021 12 : memcpy(poOptionsNew, poOptions, sizeof(GDALGridDataMetricsOptions));
3022 :
3023 12 : if (poOptionsOld->nMinPointsPerQuadrant != 0 ||
3024 7 : poOptionsOld->nMaxPointsPerQuadrant != 0)
3025 : {
3026 5 : pfnGDALGridMethod = GDALGridDataMetricMaximumPerQuadrant;
3027 5 : bCreateQuadTree = true;
3028 : }
3029 : else
3030 : {
3031 7 : pfnGDALGridMethod = GDALGridDataMetricMaximum;
3032 7 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
3033 11 : poOptionsOld->dfAngle == 0.0 &&
3034 4 : (poOptionsOld->dfRadius1 > 0.0 ||
3035 1 : poOptionsOld->dfRadius2 > 0.0));
3036 : }
3037 :
3038 12 : break;
3039 : }
3040 9 : case GGA_MetricRange:
3041 : {
3042 9 : const auto poOptionsOld =
3043 : static_cast<const GDALGridDataMetricsOptions *>(poOptions);
3044 9 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
3045 : {
3046 0 : CPLError(CE_Failure, CPLE_AppDefined,
3047 : "Wrong value of nSizeOfStructure member");
3048 0 : return nullptr;
3049 : }
3050 9 : poOptionsNew = CPLMalloc(sizeof(GDALGridDataMetricsOptions));
3051 9 : memcpy(poOptionsNew, poOptions, sizeof(GDALGridDataMetricsOptions));
3052 :
3053 9 : if (poOptionsOld->nMinPointsPerQuadrant != 0 ||
3054 4 : poOptionsOld->nMaxPointsPerQuadrant != 0)
3055 : {
3056 5 : pfnGDALGridMethod = GDALGridDataMetricRangePerQuadrant;
3057 5 : bCreateQuadTree = true;
3058 : }
3059 : else
3060 : {
3061 4 : pfnGDALGridMethod = GDALGridDataMetricRange;
3062 4 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
3063 7 : poOptionsOld->dfAngle == 0.0 &&
3064 3 : (poOptionsOld->dfRadius1 > 0.0 ||
3065 1 : poOptionsOld->dfRadius2 > 0.0));
3066 : }
3067 :
3068 9 : break;
3069 : }
3070 11 : case GGA_MetricCount:
3071 : {
3072 11 : const auto poOptionsOld =
3073 : static_cast<const GDALGridDataMetricsOptions *>(poOptions);
3074 11 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
3075 : {
3076 0 : CPLError(CE_Failure, CPLE_AppDefined,
3077 : "Wrong value of nSizeOfStructure member");
3078 0 : return nullptr;
3079 : }
3080 11 : poOptionsNew = CPLMalloc(sizeof(GDALGridDataMetricsOptions));
3081 11 : memcpy(poOptionsNew, poOptions, sizeof(GDALGridDataMetricsOptions));
3082 :
3083 11 : if (poOptionsOld->nMinPointsPerQuadrant != 0 ||
3084 6 : poOptionsOld->nMaxPointsPerQuadrant != 0)
3085 : {
3086 5 : pfnGDALGridMethod = GDALGridDataMetricCountPerQuadrant;
3087 5 : bCreateQuadTree = true;
3088 : }
3089 : else
3090 : {
3091 6 : pfnGDALGridMethod = GDALGridDataMetricCount;
3092 6 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
3093 11 : poOptionsOld->dfAngle == 0.0 &&
3094 5 : (poOptionsOld->dfRadius1 > 0.0 ||
3095 0 : poOptionsOld->dfRadius2 > 0.0));
3096 : }
3097 :
3098 11 : break;
3099 : }
3100 9 : case GGA_MetricAverageDistance:
3101 : {
3102 9 : const auto poOptionsOld =
3103 : static_cast<const GDALGridDataMetricsOptions *>(poOptions);
3104 9 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
3105 : {
3106 0 : CPLError(CE_Failure, CPLE_AppDefined,
3107 : "Wrong value of nSizeOfStructure member");
3108 0 : return nullptr;
3109 : }
3110 9 : poOptionsNew = CPLMalloc(sizeof(GDALGridDataMetricsOptions));
3111 9 : memcpy(poOptionsNew, poOptions, sizeof(GDALGridDataMetricsOptions));
3112 :
3113 9 : if (poOptionsOld->nMinPointsPerQuadrant != 0 ||
3114 4 : poOptionsOld->nMaxPointsPerQuadrant != 0)
3115 : {
3116 5 : pfnGDALGridMethod =
3117 : GDALGridDataMetricAverageDistancePerQuadrant;
3118 5 : bCreateQuadTree = true;
3119 : }
3120 : else
3121 : {
3122 4 : pfnGDALGridMethod = GDALGridDataMetricAverageDistance;
3123 4 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
3124 7 : poOptionsOld->dfAngle == 0.0 &&
3125 3 : (poOptionsOld->dfRadius1 > 0.0 ||
3126 1 : poOptionsOld->dfRadius2 > 0.0));
3127 : }
3128 :
3129 9 : break;
3130 : }
3131 4 : case GGA_MetricAverageDistancePts:
3132 : {
3133 4 : const auto poOptionsOld =
3134 : static_cast<const GDALGridDataMetricsOptions *>(poOptions);
3135 4 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
3136 : {
3137 0 : CPLError(CE_Failure, CPLE_AppDefined,
3138 : "Wrong value of nSizeOfStructure member");
3139 0 : return nullptr;
3140 : }
3141 4 : poOptionsNew = CPLMalloc(sizeof(GDALGridDataMetricsOptions));
3142 4 : memcpy(poOptionsNew, poOptions, sizeof(GDALGridDataMetricsOptions));
3143 :
3144 4 : pfnGDALGridMethod = GDALGridDataMetricAverageDistancePts;
3145 7 : bCreateQuadTree = (nPoints > nPointCountThreshold &&
3146 6 : poOptionsOld->dfAngle == 0.0 &&
3147 2 : (poOptionsOld->dfRadius1 > 0.0 ||
3148 0 : poOptionsOld->dfRadius2 > 0.0));
3149 :
3150 4 : break;
3151 : }
3152 7 : case GGA_Linear:
3153 : {
3154 7 : const auto poOptionsOld =
3155 : static_cast<const GDALGridLinearOptions *>(poOptions);
3156 7 : if (poOptionsOld->nSizeOfStructure != sizeof(*poOptionsOld))
3157 : {
3158 0 : CPLError(CE_Failure, CPLE_AppDefined,
3159 : "Wrong value of nSizeOfStructure member");
3160 0 : return nullptr;
3161 : }
3162 7 : poOptionsNew = CPLMalloc(sizeof(GDALGridLinearOptions));
3163 7 : memcpy(poOptionsNew, poOptions, sizeof(GDALGridLinearOptions));
3164 :
3165 7 : pfnGDALGridMethod = GDALGridLinear;
3166 7 : break;
3167 : }
3168 0 : default:
3169 0 : CPLError(CE_Failure, CPLE_IllegalArg,
3170 : "GDAL does not support gridding method %d", eAlgorithm);
3171 0 : return nullptr;
3172 : }
3173 :
3174 182 : if (pafXAligned == nullptr && !bCallerWillKeepPointArraysAlive)
3175 : {
3176 : double *padfXNew =
3177 0 : static_cast<double *>(VSI_MALLOC2_VERBOSE(nPoints, sizeof(double)));
3178 : double *padfYNew =
3179 0 : static_cast<double *>(VSI_MALLOC2_VERBOSE(nPoints, sizeof(double)));
3180 : double *padfZNew =
3181 0 : static_cast<double *>(VSI_MALLOC2_VERBOSE(nPoints, sizeof(double)));
3182 0 : if (padfXNew == nullptr || padfYNew == nullptr || padfZNew == nullptr)
3183 : {
3184 0 : VSIFree(padfXNew);
3185 0 : VSIFree(padfYNew);
3186 0 : VSIFree(padfZNew);
3187 0 : CPLFree(poOptionsNew);
3188 0 : return nullptr;
3189 : }
3190 0 : memcpy(padfXNew, padfX, nPoints * sizeof(double));
3191 0 : memcpy(padfYNew, padfY, nPoints * sizeof(double));
3192 0 : memcpy(padfZNew, padfZ, nPoints * sizeof(double));
3193 0 : padfX = padfXNew;
3194 0 : padfY = padfYNew;
3195 0 : padfZ = padfZNew;
3196 : }
3197 : GDALGridContext *psContext =
3198 182 : static_cast<GDALGridContext *>(CPLCalloc(1, sizeof(GDALGridContext)));
3199 182 : psContext->eAlgorithm = eAlgorithm;
3200 182 : psContext->poOptions = poOptionsNew;
3201 182 : psContext->pfnGDALGridMethod = pfnGDALGridMethod;
3202 182 : psContext->nPoints = nPoints;
3203 182 : psContext->pasGridPoints = nullptr;
3204 182 : psContext->sXYArrays.padfX = padfX;
3205 182 : psContext->sXYArrays.padfY = padfY;
3206 182 : psContext->sExtraParameters.hQuadTree = nullptr;
3207 182 : psContext->sExtraParameters.dfInitialSearchRadius = 0.0;
3208 182 : psContext->sExtraParameters.pafX = pafXAligned;
3209 182 : psContext->sExtraParameters.pafY = pafYAligned;
3210 182 : psContext->sExtraParameters.pafZ = pafZAligned;
3211 182 : psContext->sExtraParameters.psTriangulation = nullptr;
3212 182 : psContext->sExtraParameters.nInitialFacetIdx = 0;
3213 182 : psContext->padfX = pafXAligned ? nullptr : const_cast<double *>(padfX);
3214 182 : psContext->padfY = pafXAligned ? nullptr : const_cast<double *>(padfY);
3215 182 : psContext->padfZ = pafXAligned ? nullptr : const_cast<double *>(padfZ);
3216 182 : psContext->bFreePadfXYZArrays =
3217 182 : pafXAligned ? false : !bCallerWillKeepPointArraysAlive;
3218 :
3219 : /* -------------------------------------------------------------------- */
3220 : /* Create quadtree if requested and possible. */
3221 : /* -------------------------------------------------------------------- */
3222 182 : if (bCreateQuadTree)
3223 : {
3224 88 : GDALGridContextCreateQuadTree(psContext);
3225 88 : if (psContext->sExtraParameters.hQuadTree == nullptr &&
3226 0 : (eAlgorithm == GGA_InverseDistanceToAPowerNearestNeighbor ||
3227 : pfnGDALGridMethod == GDALGridMovingAveragePerQuadrant))
3228 : {
3229 : // shouldn't happen unless memory allocation failure occurs
3230 0 : GDALGridContextFree(psContext);
3231 0 : return nullptr;
3232 : }
3233 : }
3234 :
3235 : /* -------------------------------------------------------------------- */
3236 : /* Pre-compute extra parameters in GDALGridExtraParameters */
3237 : /* -------------------------------------------------------------------- */
3238 182 : if (eAlgorithm == GGA_InverseDistanceToAPowerNearestNeighbor)
3239 : {
3240 22 : const double dfPower =
3241 : static_cast<
3242 : const GDALGridInverseDistanceToAPowerNearestNeighborOptions *>(
3243 : poOptions)
3244 : ->dfPower;
3245 22 : psContext->sExtraParameters.dfPowerDiv2PreComp = dfPower / 2;
3246 :
3247 22 : const double dfRadius =
3248 : static_cast<
3249 : const GDALGridInverseDistanceToAPowerNearestNeighborOptions *>(
3250 : poOptions)
3251 : ->dfRadius;
3252 22 : psContext->sExtraParameters.dfRadiusPower2PreComp = pow(dfRadius, 2);
3253 : }
3254 :
3255 182 : if (eAlgorithm == GGA_Linear)
3256 : {
3257 7 : psContext->sExtraParameters.psTriangulation =
3258 7 : GDALTriangulationCreateDelaunay(nPoints, padfX, padfY);
3259 7 : if (psContext->sExtraParameters.psTriangulation == nullptr)
3260 : {
3261 0 : GDALGridContextFree(psContext);
3262 0 : return nullptr;
3263 : }
3264 7 : GDALTriangulationComputeBarycentricCoefficients(
3265 : psContext->sExtraParameters.psTriangulation, padfX, padfY);
3266 : }
3267 :
3268 : /* -------------------------------------------------------------------- */
3269 : /* Start thread pool. */
3270 : /* -------------------------------------------------------------------- */
3271 182 : const char *pszThreads = CPLGetConfigOption("GDAL_NUM_THREADS", "ALL_CPUS");
3272 182 : int nThreads = 0;
3273 182 : if (EQUAL(pszThreads, "ALL_CPUS"))
3274 180 : nThreads = CPLGetNumCPUs();
3275 : else
3276 2 : nThreads = atoi(pszThreads);
3277 182 : if (nThreads > 128)
3278 0 : nThreads = 128;
3279 182 : if (nThreads > 1)
3280 : {
3281 180 : psContext->poWorkerThreadPool = new CPLWorkerThreadPool();
3282 180 : if (!psContext->poWorkerThreadPool->Setup(nThreads, nullptr, nullptr))
3283 : {
3284 0 : delete psContext->poWorkerThreadPool;
3285 0 : psContext->poWorkerThreadPool = nullptr;
3286 : }
3287 : else
3288 : {
3289 180 : CPLDebug("GDAL_GRID", "Using %d threads", nThreads);
3290 : }
3291 : }
3292 : else
3293 2 : psContext->poWorkerThreadPool = nullptr;
3294 :
3295 182 : return psContext;
3296 : }
3297 :
3298 : /************************************************************************/
3299 : /* GDALGridContextCreateQuadTree() */
3300 : /************************************************************************/
3301 :
3302 94 : void GDALGridContextCreateQuadTree(GDALGridContext *psContext)
3303 : {
3304 94 : const GUInt32 nPoints = psContext->nPoints;
3305 94 : psContext->pasGridPoints = static_cast<GDALGridPoint *>(
3306 94 : VSI_MALLOC2_VERBOSE(nPoints, sizeof(GDALGridPoint)));
3307 94 : if (psContext->pasGridPoints != nullptr)
3308 : {
3309 94 : const double *const padfX = psContext->padfX;
3310 94 : const double *const padfY = psContext->padfY;
3311 :
3312 : // Determine point extents.
3313 : CPLRectObj sRect;
3314 94 : sRect.minx = padfX[0];
3315 94 : sRect.miny = padfY[0];
3316 94 : sRect.maxx = padfX[0];
3317 94 : sRect.maxy = padfY[0];
3318 28131 : for (GUInt32 i = 1; i < nPoints; i++)
3319 : {
3320 28037 : if (padfX[i] < sRect.minx)
3321 40 : sRect.minx = padfX[i];
3322 28037 : if (padfY[i] < sRect.miny)
3323 788 : sRect.miny = padfY[i];
3324 28037 : if (padfX[i] > sRect.maxx)
3325 794 : sRect.maxx = padfX[i];
3326 28037 : if (padfY[i] > sRect.maxy)
3327 30 : sRect.maxy = padfY[i];
3328 : }
3329 :
3330 : // Initial value for search radius is the typical dimension of a
3331 : // "pixel" of the point array (assuming rather uniform distribution).
3332 94 : psContext->sExtraParameters.dfInitialSearchRadius = sqrt(
3333 94 : (sRect.maxx - sRect.minx) * (sRect.maxy - sRect.miny) / nPoints);
3334 :
3335 94 : psContext->sExtraParameters.hQuadTree =
3336 94 : CPLQuadTreeCreate(&sRect, GDALGridGetPointBounds);
3337 :
3338 28225 : for (GUInt32 i = 0; i < nPoints; i++)
3339 : {
3340 28131 : psContext->pasGridPoints[i].psXYArrays = &(psContext->sXYArrays);
3341 28131 : psContext->pasGridPoints[i].i = i;
3342 28131 : CPLQuadTreeInsert(psContext->sExtraParameters.hQuadTree,
3343 28131 : psContext->pasGridPoints + i);
3344 : }
3345 : }
3346 94 : }
3347 :
3348 : /************************************************************************/
3349 : /* GDALGridContextFree() */
3350 : /************************************************************************/
3351 :
3352 : /**
3353 : * Free a context used created by GDALGridContextCreate()
3354 : *
3355 : * @param psContext the context.
3356 : *
3357 : * @since GDAL 2.1
3358 : */
3359 182 : void GDALGridContextFree(GDALGridContext *psContext)
3360 : {
3361 182 : if (psContext)
3362 : {
3363 182 : CPLFree(psContext->poOptions);
3364 182 : CPLFree(psContext->pasGridPoints);
3365 182 : if (psContext->sExtraParameters.hQuadTree != nullptr)
3366 94 : CPLQuadTreeDestroy(psContext->sExtraParameters.hQuadTree);
3367 182 : if (psContext->bFreePadfXYZArrays)
3368 : {
3369 0 : CPLFree(psContext->padfX);
3370 0 : CPLFree(psContext->padfY);
3371 0 : CPLFree(psContext->padfZ);
3372 : }
3373 182 : VSIFreeAligned(psContext->sExtraParameters.pafX);
3374 182 : VSIFreeAligned(psContext->sExtraParameters.pafY);
3375 182 : VSIFreeAligned(psContext->sExtraParameters.pafZ);
3376 182 : if (psContext->sExtraParameters.psTriangulation)
3377 7 : GDALTriangulationFree(psContext->sExtraParameters.psTriangulation);
3378 182 : delete psContext->poWorkerThreadPool;
3379 182 : CPLFree(psContext);
3380 : }
3381 182 : }
3382 :
3383 : /************************************************************************/
3384 : /* GDALGridContextProcess() */
3385 : /************************************************************************/
3386 :
3387 : /**
3388 : * Do the gridding of a window of a raster.
3389 : *
3390 : * This function takes the gridding context as input to preprare computation
3391 : * of regular grid (or call it a raster) from these scattered data.
3392 : * You should supply the extent of the output grid and allocate array
3393 : * sufficient to hold such a grid.
3394 : *
3395 : * @param psContext Gridding context.
3396 : * @param dfXMin Lowest X border of output grid.
3397 : * @param dfXMax Highest X border of output grid.
3398 : * @param dfYMin Lowest Y border of output grid.
3399 : * @param dfYMax Highest Y border of output grid.
3400 : * @param nXSize Number of columns in output grid.
3401 : * @param nYSize Number of rows in output grid.
3402 : * @param eType Data type of output array.
3403 : * @param pData Pointer to array where the computed grid will be stored.
3404 : * @param pfnProgress a GDALProgressFunc() compatible callback function for
3405 : * reporting progress or NULL.
3406 : * @param pProgressArg argument to be passed to pfnProgress. May be NULL.
3407 : *
3408 : * @return CE_None on success or CE_Failure if something goes wrong.
3409 : *
3410 : * @since GDAL 2.1
3411 : */
3412 :
3413 182 : CPLErr GDALGridContextProcess(GDALGridContext *psContext, double dfXMin,
3414 : double dfXMax, double dfYMin, double dfYMax,
3415 : GUInt32 nXSize, GUInt32 nYSize,
3416 : GDALDataType eType, void *pData,
3417 : GDALProgressFunc pfnProgress, void *pProgressArg)
3418 : {
3419 182 : CPLAssert(psContext);
3420 182 : CPLAssert(pData);
3421 :
3422 182 : if (nXSize == 0 || nYSize == 0)
3423 : {
3424 0 : CPLError(CE_Failure, CPLE_IllegalArg,
3425 : "Output raster dimensions should have non-zero size.");
3426 0 : return CE_Failure;
3427 : }
3428 :
3429 182 : const double dfDeltaX = (dfXMax - dfXMin) / nXSize;
3430 182 : const double dfDeltaY = (dfYMax - dfYMin) / nYSize;
3431 :
3432 : // For linear, check if we will need to fallback to nearest neighbour
3433 : // by sampling along the edges. If all points on edges are within
3434 : // triangles, then interior points will also be.
3435 182 : if (psContext->eAlgorithm == GGA_Linear &&
3436 7 : psContext->sExtraParameters.hQuadTree == nullptr)
3437 : {
3438 7 : bool bNeedNearest = false;
3439 7 : int nStartLeft = 0;
3440 7 : int nStartRight = 0;
3441 7 : const double dfXPointMin = dfXMin + (0 + 0.5) * dfDeltaX;
3442 7 : const double dfXPointMax = dfXMin + (nXSize - 1 + 0.5) * dfDeltaX;
3443 269 : for (GUInt32 nYPoint = 0; !bNeedNearest && nYPoint < nYSize; nYPoint++)
3444 : {
3445 262 : const double dfYPoint = dfYMin + (nYPoint + 0.5) * dfDeltaY;
3446 :
3447 262 : if (!GDALTriangulationFindFacetDirected(
3448 262 : psContext->sExtraParameters.psTriangulation, nStartLeft,
3449 : dfXPointMin, dfYPoint, &nStartLeft))
3450 : {
3451 6 : bNeedNearest = true;
3452 : }
3453 262 : if (!GDALTriangulationFindFacetDirected(
3454 262 : psContext->sExtraParameters.psTriangulation, nStartRight,
3455 : dfXPointMax, dfYPoint, &nStartRight))
3456 : {
3457 6 : bNeedNearest = true;
3458 : }
3459 : }
3460 7 : int nStartTop = 0;
3461 7 : int nStartBottom = 0;
3462 7 : const double dfYPointMin = dfYMin + (0 + 0.5) * dfDeltaY;
3463 7 : const double dfYPointMax = dfYMin + (nYSize - 1 + 0.5) * dfDeltaY;
3464 261 : for (GUInt32 nXPoint = 1; !bNeedNearest && nXPoint + 1 < nXSize;
3465 : nXPoint++)
3466 : {
3467 254 : const double dfXPoint = dfXMin + (nXPoint + 0.5) * dfDeltaX;
3468 :
3469 254 : if (!GDALTriangulationFindFacetDirected(
3470 254 : psContext->sExtraParameters.psTriangulation, nStartTop,
3471 : dfXPoint, dfYPointMin, &nStartTop))
3472 : {
3473 0 : bNeedNearest = true;
3474 : }
3475 254 : if (!GDALTriangulationFindFacetDirected(
3476 254 : psContext->sExtraParameters.psTriangulation, nStartBottom,
3477 : dfXPoint, dfYPointMax, &nStartBottom))
3478 : {
3479 0 : bNeedNearest = true;
3480 : }
3481 : }
3482 7 : if (bNeedNearest)
3483 : {
3484 6 : CPLDebug("GDAL_GRID", "Will need nearest neighbour");
3485 6 : GDALGridContextCreateQuadTree(psContext);
3486 : }
3487 : }
3488 :
3489 182 : int nCounter = 0;
3490 182 : volatile int bStop = FALSE;
3491 : GDALGridJob sJob;
3492 182 : sJob.nYStart = 0;
3493 182 : sJob.pabyData = static_cast<GByte *>(pData);
3494 182 : sJob.nYStep = 1;
3495 182 : sJob.nXSize = nXSize;
3496 182 : sJob.nYSize = nYSize;
3497 182 : sJob.dfXMin = dfXMin;
3498 182 : sJob.dfYMin = dfYMin;
3499 182 : sJob.dfDeltaX = dfDeltaX;
3500 182 : sJob.dfDeltaY = dfDeltaY;
3501 182 : sJob.nPoints = psContext->nPoints;
3502 182 : sJob.padfX = psContext->padfX;
3503 182 : sJob.padfY = psContext->padfY;
3504 182 : sJob.padfZ = psContext->padfZ;
3505 182 : sJob.poOptions = psContext->poOptions;
3506 182 : sJob.pfnGDALGridMethod = psContext->pfnGDALGridMethod;
3507 182 : sJob.psExtraParameters = &psContext->sExtraParameters;
3508 182 : sJob.pfnProgress = nullptr;
3509 182 : sJob.eType = eType;
3510 182 : sJob.pfnRealProgress = pfnProgress;
3511 182 : sJob.pRealProgressArg = pProgressArg;
3512 182 : sJob.pnCounter = &nCounter;
3513 182 : sJob.pbStop = &bStop;
3514 182 : sJob.hCond = nullptr;
3515 182 : sJob.hCondMutex = nullptr;
3516 :
3517 182 : if (psContext->poWorkerThreadPool == nullptr)
3518 : {
3519 2 : if (sJob.pfnRealProgress != nullptr &&
3520 2 : sJob.pfnRealProgress != GDALDummyProgress)
3521 : {
3522 2 : sJob.pfnProgress = GDALGridProgressMonoThread;
3523 : }
3524 :
3525 2 : GDALGridJobProcess(&sJob);
3526 : }
3527 : else
3528 : {
3529 180 : int nThreads = psContext->poWorkerThreadPool->GetThreadCount();
3530 : GDALGridJob *pasJobs = static_cast<GDALGridJob *>(
3531 180 : CPLMalloc(sizeof(GDALGridJob) * nThreads));
3532 :
3533 180 : sJob.nYStep = nThreads;
3534 180 : sJob.hCondMutex = CPLCreateMutex(); /* and implicitly take the mutex */
3535 180 : sJob.hCond = CPLCreateCond();
3536 180 : sJob.pfnProgress = GDALGridProgressMultiThread;
3537 :
3538 : /* --------------------------------------------------------------------
3539 : */
3540 : /* Start threads. */
3541 : /* --------------------------------------------------------------------
3542 : */
3543 900 : for (int i = 0; i < nThreads && !bStop; i++)
3544 : {
3545 720 : memcpy(&pasJobs[i], &sJob, sizeof(GDALGridJob));
3546 720 : pasJobs[i].nYStart = i;
3547 720 : psContext->poWorkerThreadPool->SubmitJob(GDALGridJobProcess,
3548 720 : &pasJobs[i]);
3549 : }
3550 :
3551 : /* --------------------------------------------------------------------
3552 : */
3553 : /* Report progress. */
3554 : /* --------------------------------------------------------------------
3555 : */
3556 10244 : while (*(sJob.pnCounter) < static_cast<int>(nYSize) && !bStop)
3557 : {
3558 10064 : CPLCondWait(sJob.hCond, sJob.hCondMutex);
3559 :
3560 10064 : int nLocalCounter = *(sJob.pnCounter);
3561 10064 : CPLReleaseMutex(sJob.hCondMutex);
3562 :
3563 20123 : if (pfnProgress != nullptr &&
3564 10059 : !pfnProgress(nLocalCounter / static_cast<double>(nYSize), "",
3565 : pProgressArg))
3566 : {
3567 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
3568 0 : bStop = TRUE;
3569 : }
3570 :
3571 10064 : CPLAcquireMutex(sJob.hCondMutex, 1.0);
3572 : }
3573 :
3574 : // Release mutex before joining threads, otherwise they will dead-lock
3575 : // forever in GDALGridProgressMultiThread().
3576 180 : CPLReleaseMutex(sJob.hCondMutex);
3577 :
3578 : /* --------------------------------------------------------------------
3579 : */
3580 : /* Wait for all threads to complete and finish. */
3581 : /* --------------------------------------------------------------------
3582 : */
3583 180 : psContext->poWorkerThreadPool->WaitCompletion();
3584 :
3585 180 : CPLFree(pasJobs);
3586 180 : CPLDestroyCond(sJob.hCond);
3587 180 : CPLDestroyMutex(sJob.hCondMutex);
3588 : }
3589 :
3590 182 : return bStop ? CE_Failure : CE_None;
3591 : }
3592 :
3593 : /************************************************************************/
3594 : /* GDALGridCreate() */
3595 : /************************************************************************/
3596 :
3597 : /**
3598 : * Create regular grid from the scattered data.
3599 : *
3600 : * This function takes the arrays of X and Y coordinates and corresponding Z
3601 : * values as input and computes regular grid (or call it a raster) from these
3602 : * scattered data. You should supply geometry and extent of the output grid
3603 : * and allocate array sufficient to hold such a grid.
3604 : *
3605 : * Starting with GDAL 1.10, it is possible to set the GDAL_NUM_THREADS
3606 : * configuration option to parallelize the processing. The value to set is
3607 : * the number of worker threads, or ALL_CPUS to use all the cores/CPUs of the
3608 : * computer (default value).
3609 : *
3610 : * Starting with GDAL 1.10, on Intel/AMD i386/x86_64 architectures, some
3611 : * gridding methods will be optimized with SSE instructions (provided GDAL
3612 : * has been compiled with such support, and it is available at runtime).
3613 : * Currently, only 'invdist' algorithm with default parameters has an optimized
3614 : * implementation.
3615 : * This can provide substantial speed-up, but sometimes at the expense of
3616 : * reduced floating point precision. This can be disabled by setting the
3617 : * GDAL_USE_SSE configuration option to NO.
3618 : * Starting with GDAL 1.11, a further optimized version can use the AVX
3619 : * instruction set. This can be disabled by setting the GDAL_USE_AVX
3620 : * configuration option to NO.
3621 : *
3622 : * Note: it will be more efficient to use GDALGridContextCreate(),
3623 : * GDALGridContextProcess() and GDALGridContextFree() when doing repeated
3624 : * gridding operations with the same algorithm, parameters and points, and
3625 : * moving the window in the output grid.
3626 : *
3627 : * @param eAlgorithm Gridding method.
3628 : * @param poOptions Options to control chosen gridding method.
3629 : * @param nPoints Number of elements in input arrays.
3630 : * @param padfX Input array of X coordinates.
3631 : * @param padfY Input array of Y coordinates.
3632 : * @param padfZ Input array of Z values.
3633 : * @param dfXMin Lowest X border of output grid.
3634 : * @param dfXMax Highest X border of output grid.
3635 : * @param dfYMin Lowest Y border of output grid.
3636 : * @param dfYMax Highest Y border of output grid.
3637 : * @param nXSize Number of columns in output grid.
3638 : * @param nYSize Number of rows in output grid.
3639 : * @param eType Data type of output array.
3640 : * @param pData Pointer to array where the computed grid will be stored.
3641 : * @param pfnProgress a GDALProgressFunc() compatible callback function for
3642 : * reporting progress or NULL.
3643 : * @param pProgressArg argument to be passed to pfnProgress. May be NULL.
3644 : *
3645 : * @return CE_None on success or CE_Failure if something goes wrong.
3646 : */
3647 :
3648 5 : CPLErr GDALGridCreate(GDALGridAlgorithm eAlgorithm, const void *poOptions,
3649 : GUInt32 nPoints, const double *padfX, const double *padfY,
3650 : const double *padfZ, double dfXMin, double dfXMax,
3651 : double dfYMin, double dfYMax, GUInt32 nXSize,
3652 : GUInt32 nYSize, GDALDataType eType, void *pData,
3653 : GDALProgressFunc pfnProgress, void *pProgressArg)
3654 : {
3655 5 : GDALGridContext *psContext = GDALGridContextCreate(
3656 : eAlgorithm, poOptions, nPoints, padfX, padfY, padfZ, TRUE);
3657 5 : CPLErr eErr = CE_Failure;
3658 5 : if (psContext)
3659 : {
3660 5 : eErr = GDALGridContextProcess(psContext, dfXMin, dfXMax, dfYMin, dfYMax,
3661 : nXSize, nYSize, eType, pData, pfnProgress,
3662 : pProgressArg);
3663 : }
3664 :
3665 5 : GDALGridContextFree(psContext);
3666 5 : return eErr;
3667 : }
3668 :
3669 : /************************************************************************/
3670 : /* GDALGridParseAlgorithmAndOptions() */
3671 : /************************************************************************/
3672 :
3673 : /** Translates mnemonic gridding algorithm names into GDALGridAlgorithm code,
3674 : * parse control parameters and assign defaults.
3675 : */
3676 381 : CPLErr GDALGridParseAlgorithmAndOptions(const char *pszAlgorithm,
3677 : GDALGridAlgorithm *peAlgorithm,
3678 : void **ppOptions)
3679 : {
3680 381 : CPLAssert(pszAlgorithm);
3681 381 : CPLAssert(peAlgorithm);
3682 381 : CPLAssert(ppOptions);
3683 :
3684 381 : *ppOptions = nullptr;
3685 :
3686 381 : char **papszParams = CSLTokenizeString2(pszAlgorithm, ":", FALSE);
3687 :
3688 381 : if (CSLCount(papszParams) < 1)
3689 : {
3690 0 : CSLDestroy(papszParams);
3691 0 : return CE_Failure;
3692 : }
3693 :
3694 381 : if (EQUAL(papszParams[0], szAlgNameInvDist))
3695 : {
3696 493 : if (CSLFetchNameValue(papszParams, "min_points_per_quadrant") ||
3697 245 : CSLFetchNameValue(papszParams, "max_points_per_quadrant"))
3698 : {
3699 : // Remap invdist to invdistnn if per quadrant is required
3700 4 : if (CSLFetchNameValue(papszParams, "radius") == nullptr)
3701 : {
3702 : const double dfRadius1 =
3703 0 : CPLAtofM(CSLFetchNameValueDef(papszParams, "radius1", "1"));
3704 : const double dfRadius2 =
3705 0 : CPLAtofM(CSLFetchNameValueDef(papszParams, "radius2", "1"));
3706 0 : if (dfRadius1 != dfRadius2)
3707 : {
3708 0 : CPLError(CE_Failure, CPLE_NotSupported,
3709 : "radius1 != radius2 not supported when "
3710 : "min_points_per_quadrant and/or "
3711 : "max_points_per_quadrant is specified");
3712 0 : CSLDestroy(papszParams);
3713 0 : return CE_Failure;
3714 : }
3715 : }
3716 :
3717 4 : if (CPLAtofM(CSLFetchNameValueDef(papszParams, "angle", "0")) != 0)
3718 : {
3719 0 : CPLError(CE_Failure, CPLE_NotSupported,
3720 : "angle != 0 not supported when "
3721 : "min_points_per_quadrant and/or "
3722 : "max_points_per_quadrant is specified");
3723 0 : CSLDestroy(papszParams);
3724 0 : return CE_Failure;
3725 : }
3726 :
3727 4 : char **papszNewParams = CSLAddString(nullptr, "invdistnn");
3728 4 : if (CSLFetchNameValue(papszParams, "radius") == nullptr)
3729 : {
3730 0 : papszNewParams = CSLSetNameValue(
3731 : papszNewParams, "radius",
3732 : CSLFetchNameValueDef(papszParams, "radius1", "1"));
3733 : }
3734 32 : for (const char *pszItem :
3735 : {"radius", "power", "smoothing", "max_points", "min_points",
3736 : "nodata", "min_points_per_quadrant",
3737 36 : "max_points_per_quadrant"})
3738 : {
3739 32 : const char *pszValue = CSLFetchNameValue(papszParams, pszItem);
3740 32 : if (pszValue)
3741 : papszNewParams =
3742 19 : CSLSetNameValue(papszNewParams, pszItem, pszValue);
3743 : }
3744 4 : CSLDestroy(papszParams);
3745 4 : papszParams = papszNewParams;
3746 :
3747 4 : *peAlgorithm = GGA_InverseDistanceToAPowerNearestNeighbor;
3748 : }
3749 : else
3750 : {
3751 244 : *peAlgorithm = GGA_InverseDistanceToAPower;
3752 : }
3753 : }
3754 133 : else if (EQUAL(papszParams[0], szAlgNameInvDistNearestNeighbor))
3755 : {
3756 18 : *peAlgorithm = GGA_InverseDistanceToAPowerNearestNeighbor;
3757 : }
3758 115 : else if (EQUAL(papszParams[0], szAlgNameAverage))
3759 : {
3760 25 : *peAlgorithm = GGA_MovingAverage;
3761 : }
3762 90 : else if (EQUAL(papszParams[0], szAlgNameNearest))
3763 : {
3764 19 : *peAlgorithm = GGA_NearestNeighbor;
3765 : }
3766 71 : else if (EQUAL(papszParams[0], szAlgNameMinimum))
3767 : {
3768 18 : *peAlgorithm = GGA_MetricMinimum;
3769 : }
3770 53 : else if (EQUAL(papszParams[0], szAlgNameMaximum))
3771 : {
3772 12 : *peAlgorithm = GGA_MetricMaximum;
3773 : }
3774 41 : else if (EQUAL(papszParams[0], szAlgNameRange))
3775 : {
3776 9 : *peAlgorithm = GGA_MetricRange;
3777 : }
3778 32 : else if (EQUAL(papszParams[0], szAlgNameCount))
3779 : {
3780 11 : *peAlgorithm = GGA_MetricCount;
3781 : }
3782 21 : else if (EQUAL(papszParams[0], szAlgNameAverageDistance))
3783 : {
3784 9 : *peAlgorithm = GGA_MetricAverageDistance;
3785 : }
3786 12 : else if (EQUAL(papszParams[0], szAlgNameAverageDistancePts))
3787 : {
3788 4 : *peAlgorithm = GGA_MetricAverageDistancePts;
3789 : }
3790 8 : else if (EQUAL(papszParams[0], szAlgNameLinear))
3791 : {
3792 7 : *peAlgorithm = GGA_Linear;
3793 : }
3794 : else
3795 : {
3796 1 : CPLError(CE_Failure, CPLE_IllegalArg,
3797 : "Unsupported gridding method \"%s\"", papszParams[0]);
3798 1 : CSLDestroy(papszParams);
3799 1 : return CE_Failure;
3800 : }
3801 :
3802 : /* -------------------------------------------------------------------- */
3803 : /* Parse algorithm parameters and assign defaults. */
3804 : /* -------------------------------------------------------------------- */
3805 380 : const char *const *papszKnownOptions = nullptr;
3806 :
3807 380 : switch (*peAlgorithm)
3808 : {
3809 244 : case GGA_InverseDistanceToAPower:
3810 : default:
3811 : {
3812 244 : *ppOptions =
3813 244 : CPLMalloc(sizeof(GDALGridInverseDistanceToAPowerOptions));
3814 :
3815 244 : GDALGridInverseDistanceToAPowerOptions *const poPowerOpts =
3816 : static_cast<GDALGridInverseDistanceToAPowerOptions *>(
3817 : *ppOptions);
3818 :
3819 244 : poPowerOpts->nSizeOfStructure = sizeof(*poPowerOpts);
3820 :
3821 244 : const char *pszValue = CSLFetchNameValue(papszParams, "power");
3822 244 : poPowerOpts->dfPower = pszValue ? CPLAtofM(pszValue) : 2.0;
3823 :
3824 244 : pszValue = CSLFetchNameValue(papszParams, "smoothing");
3825 244 : poPowerOpts->dfSmoothing = pszValue ? CPLAtofM(pszValue) : 0.0;
3826 :
3827 244 : pszValue = CSLFetchNameValue(papszParams, "radius");
3828 244 : if (pszValue)
3829 : {
3830 5 : poPowerOpts->dfRadius1 = CPLAtofM(pszValue);
3831 5 : poPowerOpts->dfRadius2 = poPowerOpts->dfRadius1;
3832 : }
3833 : else
3834 : {
3835 239 : pszValue = CSLFetchNameValue(papszParams, "radius1");
3836 239 : poPowerOpts->dfRadius1 = pszValue ? CPLAtofM(pszValue) : 0.0;
3837 :
3838 239 : pszValue = CSLFetchNameValue(papszParams, "radius2");
3839 239 : poPowerOpts->dfRadius2 = pszValue ? CPLAtofM(pszValue) : 0.0;
3840 : }
3841 :
3842 244 : pszValue = CSLFetchNameValue(papszParams, "angle");
3843 244 : poPowerOpts->dfAngle = pszValue ? CPLAtofM(pszValue) : 0.0;
3844 :
3845 244 : pszValue = CSLFetchNameValue(papszParams, "max_points");
3846 244 : poPowerOpts->nMaxPoints =
3847 244 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3848 :
3849 244 : pszValue = CSLFetchNameValue(papszParams, "min_points");
3850 244 : poPowerOpts->nMinPoints =
3851 244 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3852 :
3853 244 : pszValue = CSLFetchNameValue(papszParams, "nodata");
3854 244 : poPowerOpts->dfNoDataValue = pszValue ? CPLAtofM(pszValue) : 0.0;
3855 :
3856 : static const char *const apszKnownOptions[] = {
3857 : "power", "smoothing", "radius", "radius1", "radius2",
3858 : "angle", "max_points", "min_points", "nodata", nullptr};
3859 244 : papszKnownOptions = apszKnownOptions;
3860 :
3861 244 : break;
3862 : }
3863 22 : case GGA_InverseDistanceToAPowerNearestNeighbor:
3864 : {
3865 22 : *ppOptions = CPLMalloc(
3866 : sizeof(GDALGridInverseDistanceToAPowerNearestNeighborOptions));
3867 :
3868 : GDALGridInverseDistanceToAPowerNearestNeighborOptions
3869 22 : *const poPowerOpts = static_cast<
3870 : GDALGridInverseDistanceToAPowerNearestNeighborOptions *>(
3871 : *ppOptions);
3872 :
3873 22 : poPowerOpts->nSizeOfStructure = sizeof(*poPowerOpts);
3874 :
3875 22 : const char *pszValue = CSLFetchNameValue(papszParams, "power");
3876 22 : poPowerOpts->dfPower = pszValue ? CPLAtofM(pszValue) : 2.0;
3877 :
3878 22 : pszValue = CSLFetchNameValue(papszParams, "smoothing");
3879 22 : poPowerOpts->dfSmoothing = pszValue ? CPLAtofM(pszValue) : 0.0;
3880 :
3881 22 : pszValue = CSLFetchNameValue(papszParams, "radius");
3882 22 : poPowerOpts->dfRadius = pszValue ? CPLAtofM(pszValue) : 1.0;
3883 22 : if (!(poPowerOpts->dfRadius > 0))
3884 : {
3885 0 : CPLError(CE_Failure, CPLE_IllegalArg,
3886 : "Radius value should be strictly positive");
3887 0 : CSLDestroy(papszParams);
3888 0 : return CE_Failure;
3889 : }
3890 :
3891 22 : pszValue = CSLFetchNameValue(papszParams, "max_points");
3892 22 : poPowerOpts->nMaxPoints =
3893 22 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 12);
3894 :
3895 22 : pszValue = CSLFetchNameValue(papszParams, "min_points");
3896 22 : poPowerOpts->nMinPoints =
3897 22 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3898 :
3899 22 : pszValue = CSLFetchNameValue(papszParams, "nodata");
3900 22 : poPowerOpts->dfNoDataValue = pszValue ? CPLAtofM(pszValue) : 0.0;
3901 :
3902 : pszValue =
3903 22 : CSLFetchNameValue(papszParams, "min_points_per_quadrant");
3904 22 : poPowerOpts->nMinPointsPerQuadrant =
3905 22 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3906 :
3907 : pszValue =
3908 22 : CSLFetchNameValue(papszParams, "max_points_per_quadrant");
3909 22 : poPowerOpts->nMaxPointsPerQuadrant =
3910 22 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3911 :
3912 : static const char *const apszKnownOptions[] = {
3913 : "power",
3914 : "smoothing",
3915 : "radius",
3916 : "max_points",
3917 : "min_points",
3918 : "nodata",
3919 : "min_points_per_quadrant",
3920 : "max_points_per_quadrant",
3921 : nullptr};
3922 22 : papszKnownOptions = apszKnownOptions;
3923 :
3924 22 : break;
3925 : }
3926 25 : case GGA_MovingAverage:
3927 : {
3928 25 : *ppOptions = CPLMalloc(sizeof(GDALGridMovingAverageOptions));
3929 :
3930 25 : GDALGridMovingAverageOptions *const poAverageOpts =
3931 : static_cast<GDALGridMovingAverageOptions *>(*ppOptions);
3932 :
3933 25 : poAverageOpts->nSizeOfStructure = sizeof(*poAverageOpts);
3934 :
3935 25 : const char *pszValue = CSLFetchNameValue(papszParams, "radius");
3936 25 : if (pszValue)
3937 : {
3938 14 : poAverageOpts->dfRadius1 = CPLAtofM(pszValue);
3939 14 : poAverageOpts->dfRadius2 = poAverageOpts->dfRadius1;
3940 : }
3941 : else
3942 : {
3943 11 : pszValue = CSLFetchNameValue(papszParams, "radius1");
3944 11 : poAverageOpts->dfRadius1 = pszValue ? CPLAtofM(pszValue) : 0.0;
3945 :
3946 11 : pszValue = CSLFetchNameValue(papszParams, "radius2");
3947 11 : poAverageOpts->dfRadius2 = pszValue ? CPLAtofM(pszValue) : 0.0;
3948 : }
3949 :
3950 25 : pszValue = CSLFetchNameValue(papszParams, "angle");
3951 25 : poAverageOpts->dfAngle = pszValue ? CPLAtofM(pszValue) : 0.0;
3952 :
3953 25 : pszValue = CSLFetchNameValue(papszParams, "min_points");
3954 25 : poAverageOpts->nMinPoints =
3955 25 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3956 :
3957 25 : pszValue = CSLFetchNameValue(papszParams, "max_points");
3958 25 : poAverageOpts->nMaxPoints =
3959 25 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3960 :
3961 25 : pszValue = CSLFetchNameValue(papszParams, "nodata");
3962 25 : poAverageOpts->dfNoDataValue = pszValue ? CPLAtofM(pszValue) : 0.0;
3963 :
3964 : pszValue =
3965 25 : CSLFetchNameValue(papszParams, "min_points_per_quadrant");
3966 25 : poAverageOpts->nMinPointsPerQuadrant =
3967 25 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3968 :
3969 : pszValue =
3970 25 : CSLFetchNameValue(papszParams, "max_points_per_quadrant");
3971 25 : poAverageOpts->nMaxPointsPerQuadrant =
3972 25 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
3973 :
3974 25 : if (poAverageOpts->nMinPointsPerQuadrant != 0 ||
3975 17 : poAverageOpts->nMaxPointsPerQuadrant != 0)
3976 : {
3977 9 : if (!(poAverageOpts->dfRadius1 > 0) ||
3978 9 : !(poAverageOpts->dfRadius2 > 0))
3979 : {
3980 0 : CPLError(CE_Failure, CPLE_IllegalArg,
3981 : "Radius value should be strictly positive when "
3982 : "per quadrant parameters are specified");
3983 0 : CSLDestroy(papszParams);
3984 0 : return CE_Failure;
3985 : }
3986 9 : if (poAverageOpts->dfAngle != 0)
3987 : {
3988 0 : CPLError(CE_Failure, CPLE_NotSupported,
3989 : "angle != 0 not supported when "
3990 : "per quadrant parameters are specified");
3991 0 : CSLDestroy(papszParams);
3992 0 : return CE_Failure;
3993 : }
3994 : }
3995 16 : else if (poAverageOpts->nMaxPoints > 0)
3996 : {
3997 0 : CPLError(CE_Warning, CPLE_AppDefined,
3998 : "max_points is ignored unless one of "
3999 : "min_points_per_quadrant or max_points_per_quadrant "
4000 : "is >= 1");
4001 : }
4002 :
4003 : static const char *const apszKnownOptions[] = {
4004 : "radius",
4005 : "radius1",
4006 : "radius2",
4007 : "angle",
4008 : "min_points",
4009 : "max_points",
4010 : "nodata",
4011 : "min_points_per_quadrant",
4012 : "max_points_per_quadrant",
4013 : nullptr};
4014 25 : papszKnownOptions = apszKnownOptions;
4015 :
4016 25 : break;
4017 : }
4018 19 : case GGA_NearestNeighbor:
4019 : {
4020 19 : *ppOptions = CPLMalloc(sizeof(GDALGridNearestNeighborOptions));
4021 :
4022 19 : GDALGridNearestNeighborOptions *const poNeighborOpts =
4023 : static_cast<GDALGridNearestNeighborOptions *>(*ppOptions);
4024 :
4025 19 : poNeighborOpts->nSizeOfStructure = sizeof(*poNeighborOpts);
4026 :
4027 19 : const char *pszValue = CSLFetchNameValue(papszParams, "radius");
4028 19 : if (pszValue)
4029 : {
4030 2 : poNeighborOpts->dfRadius1 = CPLAtofM(pszValue);
4031 2 : poNeighborOpts->dfRadius2 = poNeighborOpts->dfRadius1;
4032 : }
4033 : else
4034 : {
4035 17 : pszValue = CSLFetchNameValue(papszParams, "radius1");
4036 17 : poNeighborOpts->dfRadius1 = pszValue ? CPLAtofM(pszValue) : 0.0;
4037 :
4038 17 : pszValue = CSLFetchNameValue(papszParams, "radius2");
4039 17 : poNeighborOpts->dfRadius2 = pszValue ? CPLAtofM(pszValue) : 0.0;
4040 : }
4041 :
4042 19 : pszValue = CSLFetchNameValue(papszParams, "angle");
4043 19 : poNeighborOpts->dfAngle = pszValue ? CPLAtofM(pszValue) : 0.0;
4044 :
4045 19 : pszValue = CSLFetchNameValue(papszParams, "nodata");
4046 19 : poNeighborOpts->dfNoDataValue = pszValue ? CPLAtofM(pszValue) : 0.0;
4047 :
4048 : static const char *const apszKnownOptions[] = {
4049 : "radius", "radius1", "radius2", "angle", "nodata", nullptr};
4050 19 : papszKnownOptions = apszKnownOptions;
4051 :
4052 19 : break;
4053 : }
4054 63 : case GGA_MetricMinimum:
4055 : case GGA_MetricMaximum:
4056 : case GGA_MetricRange:
4057 : case GGA_MetricCount:
4058 : case GGA_MetricAverageDistance:
4059 : case GGA_MetricAverageDistancePts:
4060 : {
4061 63 : *ppOptions = CPLMalloc(sizeof(GDALGridDataMetricsOptions));
4062 :
4063 63 : GDALGridDataMetricsOptions *const poMetricsOptions =
4064 : static_cast<GDALGridDataMetricsOptions *>(*ppOptions);
4065 :
4066 63 : poMetricsOptions->nSizeOfStructure = sizeof(*poMetricsOptions);
4067 :
4068 63 : const char *pszValue = CSLFetchNameValue(papszParams, "radius");
4069 63 : if (pszValue)
4070 : {
4071 31 : poMetricsOptions->dfRadius1 = CPLAtofM(pszValue);
4072 31 : poMetricsOptions->dfRadius2 = poMetricsOptions->dfRadius1;
4073 : }
4074 : else
4075 : {
4076 32 : pszValue = CSLFetchNameValue(papszParams, "radius1");
4077 32 : poMetricsOptions->dfRadius1 =
4078 32 : pszValue ? CPLAtofM(pszValue) : 0.0;
4079 :
4080 32 : pszValue = CSLFetchNameValue(papszParams, "radius2");
4081 32 : poMetricsOptions->dfRadius2 =
4082 32 : pszValue ? CPLAtofM(pszValue) : 0.0;
4083 : }
4084 :
4085 63 : pszValue = CSLFetchNameValue(papszParams, "angle");
4086 63 : poMetricsOptions->dfAngle = pszValue ? CPLAtofM(pszValue) : 0.0;
4087 :
4088 63 : pszValue = CSLFetchNameValue(papszParams, "min_points");
4089 63 : poMetricsOptions->nMinPoints = pszValue ? atoi(pszValue) : 0;
4090 :
4091 63 : pszValue = CSLFetchNameValue(papszParams, "nodata");
4092 63 : poMetricsOptions->dfNoDataValue =
4093 63 : pszValue ? CPLAtofM(pszValue) : 0.0;
4094 :
4095 : pszValue =
4096 63 : CSLFetchNameValue(papszParams, "min_points_per_quadrant");
4097 63 : poMetricsOptions->nMinPointsPerQuadrant =
4098 63 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
4099 :
4100 : pszValue =
4101 63 : CSLFetchNameValue(papszParams, "max_points_per_quadrant");
4102 63 : poMetricsOptions->nMaxPointsPerQuadrant =
4103 63 : static_cast<GUInt32>(pszValue ? CPLAtofM(pszValue) : 0);
4104 :
4105 63 : if (poMetricsOptions->nMinPointsPerQuadrant != 0 ||
4106 37 : poMetricsOptions->nMaxPointsPerQuadrant != 0)
4107 : {
4108 27 : if (*peAlgorithm == GGA_MetricAverageDistancePts)
4109 : {
4110 0 : CPLError(CE_Failure, CPLE_NotSupported,
4111 : "Algorithm %s not supported when "
4112 : "per quadrant parameters are specified",
4113 : szAlgNameAverageDistancePts);
4114 0 : CSLDestroy(papszParams);
4115 0 : return CE_Failure;
4116 : }
4117 27 : if (!(poMetricsOptions->dfRadius1 > 0) ||
4118 27 : !(poMetricsOptions->dfRadius2 > 0))
4119 : {
4120 0 : CPLError(CE_Failure, CPLE_IllegalArg,
4121 : "Radius value should be strictly positive when "
4122 : "per quadrant parameters are specified");
4123 0 : CSLDestroy(papszParams);
4124 0 : return CE_Failure;
4125 : }
4126 27 : if (poMetricsOptions->dfAngle != 0)
4127 : {
4128 0 : CPLError(CE_Failure, CPLE_NotSupported,
4129 : "angle != 0 not supported when "
4130 : "per quadrant parameters are specified");
4131 0 : CSLDestroy(papszParams);
4132 0 : return CE_Failure;
4133 : }
4134 : }
4135 :
4136 : static const char *const apszKnownOptions[] = {
4137 : "radius",
4138 : "radius1",
4139 : "radius2",
4140 : "angle",
4141 : "min_points",
4142 : "nodata",
4143 : "min_points_per_quadrant",
4144 : "max_points_per_quadrant",
4145 : nullptr};
4146 63 : papszKnownOptions = apszKnownOptions;
4147 :
4148 63 : break;
4149 : }
4150 7 : case GGA_Linear:
4151 : {
4152 7 : *ppOptions = CPLMalloc(sizeof(GDALGridLinearOptions));
4153 :
4154 7 : GDALGridLinearOptions *const poLinearOpts =
4155 : static_cast<GDALGridLinearOptions *>(*ppOptions);
4156 :
4157 7 : poLinearOpts->nSizeOfStructure = sizeof(*poLinearOpts);
4158 :
4159 7 : const char *pszValue = CSLFetchNameValue(papszParams, "radius");
4160 7 : poLinearOpts->dfRadius = pszValue ? CPLAtofM(pszValue) : -1.0;
4161 :
4162 7 : pszValue = CSLFetchNameValue(papszParams, "nodata");
4163 7 : poLinearOpts->dfNoDataValue = pszValue ? CPLAtofM(pszValue) : 0.0;
4164 :
4165 : static const char *const apszKnownOptions[] = {"radius", "nodata",
4166 : nullptr};
4167 7 : papszKnownOptions = apszKnownOptions;
4168 :
4169 7 : break;
4170 : }
4171 : }
4172 :
4173 380 : if (papszKnownOptions)
4174 : {
4175 1080 : for (int i = 1; papszParams[i] != nullptr; ++i)
4176 : {
4177 700 : char *pszKey = nullptr;
4178 700 : CPL_IGNORE_RET_VAL(CPLParseNameValue(papszParams[i], &pszKey));
4179 700 : if (pszKey)
4180 : {
4181 700 : bool bKnownKey = false;
4182 2981 : for (const char *const *papszKnownKeyIter = papszKnownOptions;
4183 2981 : *papszKnownKeyIter; ++papszKnownKeyIter)
4184 : {
4185 2981 : if (EQUAL(*papszKnownKeyIter, pszKey))
4186 : {
4187 700 : bKnownKey = true;
4188 700 : break;
4189 : }
4190 : }
4191 700 : if (!bKnownKey)
4192 : {
4193 0 : CPLError(CE_Warning, CPLE_AppDefined, "Option %s ignored",
4194 : pszKey);
4195 : }
4196 : }
4197 700 : CPLFree(pszKey);
4198 : }
4199 : }
4200 :
4201 380 : CSLDestroy(papszParams);
4202 380 : return CE_None;
4203 : }
|
