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