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