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