1 : /******************************************************************************
2 : * $Id: gdalgrid.cpp 25340 2012-12-21 20:30:21Z rouault $
3 : *
4 : * Project: GDAL Gridding API.
5 : * Purpose: Implementation of GDAL scattered data gridder.
6 : * Author: Andrey Kiselev, dron@ak4719.spb.edu
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2007, Andrey Kiselev <dron@ak4719.spb.edu>
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_vsi.h"
31 : #include "cpl_string.h"
32 : #include "gdalgrid.h"
33 : #include <float.h>
34 : #include <limits.h>
35 : #include "cpl_quad_tree.h"
36 : #include "cpl_multiproc.h"
37 :
38 : #ifdef HAVE_SSE_AT_COMPILE_TIME
39 : #include <xmmintrin.h>
40 : #endif
41 :
42 : CPL_CVSID("$Id: gdalgrid.cpp 25340 2012-12-21 20:30:21Z rouault $");
43 :
44 : #define TO_RADIANS (3.14159265358979323846 / 180.0)
45 :
46 : #ifndef DBL_MAX
47 : # ifdef __DBL_MAX__
48 : # define DBL_MAX __DBL_MAX__
49 : # else
50 : # define DBL_MAX 1.7976931348623157E+308
51 : # endif /* __DBL_MAX__ */
52 : #endif /* DBL_MAX */
53 :
54 : typedef struct
55 : {
56 : const double* padfX;
57 : const double* padfY;
58 : } GDALGridXYArrays;
59 :
60 : typedef struct
61 : {
62 : GDALGridXYArrays* psXYArrays;
63 : int i;
64 : } GDALGridPoint;
65 :
66 : typedef struct
67 : {
68 : CPLQuadTree* hQuadTree;
69 : double dfInitialSearchRadius;
70 : const float *pafX;
71 : const float *pafY;
72 : const float *pafZ;
73 : } GDALGridExtraParameters;
74 :
75 : /************************************************************************/
76 : /* CPLHaveRuntimeSSE() */
77 : /************************************************************************/
78 :
79 : #ifdef HAVE_SSE_AT_COMPILE_TIME
80 :
81 : #define CPUID_SSE_BIT 25
82 :
83 : #if (defined(_M_X64) || defined(__x86_64))
84 :
85 3 : static int CPLHaveRuntimeSSE()
86 : {
87 3 : return TRUE;
88 : }
89 :
90 : #elif defined(__GNUC__) && defined(__i386__)
91 :
92 : #define GCC_CPUID(level, a, b, c, d) \
93 : __asm__ ("xchgl %%ebx, %1\n" \
94 : "cpuid\n" \
95 : "xchgl %%ebx, %1" \
96 : : "=a" (a), "=r" (b), "=c" (c), "=d" (d) \
97 : : "0" (level))
98 :
99 : static int CPLHaveRuntimeSSE()
100 : {
101 : int cpuinfo[4] = {0,0,0,0};
102 : GCC_CPUID(1, cpuinfo[0], cpuinfo[1], cpuinfo[2], cpuinfo[3]);
103 : return (cpuinfo[3] & (1 << CPUID_SSE_BIT)) != 0;
104 : }
105 :
106 : #elif defined(_MSC_VER) && defined(_M_IX86)
107 :
108 : #if _MSC_VER <= 1310
109 : static void inline __cpuid(int cpuinfo[4], int level)
110 : {
111 : __asm
112 : {
113 : push ebx
114 : push esi
115 :
116 : mov esi,cpuinfo
117 : mov eax,level
118 : cpuid
119 : mov dword ptr [esi], eax
120 : mov dword ptr [esi+4],ebx
121 : mov dword ptr [esi+8],ecx
122 : mov dword ptr [esi+0Ch],edx
123 :
124 : pop esi
125 : pop ebx
126 : }
127 : }
128 : #else
129 : #include <intrin.h>
130 : #endif
131 :
132 : static int CPLHaveRuntimeSSE()
133 : {
134 : int cpuinfo[4] = {0,0,0,0};
135 : __cpuid(cpuinfo, 1);
136 : return (cpuinfo[3] & (1 << CPUID_SSE_BIT)) != 0;
137 : }
138 :
139 : #else
140 :
141 : static int CPLHaveRuntimeSSE()
142 : {
143 : return FALSE;
144 : }
145 :
146 : #endif
147 :
148 : #endif // HAVE_SSE_AT_COMPILE_TIME
149 :
150 : /************************************************************************/
151 : /* GDALGridGetPointBounds() */
152 : /************************************************************************/
153 :
154 292867 : static void GDALGridGetPointBounds(const void* hFeature, CPLRectObj* pBounds)
155 : {
156 292867 : GDALGridPoint* psPoint = (GDALGridPoint*) hFeature;
157 292867 : GDALGridXYArrays* psXYArrays = psPoint->psXYArrays;
158 292867 : int i = psPoint->i;
159 292867 : double dfX = psXYArrays->padfX[i];
160 292867 : double dfY = psXYArrays->padfY[i];
161 292867 : pBounds->minx = dfX;
162 292867 : pBounds->miny = dfY;
163 292867 : pBounds->maxx = dfX;
164 292867 : pBounds->maxy = dfY;
165 292867 : };
166 :
167 : /************************************************************************/
168 : /* GDALGridInverseDistanceToAPower() */
169 : /************************************************************************/
170 :
171 : /**
172 : * Inverse distance to a power.
173 : *
174 : * The Inverse Distance to a Power gridding method is a weighted average
175 : * interpolator. You should supply the input arrays with the scattered data
176 : * values including coordinates of every data point and output grid geometry.
177 : * The function will compute interpolated value for the given position in
178 : * output grid.
179 : *
180 : * For every grid node the resulting value \f$Z\f$ will be calculated using
181 : * formula:
182 : *
183 : * \f[
184 : * Z=\frac{\sum_{i=1}^n{\frac{Z_i}{r_i^p}}}{\sum_{i=1}^n{\frac{1}{r_i^p}}}
185 : * \f]
186 : *
187 : * where
188 : * <ul>
189 : * <li> \f$Z_i\f$ is a known value at point \f$i\f$,
190 : * <li> \f$r_i\f$ is an Euclidean distance from the grid node
191 : * to point \f$i\f$,
192 : * <li> \f$p\f$ is a weighting power,
193 : * <li> \f$n\f$ is a total number of points in search ellipse.
194 : * </ul>
195 : *
196 : * In this method the weighting factor \f$w\f$ is
197 : *
198 : * \f[
199 : * w=\frac{1}{r^p}
200 : * \f]
201 : *
202 : * @param poOptions Algorithm parameters. This should point to
203 : * GDALGridInverseDistanceToAPowerOptions object.
204 : * @param nPoints Number of elements in input arrays.
205 : * @param padfX Input array of X coordinates.
206 : * @param padfY Input array of Y coordinates.
207 : * @param padfZ Input array of Z values.
208 : * @param dfXPoint X coordinate of the point to compute.
209 : * @param dfYPoint Y coordinate of the point to compute.
210 : * @param pdfValue Pointer to variable where the computed grid node value
211 : * will be returned.
212 : * @param hExtraParamsIn extra parameters.
213 : *
214 : * @return CE_None on success or CE_Failure if something goes wrong.
215 : */
216 :
217 : CPLErr
218 11595 : GDALGridInverseDistanceToAPower( const void *poOptions, GUInt32 nPoints,
219 : const double *padfX, const double *padfY,
220 : const double *padfZ,
221 : double dfXPoint, double dfYPoint,
222 : double *pdfValue,
223 : void* hExtraParamsIn )
224 : {
225 : // TODO: For optimization purposes pre-computed parameters should be moved
226 : // out of this routine to the calling function.
227 :
228 : // Pre-compute search ellipse parameters
229 : double dfRadius1 =
230 11595 : ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->dfRadius1;
231 : double dfRadius2 =
232 11595 : ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->dfRadius2;
233 : double dfR12;
234 :
235 11595 : dfRadius1 *= dfRadius1;
236 11595 : dfRadius2 *= dfRadius2;
237 11595 : dfR12 = dfRadius1 * dfRadius2;
238 :
239 : // Compute coefficients for coordinate system rotation.
240 11595 : double dfCoeff1 = 0.0, dfCoeff2 = 0.0;
241 : const double dfAngle = TO_RADIANS
242 11595 : * ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->dfAngle;
243 11595 : const bool bRotated = ( dfAngle == 0.0 ) ? false : true;
244 11595 : if ( bRotated )
245 : {
246 0 : dfCoeff1 = cos(dfAngle);
247 0 : dfCoeff2 = sin(dfAngle);
248 : }
249 :
250 : const double dfPowerDiv2 =
251 11595 : ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->dfPower / 2;
252 : const double dfSmoothing =
253 11595 : ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->dfSmoothing;
254 : const GUInt32 nMaxPoints =
255 11595 : ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->nMaxPoints;
256 11595 : double dfNominator = 0.0, dfDenominator = 0.0;
257 11595 : GUInt32 i, n = 0;
258 :
259 131398 : for ( i = 0; i < nPoints; i++ )
260 : {
261 131000 : double dfRX = padfX[i] - dfXPoint;
262 131000 : double dfRY = padfY[i] - dfYPoint;
263 : const double dfR2 =
264 131000 : dfRX * dfRX + dfRY * dfRY + dfSmoothing * dfSmoothing;
265 :
266 131000 : if ( bRotated )
267 : {
268 0 : double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
269 0 : double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
270 :
271 0 : dfRX = dfRXRotated;
272 0 : dfRY = dfRYRotated;
273 : }
274 :
275 : // Is this point located inside the search ellipse?
276 131000 : if ( dfRadius2 * dfRX * dfRX + dfRadius1 * dfRY * dfRY <= dfR12 )
277 : {
278 : // If the test point is close to the grid node, use the point
279 : // value directly as a node value to avoid singularity.
280 14556 : if ( dfR2 < 0.0000000000001 )
281 : {
282 11197 : (*pdfValue) = padfZ[i];
283 11197 : return CE_None;
284 : }
285 : else
286 : {
287 3359 : const double dfW = pow( dfR2, dfPowerDiv2 );
288 3359 : double dfInvW = 1.0 / dfW;
289 3359 : dfNominator += dfInvW * padfZ[i];
290 3359 : dfDenominator += dfInvW;
291 3359 : n++;
292 3359 : if ( nMaxPoints > 0 && n > nMaxPoints )
293 0 : break;
294 : }
295 : }
296 : }
297 :
298 473 : if ( n < ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->nMinPoints
299 : || dfDenominator == 0.0 )
300 : {
301 : (*pdfValue) =
302 75 : ((GDALGridInverseDistanceToAPowerOptions*)poOptions)->dfNoDataValue;
303 : }
304 : else
305 323 : (*pdfValue) = dfNominator / dfDenominator;
306 :
307 398 : return CE_None;
308 : }
309 :
310 : /************************************************************************/
311 : /* GDALGridInverseDistanceToAPowerNoSearch() */
312 : /************************************************************************/
313 :
314 : /**
315 : * Inverse distance to a power for whole data set.
316 : *
317 : * This is somewhat optimized version of the Inverse Distance to a Power
318 : * method. It is used when the search ellips is not set. The algorithm and
319 : * parameters are the same as in GDALGridInverseDistanceToAPower(), but this
320 : * implementation works faster, because of no search.
321 : *
322 : * @see GDALGridInverseDistanceToAPower()
323 : */
324 :
325 : CPLErr
326 399 : GDALGridInverseDistanceToAPowerNoSearch( const void *poOptions, GUInt32 nPoints,
327 : const double *padfX, const double *padfY,
328 : const double *padfZ,
329 : double dfXPoint, double dfYPoint,
330 : double *pdfValue,
331 : void* hExtraParamsIn )
332 : {
333 : const double dfPowerDiv2 =
334 399 : ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->dfPower / 2;
335 : const double dfSmoothing =
336 399 : ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->dfSmoothing;
337 399 : const double dfSmoothing2 = dfSmoothing * dfSmoothing;
338 399 : double dfNominator = 0.0, dfDenominator = 0.0;
339 399 : GUInt32 i = 0;
340 399 : int bPower2 = (dfPowerDiv2 == 1.0);
341 :
342 399 : if( bPower2 )
343 : {
344 399 : if( dfSmoothing2 > 0 )
345 : {
346 1 : for ( i = 0; i < nPoints; i++ )
347 : {
348 0 : const double dfRX = padfX[i] - dfXPoint;
349 0 : const double dfRY = padfY[i] - dfYPoint;
350 : const double dfR2 =
351 0 : dfRX * dfRX + dfRY * dfRY + dfSmoothing2;
352 :
353 0 : double dfInvR2 = 1.0 / dfR2;
354 0 : dfNominator += dfInvR2 * padfZ[i];
355 0 : dfDenominator += dfInvR2;
356 : }
357 : }
358 : else
359 : {
360 49413 : for ( i = 0; i < nPoints; i++ )
361 : {
362 62516 : const double dfRX = padfX[i] - dfXPoint;
363 62516 : const double dfRY = padfY[i] - dfYPoint;
364 : const double dfR2 =
365 62516 : dfRX * dfRX + dfRY * dfRY;
366 :
367 : // If the test point is close to the grid node, use the point
368 : // value directly as a node value to avoid singularity.
369 62516 : if ( dfR2 < 0.0000000000001 )
370 : {
371 13501 : break;
372 : }
373 : else
374 : {
375 49015 : double dfInvR2 = 1.0 / dfR2;
376 49015 : dfNominator += dfInvR2 * padfZ[i];
377 49015 : dfDenominator += dfInvR2;
378 : }
379 : }
380 : }
381 : }
382 : else
383 : {
384 0 : for ( i = 0; i < nPoints; i++ )
385 : {
386 0 : const double dfRX = padfX[i] - dfXPoint;
387 0 : const double dfRY = padfY[i] - dfYPoint;
388 : const double dfR2 =
389 0 : dfRX * dfRX + dfRY * dfRY + dfSmoothing2;
390 :
391 : // If the test point is close to the grid node, use the point
392 : // value directly as a node value to avoid singularity.
393 0 : if ( dfR2 < 0.0000000000001 )
394 : {
395 0 : break;
396 : }
397 : else
398 : {
399 0 : const double dfW = pow( dfR2, dfPowerDiv2 );
400 0 : double dfInvW = 1.0 / dfW;
401 0 : dfNominator += dfInvW * padfZ[i];
402 0 : dfDenominator += dfInvW;
403 : }
404 : }
405 : }
406 :
407 399 : if( i != nPoints )
408 : {
409 399 : (*pdfValue) = padfZ[i];
410 : }
411 : else
412 0 : if ( dfDenominator == 0.0 )
413 : {
414 : (*pdfValue) =
415 0 : ((GDALGridInverseDistanceToAPowerOptions*)poOptions)->dfNoDataValue;
416 : }
417 : else
418 0 : (*pdfValue) = dfNominator / dfDenominator;
419 :
420 399 : return CE_None;
421 : }
422 :
423 : /************************************************************************/
424 : /* GDALGridInverseDistanceToAPower2NoSmoothingNoSearchSSE() */
425 : /************************************************************************/
426 :
427 : #ifdef HAVE_SSE_AT_COMPILE_TIME
428 :
429 : static CPLErr
430 1195 : GDALGridInverseDistanceToAPower2NoSmoothingNoSearchSSE(
431 : const void *poOptions,
432 : GUInt32 nPoints,
433 : const double *unused_padfX,
434 : const double *unused_padfY,
435 : const double *unused_padfZ,
436 : double dfXPoint, double dfYPoint,
437 : double *pdfValue,
438 : void* hExtraParamsIn )
439 : {
440 1195 : size_t i = 0;
441 1195 : GDALGridExtraParameters* psExtraParams = (GDALGridExtraParameters*) hExtraParamsIn;
442 1195 : const float* pafX = psExtraParams->pafX;
443 1195 : const float* pafY = psExtraParams->pafY;
444 1195 : const float* pafZ = psExtraParams->pafZ;
445 :
446 1195 : const float fEpsilon = 0.0000000000001f;
447 1195 : const float fXPoint = (float)dfXPoint;
448 1195 : const float fYPoint = (float)dfYPoint;
449 1195 : const __m128 xmm_small = _mm_load1_ps(&fEpsilon);
450 1195 : const __m128 xmm_x = _mm_load1_ps(&fXPoint);
451 1195 : const __m128 xmm_y = _mm_load1_ps(&fYPoint);
452 1195 : __m128 xmm_nominator = _mm_setzero_ps();
453 1195 : __m128 xmm_denominator = _mm_setzero_ps();
454 1195 : int mask = 0;
455 :
456 : #if defined(__x86_64) || defined(_M_X64)
457 : /* This would also work in 32bit mode, but there are only 8 XMM registers */
458 : /* whereas we have 16 for 64bit */
459 : #define LOOP_SIZE 8
460 1195 : size_t nPointsRound = (nPoints / LOOP_SIZE) * LOOP_SIZE;
461 29182 : for ( i = 0; i < nPointsRound; i += LOOP_SIZE )
462 : {
463 58582 : __m128 xmm_rx = _mm_sub_ps(_mm_load_ps(pafX + i), xmm_x); /* rx = pafX[i] - fXPoint */
464 58287 : __m128 xmm_rx_4 = _mm_sub_ps(_mm_load_ps(pafX + i + 4), xmm_x);
465 57721 : __m128 xmm_ry = _mm_sub_ps(_mm_load_ps(pafY + i), xmm_y); /* ry = pafY[i] - fYPoint */
466 56756 : __m128 xmm_ry_4 = _mm_sub_ps(_mm_load_ps(pafY + i + 4), xmm_y);
467 : __m128 xmm_r2 = _mm_add_ps(_mm_mul_ps(xmm_rx, xmm_rx), /* r2 = rx * rx + ry * ry */
468 28697 : _mm_mul_ps(xmm_ry, xmm_ry));
469 : __m128 xmm_r2_4 = _mm_add_ps(_mm_mul_ps(xmm_rx_4, xmm_rx_4),
470 29363 : _mm_mul_ps(xmm_ry_4, xmm_ry_4));
471 28603 : __m128 xmm_invr2 = _mm_rcp_ps(xmm_r2); /* invr2 = 1.0f / r2 */
472 28837 : __m128 xmm_invr2_4 = _mm_rcp_ps(xmm_r2_4);
473 : xmm_nominator = _mm_add_ps(xmm_nominator, /* nominator += invr2 * pafZ[i] */
474 57749 : _mm_mul_ps(xmm_invr2, _mm_load_ps(pafZ + i)));
475 : xmm_nominator = _mm_add_ps(xmm_nominator,
476 58320 : _mm_mul_ps(xmm_invr2_4, _mm_load_ps(pafZ + i + 4)));
477 29718 : xmm_denominator = _mm_add_ps(xmm_denominator, xmm_invr2); /* denominator += invr2 */
478 28933 : xmm_denominator = _mm_add_ps(xmm_denominator, xmm_invr2_4);
479 : mask = _mm_movemask_ps(_mm_cmplt_ps(xmm_r2, xmm_small)) | /* if( r2 < fEpsilon) */
480 29186 : (_mm_movemask_ps(_mm_cmplt_ps(xmm_r2_4, xmm_small)) << 4);
481 29186 : if( mask )
482 1199 : break;
483 : }
484 : #else
485 : #define LOOP_SIZE 4
486 : size_t nPointsRound = (nPoints / LOOP_SIZE) * LOOP_SIZE;
487 : for ( i = 0; i < nPointsRound; i += LOOP_SIZE )
488 : {
489 : __m128 xmm_rx = _mm_sub_ps(_mm_load_ps(pafX + i), xmm_x); /* rx = pafX[i] - fXPoint */
490 : __m128 xmm_ry = _mm_sub_ps(_mm_load_ps(pafY + i), xmm_y); /* ry = pafY[i] - fYPoint */
491 : __m128 xmm_r2 = _mm_add_ps(_mm_mul_ps(xmm_rx, xmm_rx), /* r2 = rx * rx + ry * ry */
492 : _mm_mul_ps(xmm_ry, xmm_ry));
493 : __m128 xmm_invr2 = _mm_rcp_ps(xmm_r2); /* invr2 = 1.0f / r2 */
494 : xmm_nominator = _mm_add_ps(xmm_nominator, /* nominator += invr2 * pafZ[i] */
495 : _mm_mul_ps(xmm_invr2, _mm_load_ps(pafZ + i)));
496 : xmm_denominator = _mm_add_ps(xmm_denominator, xmm_invr2); /* denominator += invr2 */
497 : mask = _mm_movemask_ps(_mm_cmplt_ps(xmm_r2, xmm_small)); /* if( r2 < fEpsilon) */
498 : if( mask )
499 : break;
500 : }
501 : #endif
502 :
503 : /* Find which i triggered r2 < fEpsilon */
504 994 : if( mask )
505 : {
506 5373 : for(int j = 0; j < LOOP_SIZE; j++ )
507 : {
508 5373 : if( mask & (1 << j) )
509 : {
510 1199 : (*pdfValue) = (pafZ)[i + j];
511 1199 : return CE_None;
512 : }
513 : }
514 : }
515 :
516 : /* Get back nominator and denominator values for XMM registers */
517 : float afNominator[4], afDenominator[4];
518 : _mm_storeu_ps(afNominator, xmm_nominator);
519 : _mm_storeu_ps(afDenominator, xmm_denominator);
520 :
521 0 : float fNominator = afNominator[0] + afNominator[1] +
522 0 : afNominator[2] + afNominator[3];
523 0 : float fDenominator = afDenominator[0] + afDenominator[1] +
524 0 : afDenominator[2] + afDenominator[3];
525 :
526 : /* Do the few remaining loop iterations */
527 0 : for ( ; i < nPoints; i++ )
528 : {
529 0 : const float fRX = pafX[i] - fXPoint;
530 0 : const float fRY = pafY[i] - fYPoint;
531 : const float fR2 =
532 0 : fRX * fRX + fRY * fRY;
533 :
534 : // If the test point is close to the grid node, use the point
535 : // value directly as a node value to avoid singularity.
536 0 : if ( fR2 < 0.0000000000001 )
537 : {
538 0 : break;
539 : }
540 : else
541 : {
542 0 : const float fInvR2 = 1.0f / fR2;
543 0 : fNominator += fInvR2 * pafZ[i];
544 0 : fDenominator += fInvR2;
545 : }
546 : }
547 :
548 0 : if( i != nPoints )
549 : {
550 0 : (*pdfValue) = pafZ[i];
551 : }
552 : else
553 0 : if ( fDenominator == 0.0 )
554 : {
555 : (*pdfValue) =
556 0 : ((GDALGridInverseDistanceToAPowerOptions*)poOptions)->dfNoDataValue;
557 : }
558 : else
559 0 : (*pdfValue) = fNominator / fDenominator;
560 :
561 0 : return CE_None;
562 : }
563 : #endif // HAVE_SSE_AT_COMPILE_TIME
564 :
565 : /************************************************************************/
566 : /* GDALGridMovingAverage() */
567 : /************************************************************************/
568 :
569 : /**
570 : * Moving average.
571 : *
572 : * The Moving Average is a simple data averaging algorithm. It uses a moving
573 : * window of elliptic form to search values and averages all data points
574 : * within the window. Search ellipse can be rotated by specified angle, the
575 : * center of ellipse located at the grid node. Also the minimum number of data
576 : * points to average can be set, if there are not enough points in window, the
577 : * grid node considered empty and will be filled with specified NODATA value.
578 : *
579 : * Mathematically it can be expressed with the formula:
580 : *
581 : * \f[
582 : * Z=\frac{\sum_{i=1}^n{Z_i}}{n}
583 : * \f]
584 : *
585 : * where
586 : * <ul>
587 : * <li> \f$Z\f$ is a resulting value at the grid node,
588 : * <li> \f$Z_i\f$ is a known value at point \f$i\f$,
589 : * <li> \f$n\f$ is a total number of points in search ellipse.
590 : * </ul>
591 : *
592 : * @param poOptions Algorithm parameters. This should point to
593 : * GDALGridMovingAverageOptions object.
594 : * @param nPoints Number of elements in input arrays.
595 : * @param padfX Input array of X coordinates.
596 : * @param padfY Input array of Y coordinates.
597 : * @param padfZ Input array of Z values.
598 : * @param dfXPoint X coordinate of the point to compute.
599 : * @param dfYPoint Y coordinate of the point to compute.
600 : * @param pdfValue Pointer to variable where the computed grid node value
601 : * will be returned.
602 : *
603 : * @return CE_None on success or CE_Failure if something goes wrong.
604 : */
605 :
606 : CPLErr
607 1589 : GDALGridMovingAverage( const void *poOptions, GUInt32 nPoints,
608 : const double *padfX, const double *padfY,
609 : const double *padfZ,
610 : double dfXPoint, double dfYPoint, double *pdfValue,
611 : void* hExtraParamsIn )
612 : {
613 : // TODO: For optimization purposes pre-computed parameters should be moved
614 : // out of this routine to the calling function.
615 :
616 : // Pre-compute search ellipse parameters
617 1589 : double dfRadius1 = ((GDALGridMovingAverageOptions *)poOptions)->dfRadius1;
618 1589 : double dfRadius2 = ((GDALGridMovingAverageOptions *)poOptions)->dfRadius2;
619 : double dfR12;
620 :
621 1589 : dfRadius1 *= dfRadius1;
622 1589 : dfRadius2 *= dfRadius2;
623 1589 : dfR12 = dfRadius1 * dfRadius2;
624 :
625 : // Compute coefficients for coordinate system rotation.
626 1589 : double dfCoeff1 = 0.0, dfCoeff2 = 0.0;
627 : const double dfAngle =
628 1589 : TO_RADIANS * ((GDALGridMovingAverageOptions *)poOptions)->dfAngle;
629 1589 : const bool bRotated = ( dfAngle == 0.0 ) ? false : true;
630 1589 : if ( bRotated )
631 : {
632 400 : dfCoeff1 = cos(dfAngle);
633 400 : dfCoeff2 = sin(dfAngle);
634 : }
635 :
636 1589 : double dfAccumulator = 0.0;
637 1589 : GUInt32 i = 0, n = 0;
638 :
639 420158 : while ( i < nPoints )
640 : {
641 416980 : double dfRX = padfX[i] - dfXPoint;
642 416980 : double dfRY = padfY[i] - dfYPoint;
643 :
644 416980 : if ( bRotated )
645 : {
646 116694 : double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
647 116694 : double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
648 :
649 116694 : dfRX = dfRXRotated;
650 116694 : dfRY = dfRYRotated;
651 : }
652 :
653 : // Is this point located inside the search ellipse?
654 416980 : if ( dfRadius2 * dfRX * dfRX + dfRadius1 * dfRY * dfRY <= dfR12 )
655 : {
656 138478 : dfAccumulator += padfZ[i];
657 138478 : n++;
658 : }
659 :
660 416980 : i++;
661 : }
662 :
663 1665 : if ( n < ((GDALGridMovingAverageOptions *)poOptions)->nMinPoints
664 : || n == 0 )
665 : {
666 : (*pdfValue) =
667 76 : ((GDALGridMovingAverageOptions *)poOptions)->dfNoDataValue;
668 : }
669 : else
670 1513 : (*pdfValue) = dfAccumulator / n;
671 :
672 1589 : return CE_None;
673 : }
674 :
675 : /************************************************************************/
676 : /* GDALGridNearestNeighbor() */
677 : /************************************************************************/
678 :
679 : /**
680 : * Nearest neighbor.
681 : *
682 : * The Nearest Neighbor method doesn't perform any interpolation or smoothing,
683 : * it just takes the value of nearest point found in grid node search ellipse
684 : * and returns it as a result. If there are no points found, the specified
685 : * NODATA value will be returned.
686 : *
687 : * @param poOptions Algorithm parameters. This should point to
688 : * GDALGridNearestNeighborOptions object.
689 : * @param nPoints Number of elements in input arrays.
690 : * @param padfX Input array of X coordinates.
691 : * @param padfY Input array of Y coordinates.
692 : * @param padfZ Input array of Z values.
693 : * @param dfXPoint X coordinate of the point to compute.
694 : * @param dfYPoint Y coordinate of the point to compute.
695 : * @param pdfValue Pointer to variable where the computed grid node value
696 : * will be returned.
697 : *
698 : * @return CE_None on success or CE_Failure if something goes wrong.
699 : */
700 :
701 : CPLErr
702 17007 : GDALGridNearestNeighbor( const void *poOptions, GUInt32 nPoints,
703 : const double *padfX, const double *padfY,
704 : const double *padfZ,
705 : double dfXPoint, double dfYPoint, double *pdfValue,
706 : void* hExtraParamsIn)
707 : {
708 : // TODO: For optimization purposes pre-computed parameters should be moved
709 : // out of this routine to the calling function.
710 :
711 : // Pre-compute search ellipse parameters
712 : double dfRadius1 =
713 17007 : ((GDALGridNearestNeighborOptions *)poOptions)->dfRadius1;
714 : double dfRadius2 =
715 17007 : ((GDALGridNearestNeighborOptions *)poOptions)->dfRadius2;
716 : double dfR12;
717 17007 : GDALGridExtraParameters* psExtraParams = (GDALGridExtraParameters*) hExtraParamsIn;
718 17007 : CPLQuadTree* hQuadTree = psExtraParams->hQuadTree;
719 :
720 17007 : dfRadius1 *= dfRadius1;
721 17007 : dfRadius2 *= dfRadius2;
722 17007 : dfR12 = dfRadius1 * dfRadius2;
723 :
724 : // Compute coefficients for coordinate system rotation.
725 17007 : double dfCoeff1 = 0.0, dfCoeff2 = 0.0;
726 : const double dfAngle =
727 17007 : TO_RADIANS * ((GDALGridNearestNeighborOptions *)poOptions)->dfAngle;
728 17007 : const bool bRotated = ( dfAngle == 0.0 ) ? false : true;
729 17007 : if ( bRotated )
730 : {
731 0 : dfCoeff1 = cos(dfAngle);
732 0 : dfCoeff2 = sin(dfAngle);
733 : }
734 :
735 : // If the nearest point will not be found, its value remains as NODATA.
736 : double dfNearestValue =
737 17007 : ((GDALGridNearestNeighborOptions *)poOptions)->dfNoDataValue;
738 : // Nearest distance will be initialized with the distance to the first
739 : // point in array.
740 17007 : double dfNearestR = DBL_MAX;
741 17007 : GUInt32 i = 0;
742 :
743 17007 : double dfSearchRadius = psExtraParams->dfInitialSearchRadius;
744 34033 : if( hQuadTree != NULL && dfRadius1 == dfRadius2 && dfSearchRadius > 0 )
745 : {
746 : CPLRectObj sAoi;
747 17007 : if( dfRadius1 > 0 )
748 1598 : dfSearchRadius = ((GDALGridNearestNeighborOptions *)poOptions)->dfRadius1;
749 34014 : while(dfSearchRadius > 0)
750 : {
751 17007 : sAoi.minx = dfXPoint - dfSearchRadius;
752 17007 : sAoi.miny = dfYPoint - dfSearchRadius;
753 17007 : sAoi.maxx = dfXPoint + dfSearchRadius;
754 17007 : sAoi.maxy = dfYPoint + dfSearchRadius;
755 17007 : int nFeatureCount = 0;
756 : GDALGridPoint** papsPoints = (GDALGridPoint**)
757 17007 : CPLQuadTreeSearch(hQuadTree, &sAoi, &nFeatureCount);
758 16365 : if( nFeatureCount != 0 )
759 : {
760 16365 : if( dfRadius1 > 0 )
761 1542 : dfNearestR = dfRadius1;
762 59720 : for(int k=0; k<nFeatureCount; k++)
763 : {
764 43355 : int i = papsPoints[k]->i;
765 43355 : double dfRX = padfX[i] - dfXPoint;
766 43355 : double dfRY = padfY[i] - dfYPoint;
767 :
768 43355 : const double dfR2 = dfRX * dfRX + dfRY * dfRY;
769 43355 : if( dfR2 <= dfNearestR )
770 : {
771 20133 : dfNearestR = dfR2;
772 20133 : dfNearestValue = padfZ[i];
773 : }
774 : }
775 :
776 16365 : CPLFree(papsPoints);
777 17026 : break;
778 : }
779 : else
780 : {
781 0 : CPLFree(papsPoints);
782 0 : if( dfRadius1 > 0 )
783 0 : break;
784 0 : dfSearchRadius *= 2;
785 : //CPLDebug("GDAL_GRID", "Increasing search radius to %.16g", dfSearchRadius);
786 : }
787 : }
788 : }
789 : else
790 : {
791 0 : while ( i < nPoints )
792 : {
793 0 : double dfRX = padfX[i] - dfXPoint;
794 0 : double dfRY = padfY[i] - dfYPoint;
795 :
796 0 : if ( bRotated )
797 : {
798 0 : double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
799 0 : double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
800 :
801 0 : dfRX = dfRXRotated;
802 0 : dfRY = dfRYRotated;
803 : }
804 :
805 : // Is this point located inside the search ellipse?
806 0 : if ( dfRadius2 * dfRX * dfRX + dfRadius1 * dfRY * dfRY <= dfR12 )
807 : {
808 0 : const double dfR2 = dfRX * dfRX + dfRY * dfRY;
809 0 : if ( dfR2 <= dfNearestR )
810 : {
811 0 : dfNearestR = dfR2;
812 0 : dfNearestValue = padfZ[i];
813 : }
814 : }
815 :
816 0 : i++;
817 : }
818 : }
819 :
820 17026 : (*pdfValue) = dfNearestValue;
821 :
822 17026 : return CE_None;
823 : }
824 :
825 : /************************************************************************/
826 : /* GDALGridDataMetricMinimum() */
827 : /************************************************************************/
828 :
829 : /**
830 : * Minimum data value (data metric).
831 : *
832 : * Minimum value found in grid node search ellipse. If there are no points
833 : * found, the specified NODATA value will be returned.
834 : *
835 : * \f[
836 : * Z=\min{(Z_1,Z_2,\ldots,Z_n)}
837 : * \f]
838 : *
839 : * where
840 : * <ul>
841 : * <li> \f$Z\f$ is a resulting value at the grid node,
842 : * <li> \f$Z_i\f$ is a known value at point \f$i\f$,
843 : * <li> \f$n\f$ is a total number of points in search ellipse.
844 : * </ul>
845 : *
846 : * @param poOptions Algorithm parameters. This should point to
847 : * GDALGridDataMetricsOptions object.
848 : * @param nPoints Number of elements in input arrays.
849 : * @param padfX Input array of X coordinates.
850 : * @param padfY Input array of Y coordinates.
851 : * @param padfZ Input array of Z values.
852 : * @param dfXPoint X coordinate of the point to compute.
853 : * @param dfYPoint Y coordinate of the point to compute.
854 : * @param pdfValue Pointer to variable where the computed grid node value
855 : * will be returned.
856 : *
857 : * @return CE_None on success or CE_Failure if something goes wrong.
858 : */
859 :
860 : CPLErr
861 800 : GDALGridDataMetricMinimum( const void *poOptions, GUInt32 nPoints,
862 : const double *padfX, const double *padfY,
863 : const double *padfZ,
864 : double dfXPoint, double dfYPoint, double *pdfValue,
865 : void* hExtraParamsIn )
866 : {
867 : // TODO: For optimization purposes pre-computed parameters should be moved
868 : // out of this routine to the calling function.
869 :
870 : // Pre-compute search ellipse parameters
871 : double dfRadius1 =
872 800 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius1;
873 : double dfRadius2 =
874 800 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius2;
875 : double dfR12;
876 :
877 800 : dfRadius1 *= dfRadius1;
878 800 : dfRadius2 *= dfRadius2;
879 800 : dfR12 = dfRadius1 * dfRadius2;
880 :
881 : // Compute coefficients for coordinate system rotation.
882 800 : double dfCoeff1 = 0.0, dfCoeff2 = 0.0;
883 : const double dfAngle =
884 800 : TO_RADIANS * ((GDALGridDataMetricsOptions *)poOptions)->dfAngle;
885 800 : const bool bRotated = ( dfAngle == 0.0 ) ? false : true;
886 800 : if ( bRotated )
887 : {
888 400 : dfCoeff1 = cos(dfAngle);
889 400 : dfCoeff2 = sin(dfAngle);
890 : }
891 :
892 800 : double dfMinimumValue=0.0;
893 800 : GUInt32 i = 0, n = 0;
894 :
895 296275 : while ( i < nPoints )
896 : {
897 294675 : double dfRX = padfX[i] - dfXPoint;
898 294675 : double dfRY = padfY[i] - dfYPoint;
899 :
900 294675 : if ( bRotated )
901 : {
902 130657 : double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
903 130657 : double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
904 :
905 130657 : dfRX = dfRXRotated;
906 130657 : dfRY = dfRYRotated;
907 : }
908 :
909 : // Is this point located inside the search ellipse?
910 294675 : if ( dfRadius2 * dfRX * dfRX + dfRadius1 * dfRY * dfRY <= dfR12 )
911 : {
912 140833 : if ( n )
913 : {
914 140034 : if ( dfMinimumValue > padfZ[i] )
915 2069 : dfMinimumValue = padfZ[i];
916 : }
917 : else
918 799 : dfMinimumValue = padfZ[i];
919 140833 : n++;
920 : }
921 :
922 294675 : i++;
923 : }
924 :
925 800 : if ( n < ((GDALGridDataMetricsOptions *)poOptions)->nMinPoints
926 : || n == 0 )
927 : {
928 : (*pdfValue) =
929 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfNoDataValue;
930 : }
931 : else
932 800 : (*pdfValue) = dfMinimumValue;
933 :
934 800 : return CE_None;
935 : }
936 :
937 : /************************************************************************/
938 : /* GDALGridDataMetricMaximum() */
939 : /************************************************************************/
940 :
941 : /**
942 : * Maximum data value (data metric).
943 : *
944 : * Maximum value found in grid node search ellipse. If there are no points
945 : * found, the specified NODATA value will be returned.
946 : *
947 : * \f[
948 : * Z=\max{(Z_1,Z_2,\ldots,Z_n)}
949 : * \f]
950 : *
951 : * where
952 : * <ul>
953 : * <li> \f$Z\f$ is a resulting value at the grid node,
954 : * <li> \f$Z_i\f$ is a known value at point \f$i\f$,
955 : * <li> \f$n\f$ is a total number of points in search ellipse.
956 : * </ul>
957 : *
958 : * @param poOptions Algorithm parameters. This should point to
959 : * GDALGridDataMetricsOptions object.
960 : * @param nPoints Number of elements in input arrays.
961 : * @param padfX Input array of X coordinates.
962 : * @param padfY Input array of Y coordinates.
963 : * @param padfZ Input array of Z values.
964 : * @param dfXPoint X coordinate of the point to compute.
965 : * @param dfYPoint Y coordinate of the point to compute.
966 : * @param pdfValue Pointer to variable where the computed grid node value
967 : * will be returned.
968 : *
969 : * @return CE_None on success or CE_Failure if something goes wrong.
970 : */
971 :
972 : CPLErr
973 795 : GDALGridDataMetricMaximum( const void *poOptions, GUInt32 nPoints,
974 : const double *padfX, const double *padfY,
975 : const double *padfZ,
976 : double dfXPoint, double dfYPoint, double *pdfValue,
977 : void* hExtraParamsIn )
978 : {
979 : // TODO: For optimization purposes pre-computed parameters should be moved
980 : // out of this routine to the calling function.
981 :
982 : // Pre-compute search ellipse parameters
983 : double dfRadius1 =
984 795 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius1;
985 : double dfRadius2 =
986 795 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius2;
987 : double dfR12;
988 :
989 795 : dfRadius1 *= dfRadius1;
990 795 : dfRadius2 *= dfRadius2;
991 795 : dfR12 = dfRadius1 * dfRadius2;
992 :
993 : // Compute coefficients for coordinate system rotation.
994 795 : double dfCoeff1 = 0.0, dfCoeff2 = 0.0;
995 : const double dfAngle =
996 795 : TO_RADIANS * ((GDALGridDataMetricsOptions *)poOptions)->dfAngle;
997 795 : const bool bRotated = ( dfAngle == 0.0 ) ? false : true;
998 795 : if ( bRotated )
999 : {
1000 0 : dfCoeff1 = cos(dfAngle);
1001 0 : dfCoeff2 = sin(dfAngle);
1002 : }
1003 :
1004 795 : double dfMaximumValue=0.0;
1005 795 : GUInt32 i = 0, n = 0;
1006 :
1007 255408 : while ( i < nPoints )
1008 : {
1009 253818 : double dfRX = padfX[i] - dfXPoint;
1010 253818 : double dfRY = padfY[i] - dfYPoint;
1011 :
1012 253818 : if ( bRotated )
1013 : {
1014 0 : double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
1015 0 : double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
1016 :
1017 0 : dfRX = dfRXRotated;
1018 0 : dfRY = dfRYRotated;
1019 : }
1020 :
1021 : // Is this point located inside the search ellipse?
1022 253818 : if ( dfRadius2 * dfRX * dfRX + dfRadius1 * dfRY * dfRY <= dfR12 )
1023 : {
1024 143166 : if ( n )
1025 : {
1026 142371 : if ( dfMaximumValue < padfZ[i] )
1027 4074 : dfMaximumValue = padfZ[i];
1028 : }
1029 : else
1030 795 : dfMaximumValue = padfZ[i];
1031 143166 : n++;
1032 : }
1033 :
1034 253818 : i++;
1035 : }
1036 :
1037 795 : if ( n < ((GDALGridDataMetricsOptions *)poOptions)->nMinPoints
1038 : || n == 0 )
1039 : {
1040 : (*pdfValue) =
1041 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfNoDataValue;
1042 : }
1043 : else
1044 795 : (*pdfValue) = dfMaximumValue;
1045 :
1046 795 : return CE_None;
1047 : }
1048 :
1049 : /************************************************************************/
1050 : /* GDALGridDataMetricRange() */
1051 : /************************************************************************/
1052 :
1053 : /**
1054 : * Data range (data metric).
1055 : *
1056 : * A difference between the minimum and maximum values found in grid node
1057 : * search ellipse. If there are no points found, the specified NODATA
1058 : * value will be returned.
1059 : *
1060 : * \f[
1061 : * Z=\max{(Z_1,Z_2,\ldots,Z_n)}-\min{(Z_1,Z_2,\ldots,Z_n)}
1062 : * \f]
1063 : *
1064 : * where
1065 : * <ul>
1066 : * <li> \f$Z\f$ is a resulting value at the grid node,
1067 : * <li> \f$Z_i\f$ is a known value at point \f$i\f$,
1068 : * <li> \f$n\f$ is a total number of points in search ellipse.
1069 : * </ul>
1070 : *
1071 : * @param poOptions Algorithm parameters. This should point to
1072 : * GDALGridDataMetricsOptions object.
1073 : * @param nPoints Number of elements in input arrays.
1074 : * @param padfX Input array of X coordinates.
1075 : * @param padfY Input array of Y coordinates.
1076 : * @param padfZ Input array of Z values.
1077 : * @param dfXPoint X coordinate of the point to compute.
1078 : * @param dfYPoint Y coordinate of the point to compute.
1079 : * @param pdfValue Pointer to variable where the computed grid node value
1080 : * will be returned.
1081 : *
1082 : * @return CE_None on success or CE_Failure if something goes wrong.
1083 : */
1084 :
1085 : CPLErr
1086 800 : GDALGridDataMetricRange( const void *poOptions, GUInt32 nPoints,
1087 : const double *padfX, const double *padfY,
1088 : const double *padfZ,
1089 : double dfXPoint, double dfYPoint, double *pdfValue,
1090 : void* hExtraParamsIn )
1091 : {
1092 : // TODO: For optimization purposes pre-computed parameters should be moved
1093 : // out of this routine to the calling function.
1094 :
1095 : // Pre-compute search ellipse parameters
1096 : double dfRadius1 =
1097 800 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius1;
1098 : double dfRadius2 =
1099 800 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius2;
1100 : double dfR12;
1101 :
1102 800 : dfRadius1 *= dfRadius1;
1103 800 : dfRadius2 *= dfRadius2;
1104 800 : dfR12 = dfRadius1 * dfRadius2;
1105 :
1106 : // Compute coefficients for coordinate system rotation.
1107 800 : double dfCoeff1 = 0.0, dfCoeff2 = 0.0;
1108 : const double dfAngle =
1109 800 : TO_RADIANS * ((GDALGridDataMetricsOptions *)poOptions)->dfAngle;
1110 800 : const bool bRotated = ( dfAngle == 0.0 ) ? false : true;
1111 800 : if ( bRotated )
1112 : {
1113 0 : dfCoeff1 = cos(dfAngle);
1114 0 : dfCoeff2 = sin(dfAngle);
1115 : }
1116 :
1117 800 : double dfMaximumValue=0.0, dfMinimumValue=0.0;
1118 800 : GUInt32 i = 0, n = 0;
1119 :
1120 233685 : while ( i < nPoints )
1121 : {
1122 232085 : double dfRX = padfX[i] - dfXPoint;
1123 232085 : double dfRY = padfY[i] - dfYPoint;
1124 :
1125 232085 : if ( bRotated )
1126 : {
1127 0 : double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
1128 0 : double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
1129 :
1130 0 : dfRX = dfRXRotated;
1131 0 : dfRY = dfRYRotated;
1132 : }
1133 :
1134 : // Is this point located inside the search ellipse?
1135 232085 : if ( dfRadius2 * dfRX * dfRX + dfRadius1 * dfRY * dfRY <= dfR12 )
1136 : {
1137 134265 : if ( n )
1138 : {
1139 133467 : if ( dfMinimumValue > padfZ[i] )
1140 1746 : dfMinimumValue = padfZ[i];
1141 133467 : if ( dfMaximumValue < padfZ[i] )
1142 4137 : dfMaximumValue = padfZ[i];
1143 : }
1144 : else
1145 798 : dfMinimumValue = dfMaximumValue = padfZ[i];
1146 134265 : n++;
1147 : }
1148 :
1149 232085 : i++;
1150 : }
1151 :
1152 876 : if ( n < ((GDALGridDataMetricsOptions *)poOptions)->nMinPoints
1153 : || n == 0 )
1154 : {
1155 : (*pdfValue) =
1156 76 : ((GDALGridDataMetricsOptions *)poOptions)->dfNoDataValue;
1157 : }
1158 : else
1159 724 : (*pdfValue) = dfMaximumValue - dfMinimumValue;
1160 :
1161 800 : return CE_None;
1162 : }
1163 :
1164 : /************************************************************************/
1165 : /* GDALGridDataMetricCount() */
1166 : /************************************************************************/
1167 :
1168 : /**
1169 : * Number of data points (data metric).
1170 : *
1171 : * A number of data points found in grid node search ellipse.
1172 : *
1173 : * \f[
1174 : * Z=n
1175 : * \f]
1176 : *
1177 : * where
1178 : * <ul>
1179 : * <li> \f$Z\f$ is a resulting value at the grid node,
1180 : * <li> \f$n\f$ is a total number of points in search ellipse.
1181 : * </ul>
1182 : *
1183 : * @param poOptions Algorithm parameters. This should point to
1184 : * GDALGridDataMetricsOptions object.
1185 : * @param nPoints Number of elements in input arrays.
1186 : * @param padfX Input array of X coordinates.
1187 : * @param padfY Input array of Y coordinates.
1188 : * @param padfZ Input array of Z values.
1189 : * @param dfXPoint X coordinate of the point to compute.
1190 : * @param dfYPoint Y coordinate of the point to compute.
1191 : * @param pdfValue Pointer to variable where the computed grid node value
1192 : * will be returned.
1193 : *
1194 : * @return CE_None on success or CE_Failure if something goes wrong.
1195 : */
1196 :
1197 : CPLErr
1198 798 : GDALGridDataMetricCount( const void *poOptions, GUInt32 nPoints,
1199 : const double *padfX, const double *padfY,
1200 : const double *padfZ,
1201 : double dfXPoint, double dfYPoint, double *pdfValue,
1202 : void* hExtraParamsIn )
1203 : {
1204 : // TODO: For optimization purposes pre-computed parameters should be moved
1205 : // out of this routine to the calling function.
1206 :
1207 : // Pre-compute search ellipse parameters
1208 : double dfRadius1 =
1209 798 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius1;
1210 : double dfRadius2 =
1211 798 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius2;
1212 : double dfR12;
1213 :
1214 798 : dfRadius1 *= dfRadius1;
1215 798 : dfRadius2 *= dfRadius2;
1216 798 : dfR12 = dfRadius1 * dfRadius2;
1217 :
1218 : // Compute coefficients for coordinate system rotation.
1219 798 : double dfCoeff1 = 0.0, dfCoeff2 = 0.0;
1220 : const double dfAngle =
1221 798 : TO_RADIANS * ((GDALGridDataMetricsOptions *)poOptions)->dfAngle;
1222 798 : const bool bRotated = ( dfAngle == 0.0 ) ? false : true;
1223 798 : if ( bRotated )
1224 : {
1225 0 : dfCoeff1 = cos(dfAngle);
1226 0 : dfCoeff2 = sin(dfAngle);
1227 : }
1228 :
1229 798 : GUInt32 i = 0, n = 0;
1230 :
1231 150694 : while ( i < nPoints )
1232 : {
1233 149098 : double dfRX = padfX[i] - dfXPoint;
1234 149098 : double dfRY = padfY[i] - dfYPoint;
1235 :
1236 149098 : if ( bRotated )
1237 : {
1238 0 : double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
1239 0 : double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
1240 :
1241 0 : dfRX = dfRXRotated;
1242 0 : dfRY = dfRYRotated;
1243 : }
1244 :
1245 : // Is this point located inside the search ellipse?
1246 149098 : if ( dfRadius2 * dfRX * dfRX + dfRadius1 * dfRY * dfRY <= dfR12 )
1247 22733 : n++;
1248 :
1249 149098 : i++;
1250 : }
1251 :
1252 798 : if ( n < ((GDALGridDataMetricsOptions *)poOptions)->nMinPoints )
1253 : {
1254 : (*pdfValue) =
1255 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfNoDataValue;
1256 : }
1257 : else
1258 798 : (*pdfValue) = (double)n;
1259 :
1260 798 : return CE_None;
1261 : }
1262 :
1263 : /************************************************************************/
1264 : /* GDALGridDataMetricAverageDistance() */
1265 : /************************************************************************/
1266 :
1267 : /**
1268 : * Average distance (data metric).
1269 : *
1270 : * An average distance between the grid node (center of the search ellipse)
1271 : * and all of the data points found in grid node search ellipse. If there are
1272 : * no points found, the specified NODATA value will be returned.
1273 : *
1274 : * \f[
1275 : * Z=\frac{\sum_{i = 1}^n r_i}{n}
1276 : * \f]
1277 : *
1278 : * where
1279 : * <ul>
1280 : * <li> \f$Z\f$ is a resulting value at the grid node,
1281 : * <li> \f$r_i\f$ is an Euclidean distance from the grid node
1282 : * to point \f$i\f$,
1283 : * <li> \f$n\f$ is a total number of points in search ellipse.
1284 : * </ul>
1285 : *
1286 : * @param poOptions Algorithm parameters. This should point to
1287 : * GDALGridDataMetricsOptions object.
1288 : * @param nPoints Number of elements in input arrays.
1289 : * @param padfX Input array of X coordinates.
1290 : * @param padfY Input array of Y coordinates.
1291 : * @param padfZ Input array of Z values.
1292 : * @param dfXPoint X coordinate of the point to compute.
1293 : * @param dfYPoint Y coordinate of the point to compute.
1294 : * @param pdfValue Pointer to variable where the computed grid node value
1295 : * will be returned.
1296 : *
1297 : * @return CE_None on success or CE_Failure if something goes wrong.
1298 : */
1299 :
1300 : CPLErr
1301 798 : GDALGridDataMetricAverageDistance( const void *poOptions, GUInt32 nPoints,
1302 : const double *padfX, const double *padfY,
1303 : const double *padfZ,
1304 : double dfXPoint, double dfYPoint,
1305 : double *pdfValue,
1306 : void* hExtraParamsIn )
1307 : {
1308 : // TODO: For optimization purposes pre-computed parameters should be moved
1309 : // out of this routine to the calling function.
1310 :
1311 : // Pre-compute search ellipse parameters
1312 : double dfRadius1 =
1313 798 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius1;
1314 : double dfRadius2 =
1315 798 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius2;
1316 : double dfR12;
1317 :
1318 798 : dfRadius1 *= dfRadius1;
1319 798 : dfRadius2 *= dfRadius2;
1320 798 : dfR12 = dfRadius1 * dfRadius2;
1321 :
1322 : // Compute coefficients for coordinate system rotation.
1323 798 : double dfCoeff1 = 0.0, dfCoeff2 = 0.0;
1324 : const double dfAngle =
1325 798 : TO_RADIANS * ((GDALGridDataMetricsOptions *)poOptions)->dfAngle;
1326 798 : const bool bRotated = ( dfAngle == 0.0 ) ? false : true;
1327 798 : if ( bRotated )
1328 : {
1329 0 : dfCoeff1 = cos(dfAngle);
1330 0 : dfCoeff2 = sin(dfAngle);
1331 : }
1332 :
1333 798 : double dfAccumulator = 0.0;
1334 798 : GUInt32 i = 0, n = 0;
1335 :
1336 232274 : while ( i < nPoints )
1337 : {
1338 230678 : double dfRX = padfX[i] - dfXPoint;
1339 230678 : double dfRY = padfY[i] - dfYPoint;
1340 :
1341 230678 : if ( bRotated )
1342 : {
1343 0 : double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
1344 0 : double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
1345 :
1346 0 : dfRX = dfRXRotated;
1347 0 : dfRY = dfRYRotated;
1348 : }
1349 :
1350 : // Is this point located inside the search ellipse?
1351 230678 : if ( dfRadius2 * dfRX * dfRX + dfRadius1 * dfRY * dfRY <= dfR12 )
1352 : {
1353 142793 : dfAccumulator += sqrt( dfRX * dfRX + dfRY * dfRY );
1354 142793 : n++;
1355 : }
1356 :
1357 230678 : i++;
1358 : }
1359 :
1360 798 : if ( n < ((GDALGridDataMetricsOptions *)poOptions)->nMinPoints
1361 : || n == 0 )
1362 : {
1363 : (*pdfValue) =
1364 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfNoDataValue;
1365 : }
1366 : else
1367 798 : (*pdfValue) = dfAccumulator / n;
1368 :
1369 798 : return CE_None;
1370 : }
1371 :
1372 : /************************************************************************/
1373 : /* GDALGridDataMetricAverageDistance() */
1374 : /************************************************************************/
1375 :
1376 : /**
1377 : * Average distance between points (data metric).
1378 : *
1379 : * An average distance between the data points found in grid node search
1380 : * ellipse. The distance between each pair of points within ellipse is
1381 : * calculated and average of all distances is set as a grid node value. If
1382 : * there are no points found, the specified NODATA value will be returned.
1383 :
1384 : *
1385 : * \f[
1386 : * Z=\frac{\sum_{i = 1}^{n-1}\sum_{j=i+1}^{n} r_{ij}}{\left(n-1\right)\,n-\frac{n+{\left(n-1\right)}^{2}-1}{2}}
1387 : * \f]
1388 : *
1389 : * where
1390 : * <ul>
1391 : * <li> \f$Z\f$ is a resulting value at the grid node,
1392 : * <li> \f$r_{ij}\f$ is an Euclidean distance between points
1393 : * \f$i\f$ and \f$j\f$,
1394 : * <li> \f$n\f$ is a total number of points in search ellipse.
1395 : * </ul>
1396 : *
1397 : * @param poOptions Algorithm parameters. This should point to
1398 : * GDALGridDataMetricsOptions object.
1399 : * @param nPoints Number of elements in input arrays.
1400 : * @param padfX Input array of X coordinates.
1401 : * @param padfY Input array of Y coordinates.
1402 : * @param padfZ Input array of Z values.
1403 : * @param dfXPoint X coordinate of the point to compute.
1404 : * @param dfYPoint Y coordinate of the point to compute.
1405 : * @param pdfValue Pointer to variable where the computed grid node value
1406 : * will be returned.
1407 : *
1408 : * @return CE_None on success or CE_Failure if something goes wrong.
1409 : */
1410 :
1411 : CPLErr
1412 400 : GDALGridDataMetricAverageDistancePts( const void *poOptions, GUInt32 nPoints,
1413 : const double *padfX, const double *padfY,
1414 : const double *padfZ,
1415 : double dfXPoint, double dfYPoint,
1416 : double *pdfValue,
1417 : void* hExtraParamsIn )
1418 : {
1419 : // TODO: For optimization purposes pre-computed parameters should be moved
1420 : // out of this routine to the calling function.
1421 :
1422 : // Pre-compute search ellipse parameters
1423 : double dfRadius1 =
1424 400 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius1;
1425 : double dfRadius2 =
1426 400 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius2;
1427 : double dfR12;
1428 :
1429 400 : dfRadius1 *= dfRadius1;
1430 400 : dfRadius2 *= dfRadius2;
1431 400 : dfR12 = dfRadius1 * dfRadius2;
1432 :
1433 : // Compute coefficients for coordinate system rotation.
1434 400 : double dfCoeff1 = 0.0, dfCoeff2 = 0.0;
1435 : const double dfAngle =
1436 400 : TO_RADIANS * ((GDALGridDataMetricsOptions *)poOptions)->dfAngle;
1437 400 : const bool bRotated = ( dfAngle == 0.0 ) ? false : true;
1438 400 : if ( bRotated )
1439 : {
1440 400 : dfCoeff1 = cos(dfAngle);
1441 400 : dfCoeff2 = sin(dfAngle);
1442 : }
1443 :
1444 400 : double dfAccumulator = 0.0;
1445 400 : GUInt32 i = 0, n = 0;
1446 :
1447 : // Search for the first point within the search ellipse
1448 147484 : while ( i < nPoints - 1 )
1449 : {
1450 146684 : double dfRX1 = padfX[i] - dfXPoint;
1451 146684 : double dfRY1 = padfY[i] - dfYPoint;
1452 :
1453 146684 : if ( bRotated )
1454 : {
1455 151124 : double dfRXRotated = dfRX1 * dfCoeff1 + dfRY1 * dfCoeff2;
1456 151124 : double dfRYRotated = dfRY1 * dfCoeff1 - dfRX1 * dfCoeff2;
1457 :
1458 151124 : dfRX1 = dfRXRotated;
1459 151124 : dfRY1 = dfRYRotated;
1460 : }
1461 :
1462 : // Is this point located inside the search ellipse?
1463 146684 : if ( dfRadius2 * dfRX1 * dfRX1 + dfRadius1 * dfRY1 * dfRY1 <= dfR12 )
1464 : {
1465 : GUInt32 j;
1466 :
1467 : // Search all the remaining points within the ellipse and compute
1468 : // distances between them and the first point
1469 466634 : for ( j = i + 1; j < nPoints; j++ )
1470 : {
1471 464036 : double dfRX2 = padfX[j] - dfXPoint;
1472 464036 : double dfRY2 = padfY[j] - dfYPoint;
1473 :
1474 464036 : if ( bRotated )
1475 : {
1476 463429 : double dfRXRotated = dfRX2 * dfCoeff1 + dfRY2 * dfCoeff2;
1477 463429 : double dfRYRotated = dfRY2 * dfCoeff1 - dfRX2 * dfCoeff2;
1478 :
1479 463429 : dfRX2 = dfRXRotated;
1480 463429 : dfRY2 = dfRYRotated;
1481 : }
1482 :
1483 464036 : if ( dfRadius2 * dfRX2 * dfRX2 + dfRadius1 * dfRY2 * dfRY2 <= dfR12 )
1484 : {
1485 7264 : const double dfRX = padfX[j] - padfX[i];
1486 7264 : const double dfRY = padfY[j] - padfY[i];
1487 :
1488 7264 : dfAccumulator += sqrt( dfRX * dfRX + dfRY * dfRY );
1489 7264 : n++;
1490 : }
1491 : }
1492 : }
1493 :
1494 146684 : i++;
1495 : }
1496 :
1497 400 : if ( n < ((GDALGridDataMetricsOptions *)poOptions)->nMinPoints
1498 : || n == 0 )
1499 : {
1500 : (*pdfValue) =
1501 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfNoDataValue;
1502 : }
1503 : else
1504 400 : (*pdfValue) = dfAccumulator / n;
1505 :
1506 400 : return CE_None;
1507 : }
1508 :
1509 : /************************************************************************/
1510 : /* GDALGridJob */
1511 : /************************************************************************/
1512 :
1513 : typedef struct _GDALGridJob GDALGridJob;
1514 :
1515 : struct _GDALGridJob
1516 : {
1517 : GUInt32 nYStart;
1518 :
1519 : GByte *pabyData;
1520 : GUInt32 nYStep;
1521 : GUInt32 nXSize;
1522 : GUInt32 nYSize;
1523 : double dfXMin;
1524 : double dfYMin;
1525 : double dfDeltaX;
1526 : double dfDeltaY;
1527 : GUInt32 nPoints;
1528 : const double *padfX;
1529 : const double *padfY;
1530 : const double *padfZ;
1531 : const void *poOptions;
1532 : GDALGridFunction pfnGDALGridMethod;
1533 : GDALGridExtraParameters* psExtraParameters;
1534 : int (*pfnProgress)(GDALGridJob* psJob);
1535 : GDALDataType eType;
1536 :
1537 : void *hThread;
1538 : volatile int *pnCounter;
1539 : volatile int *pbStop;
1540 : void *hCond;
1541 : void *hCondMutex;
1542 :
1543 : GDALProgressFunc pfnRealProgress;
1544 : void *pRealProgressArg;
1545 : };
1546 :
1547 : /************************************************************************/
1548 : /* GDALGridProgressMultiThread() */
1549 : /************************************************************************/
1550 :
1551 : /* Return TRUE if the computation must be interrupted */
1552 618 : static int GDALGridProgressMultiThread(GDALGridJob* psJob)
1553 : {
1554 618 : CPLAcquireMutex(psJob->hCondMutex, 1.0);
1555 621 : (*(psJob->pnCounter)) ++;
1556 621 : CPLCondSignal(psJob->hCond);
1557 621 : int bStop = *(psJob->pbStop);
1558 621 : CPLReleaseMutex(psJob->hCondMutex);
1559 :
1560 621 : return bStop;
1561 : }
1562 :
1563 : /************************************************************************/
1564 : /* GDALGridProgressMonoThread() */
1565 : /************************************************************************/
1566 :
1567 : /* Return TRUE if the computation must be interrupted */
1568 20 : static int GDALGridProgressMonoThread(GDALGridJob* psJob)
1569 : {
1570 20 : int nCounter = ++(*(psJob->pnCounter));
1571 20 : if( !psJob->pfnRealProgress( (nCounter / (double) psJob->nYSize),
1572 : "", psJob->pRealProgressArg ) )
1573 : {
1574 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
1575 0 : *(psJob->pbStop) = TRUE;
1576 0 : return TRUE;
1577 : }
1578 20 : return FALSE;
1579 : }
1580 :
1581 : /************************************************************************/
1582 : /* GDALGridJobProcess() */
1583 : /************************************************************************/
1584 :
1585 103 : static void GDALGridJobProcess(void* user_data)
1586 : {
1587 103 : GDALGridJob* psJob = (GDALGridJob*)user_data;
1588 : GUInt32 nXPoint, nYPoint;
1589 :
1590 103 : const GUInt32 nYStart = psJob->nYStart;
1591 103 : const GUInt32 nYStep = psJob->nYStep;
1592 103 : GByte *pabyData = psJob->pabyData;
1593 :
1594 103 : const GUInt32 nXSize = psJob->nXSize;
1595 103 : const GUInt32 nYSize = psJob->nYSize;
1596 103 : const double dfXMin = psJob->dfXMin;
1597 103 : const double dfYMin = psJob->dfYMin;
1598 103 : const double dfDeltaX = psJob->dfDeltaX;
1599 103 : const double dfDeltaY = psJob->dfDeltaY;
1600 103 : GUInt32 nPoints = psJob->nPoints;
1601 103 : const double* padfX = psJob->padfX;
1602 103 : const double* padfY = psJob->padfY;
1603 103 : const double* padfZ = psJob->padfZ;
1604 103 : const void *poOptions = psJob->poOptions;
1605 103 : GDALGridFunction pfnGDALGridMethod = psJob->pfnGDALGridMethod;
1606 103 : GDALGridExtraParameters *psExtraParameters = psJob->psExtraParameters;
1607 103 : GDALDataType eType = psJob->eType;
1608 103 : int (*pfnProgress)(GDALGridJob* psJob) = psJob->pfnProgress;
1609 :
1610 103 : int nDataTypeSize = GDALGetDataTypeSize(eType) / 8;
1611 101 : int nLineSpace = nXSize * nDataTypeSize;
1612 :
1613 : /* -------------------------------------------------------------------- */
1614 : /* Allocate a buffer of scanline size, fill it with gridded values */
1615 : /* and use GDALCopyWords() to copy values into output data array with */
1616 : /* appropriate data type conversion. */
1617 : /* -------------------------------------------------------------------- */
1618 101 : double *padfValues = (double *)VSIMalloc2( sizeof(double), nXSize );
1619 103 : if( padfValues == NULL )
1620 : {
1621 0 : *(psJob->pbStop) = TRUE;
1622 0 : pfnProgress(psJob); /* to notify the main thread */
1623 0 : return;
1624 : }
1625 :
1626 744 : for ( nYPoint = nYStart; nYPoint < nYSize; nYPoint += nYStep )
1627 : {
1628 641 : const double dfYPoint = dfYMin + ( nYPoint + 0.5 ) * dfDeltaY;
1629 :
1630 25618 : for ( nXPoint = 0; nXPoint < nXSize; nXPoint++ )
1631 : {
1632 24978 : const double dfXPoint = dfXMin + ( nXPoint + 0.5 ) * dfDeltaX;
1633 :
1634 24978 : if ( (*pfnGDALGridMethod)( poOptions, nPoints, padfX, padfY, padfZ,
1635 : dfXPoint, dfYPoint,
1636 : padfValues + nXPoint, psExtraParameters ) != CE_None )
1637 : {
1638 : CPLError( CE_Failure, CPLE_AppDefined,
1639 : "Gridding failed at X position %lu, Y position %lu",
1640 : (long unsigned int)nXPoint,
1641 0 : (long unsigned int)nYPoint );
1642 0 : *(psJob->pbStop) = TRUE;
1643 0 : pfnProgress(psJob); /* to notify the main thread */
1644 0 : break;
1645 : }
1646 : }
1647 :
1648 : GDALCopyWords( padfValues, GDT_Float64, sizeof(double),
1649 : pabyData + nYPoint * nLineSpace, eType, nDataTypeSize,
1650 640 : nXSize );
1651 :
1652 639 : if( *(psJob->pbStop) || pfnProgress(psJob) )
1653 0 : break;
1654 : }
1655 :
1656 103 : CPLFree(padfValues);
1657 : }
1658 :
1659 : /************************************************************************/
1660 : /* GDALGridCreate() */
1661 : /************************************************************************/
1662 :
1663 : /**
1664 : * Create regular grid from the scattered data.
1665 : *
1666 : * This function takes the arrays of X and Y coordinates and corresponding Z
1667 : * values as input and computes regular grid (or call it a raster) from these
1668 : * scattered data. You should supply geometry and extent of the output grid
1669 : * and allocate array sufficient to hold such a grid.
1670 : *
1671 : * Starting with GDAL 1.10, it is possible to set the GDAL_NUM_THREADS
1672 : * configuration option to parallelize the processing. The value to set is
1673 : * the number of worker threads, or ALL_CPUS to use all the cores/CPUs of the
1674 : * computer (default value).
1675 : *
1676 : * Starting with GDAL 1.10, on Intel/AMD i386/x86_64 architectures, some
1677 : * gridding methods will be optimized with SSE instructions (provided GDAL
1678 : * has been compiled with such support, and it is availabable at runtime).
1679 : * Currently, only 'invdist' algorithm with default parameters has an optimized
1680 : * implementation.
1681 : * This can provide substantial speed-up, but sometimes at the expense of
1682 : * reduced floating point precision. This can be disabled by setting the
1683 : * GDAL_USE_SSE configuration option to NO.
1684 : *
1685 : * @param eAlgorithm Gridding method.
1686 : * @param poOptions Options to control choosen gridding method.
1687 : * @param nPoints Number of elements in input arrays.
1688 : * @param padfX Input array of X coordinates.
1689 : * @param padfY Input array of Y coordinates.
1690 : * @param padfZ Input array of Z values.
1691 : * @param dfXMin Lowest X border of output grid.
1692 : * @param dfXMax Highest X border of output grid.
1693 : * @param dfYMin Lowest Y border of output grid.
1694 : * @param dfYMax Highest Y border of output grid.
1695 : * @param nXSize Number of columns in output grid.
1696 : * @param nYSize Number of rows in output grid.
1697 : * @param eType Data type of output array.
1698 : * @param pData Pointer to array where the computed grid will be stored.
1699 : * @param pfnProgress a GDALProgressFunc() compatible callback function for
1700 : * reporting progress or NULL.
1701 : * @param pProgressArg argument to be passed to pfnProgress. May be NULL.
1702 : *
1703 : * @return CE_None on success or CE_Failure if something goes wrong.
1704 : */
1705 :
1706 : CPLErr
1707 27 : GDALGridCreate( GDALGridAlgorithm eAlgorithm, const void *poOptions,
1708 : GUInt32 nPoints,
1709 : const double *padfX, const double *padfY, const double *padfZ,
1710 : double dfXMin, double dfXMax, double dfYMin, double dfYMax,
1711 : GUInt32 nXSize, GUInt32 nYSize, GDALDataType eType, void *pData,
1712 : GDALProgressFunc pfnProgress, void *pProgressArg )
1713 : {
1714 27 : CPLAssert( poOptions );
1715 27 : CPLAssert( padfX );
1716 27 : CPLAssert( padfY );
1717 27 : CPLAssert( padfZ );
1718 27 : CPLAssert( pData );
1719 :
1720 27 : if ( pfnProgress == NULL )
1721 0 : pfnProgress = GDALDummyProgress;
1722 :
1723 27 : if ( nXSize == 0 || nYSize == 0 )
1724 : {
1725 : CPLError( CE_Failure, CPLE_IllegalArg,
1726 0 : "Output raster dimesions should have non-zero size." );
1727 0 : return CE_Failure;
1728 : }
1729 :
1730 : GDALGridFunction pfnGDALGridMethod;
1731 27 : int bCreateQuadTree = FALSE;
1732 :
1733 : /* Potentially unaligned pointers */
1734 27 : void* pabyX = NULL;
1735 27 : void* pabyY = NULL;
1736 27 : void* pabyZ = NULL;
1737 :
1738 : /* Starting address aligned on 16-byte boundary */
1739 27 : float* pafXAligned = NULL;
1740 27 : float* pafYAligned = NULL;
1741 27 : float* pafZAligned = NULL;
1742 :
1743 27 : switch ( eAlgorithm )
1744 : {
1745 : case GGA_InverseDistanceToAPower:
1746 9 : if ( ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->
1747 : dfRadius1 == 0.0 &&
1748 : ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->
1749 : dfRadius2 == 0.0 )
1750 : {
1751 : const double dfPower =
1752 4 : ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->dfPower;
1753 : const double dfSmoothing =
1754 4 : ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->dfSmoothing;
1755 :
1756 4 : pfnGDALGridMethod = GDALGridInverseDistanceToAPowerNoSearch;
1757 4 : if( dfPower == 2.0 && dfSmoothing == 0.0 )
1758 : {
1759 : #ifdef HAVE_SSE_AT_COMPILE_TIME
1760 :
1761 : #define ALIGN16(x) (((char*)(x)) + ((16 - (((size_t)(x)) % 16)) % 16))
1762 :
1763 4 : if( CSLTestBoolean(CPLGetConfigOption("GDAL_USE_SSE", "YES")) &&
1764 : CPLHaveRuntimeSSE() )
1765 : {
1766 3 : pabyX = (float*)VSIMalloc(sizeof(float) * nPoints + 15);
1767 3 : pabyY = (float*)VSIMalloc(sizeof(float) * nPoints + 15);
1768 3 : pabyZ = (float*)VSIMalloc(sizeof(float) * nPoints + 15);
1769 6 : if( pabyX != NULL && pabyY != NULL && pabyZ != NULL)
1770 : {
1771 3 : CPLDebug("GDAL_GRID", "Using SSE optimized version");
1772 3 : pafXAligned = (float*) ALIGN16(pabyX);
1773 3 : pafYAligned = (float*) ALIGN16(pabyY);
1774 3 : pafZAligned = (float*) ALIGN16(pabyZ);
1775 3 : pfnGDALGridMethod = GDALGridInverseDistanceToAPower2NoSmoothingNoSearchSSE;
1776 : GUInt32 i;
1777 1203 : for(i=0;i<nPoints;i++)
1778 : {
1779 1200 : pafXAligned[i] = (float) padfX[i];
1780 1200 : pafYAligned[i] = (float) padfY[i];
1781 1200 : pafZAligned[i] = (float) padfZ[i];
1782 : }
1783 : }
1784 : else
1785 : {
1786 0 : VSIFree(pabyX);
1787 0 : VSIFree(pabyY);
1788 0 : VSIFree(pabyZ);
1789 0 : pabyX = pabyY = pabyZ = NULL;
1790 : }
1791 : }
1792 : #endif // HAVE_SSE_AT_COMPILE_TIME
1793 : }
1794 : }
1795 : else
1796 1 : pfnGDALGridMethod = GDALGridInverseDistanceToAPower;
1797 5 : break;
1798 :
1799 : case GGA_MovingAverage:
1800 4 : pfnGDALGridMethod = GDALGridMovingAverage;
1801 4 : break;
1802 :
1803 : case GGA_NearestNeighbor:
1804 7 : pfnGDALGridMethod = GDALGridNearestNeighbor;
1805 : bCreateQuadTree = (nPoints > 100 &&
1806 : (((GDALGridNearestNeighborOptions *)poOptions)->dfRadius1 ==
1807 7 : ((GDALGridNearestNeighborOptions *)poOptions)->dfRadius2));
1808 7 : break;
1809 :
1810 : case GGA_MetricMinimum:
1811 2 : pfnGDALGridMethod = GDALGridDataMetricMinimum;
1812 2 : break;
1813 :
1814 : case GGA_MetricMaximum:
1815 2 : pfnGDALGridMethod = GDALGridDataMetricMaximum;
1816 2 : break;
1817 :
1818 : case GGA_MetricRange:
1819 2 : pfnGDALGridMethod = GDALGridDataMetricRange;
1820 2 : break;
1821 :
1822 : case GGA_MetricCount:
1823 2 : pfnGDALGridMethod = GDALGridDataMetricCount;
1824 2 : break;
1825 :
1826 : case GGA_MetricAverageDistance:
1827 2 : pfnGDALGridMethod = GDALGridDataMetricAverageDistance;
1828 2 : break;
1829 :
1830 : case GGA_MetricAverageDistancePts:
1831 1 : pfnGDALGridMethod = GDALGridDataMetricAverageDistancePts;
1832 1 : break;
1833 :
1834 : default:
1835 : CPLError( CE_Failure, CPLE_IllegalArg,
1836 0 : "GDAL does not support gridding method %d", eAlgorithm );
1837 0 : return CE_Failure;
1838 : }
1839 :
1840 27 : const double dfDeltaX = ( dfXMax - dfXMin ) / nXSize;
1841 27 : const double dfDeltaY = ( dfYMax - dfYMin ) / nYSize;
1842 :
1843 : /* -------------------------------------------------------------------- */
1844 : /* Create quadtree if requested and possible. */
1845 : /* -------------------------------------------------------------------- */
1846 27 : CPLQuadTree* hQuadTree = NULL;
1847 27 : double dfInitialSearchRadius = 0;
1848 27 : GDALGridPoint* pasGridPoints = NULL;
1849 : GDALGridXYArrays sXYArrays;
1850 27 : sXYArrays.padfX = padfX;
1851 27 : sXYArrays.padfY = padfY;
1852 :
1853 27 : if( bCreateQuadTree )
1854 : {
1855 7 : pasGridPoints = (GDALGridPoint*) VSIMalloc2(nPoints, sizeof(GDALGridPoint));
1856 7 : if( pasGridPoints != NULL )
1857 : {
1858 : CPLRectObj sRect;
1859 : GUInt32 i;
1860 :
1861 : /* Determine point extents */
1862 7 : sRect.minx = padfX[0];
1863 7 : sRect.miny = padfY[0];
1864 7 : sRect.maxx = padfX[0];
1865 7 : sRect.maxy = padfY[0];
1866 17041 : for(i = 1; i < nPoints; i++)
1867 : {
1868 17034 : if( padfX[i] < sRect.minx ) sRect.minx = padfX[i];
1869 17034 : if( padfY[i] < sRect.miny ) sRect.miny = padfY[i];
1870 17034 : if( padfX[i] > sRect.maxx ) sRect.maxx = padfX[i];
1871 17034 : if( padfY[i] > sRect.maxy ) sRect.maxy = padfY[i];
1872 : }
1873 :
1874 : /* Initial value for search radius is the typical dimension of a */
1875 : /* "pixel" of the point array (assuming rather uniform distribution) */
1876 : dfInitialSearchRadius = sqrt((sRect.maxx - sRect.minx) *
1877 7 : (sRect.maxy - sRect.miny) / nPoints);
1878 :
1879 7 : hQuadTree = CPLQuadTreeCreate(&sRect, GDALGridGetPointBounds );
1880 :
1881 17048 : for(i = 0; i < nPoints; i++)
1882 : {
1883 17041 : pasGridPoints[i].psXYArrays = &sXYArrays;
1884 17041 : pasGridPoints[i].i = i;
1885 17041 : CPLQuadTreeInsert(hQuadTree, pasGridPoints + i);
1886 : }
1887 : }
1888 : }
1889 :
1890 :
1891 : GDALGridExtraParameters sExtraParameters;
1892 :
1893 27 : sExtraParameters.hQuadTree = hQuadTree;
1894 27 : sExtraParameters.dfInitialSearchRadius = dfInitialSearchRadius;
1895 27 : sExtraParameters.pafX = pafXAligned;
1896 27 : sExtraParameters.pafY = pafYAligned;
1897 27 : sExtraParameters.pafZ = pafZAligned;
1898 :
1899 27 : const char* pszThreads = CPLGetConfigOption("GDAL_NUM_THREADS", "ALL_CPUS");
1900 : int nThreads;
1901 27 : if (EQUAL(pszThreads, "ALL_CPUS"))
1902 25 : nThreads = CPLGetNumCPUs();
1903 : else
1904 2 : nThreads = atoi(pszThreads);
1905 27 : if (nThreads > 128)
1906 0 : nThreads = 128;
1907 27 : if (nThreads >= (int)nYSize / 2)
1908 0 : nThreads = (int)nYSize / 2;
1909 :
1910 27 : volatile int nCounter = 0;
1911 27 : volatile int bStop = FALSE;
1912 :
1913 : GDALGridJob sJob;
1914 27 : sJob.nYStart = 0;
1915 27 : sJob.pabyData = (GByte*) pData;
1916 27 : sJob.nYStep = 1;
1917 27 : sJob.nXSize = nXSize;
1918 27 : sJob.nYSize = nYSize;
1919 27 : sJob.dfXMin = dfXMin;
1920 27 : sJob.dfYMin = dfYMin;
1921 27 : sJob.dfDeltaX = dfDeltaX;
1922 27 : sJob.dfDeltaY = dfDeltaY;
1923 27 : sJob.nPoints = nPoints;
1924 27 : sJob.padfX = padfX;
1925 27 : sJob.padfY = padfY;
1926 27 : sJob.padfZ = padfZ;
1927 27 : sJob.poOptions = poOptions;
1928 27 : sJob.pfnGDALGridMethod = pfnGDALGridMethod;
1929 27 : sJob.psExtraParameters = &sExtraParameters;
1930 27 : sJob.pfnProgress = NULL;
1931 27 : sJob.eType = eType;
1932 27 : sJob.pfnRealProgress = pfnProgress;
1933 27 : sJob.pRealProgressArg = pProgressArg;
1934 27 : sJob.pnCounter = &nCounter;
1935 27 : sJob.pbStop = &bStop;
1936 27 : sJob.hCond = NULL;
1937 27 : sJob.hCondMutex = NULL;
1938 27 : sJob.hThread = NULL;
1939 :
1940 27 : if( nThreads > 1 )
1941 : {
1942 26 : sJob.hCond = CPLCreateCond();
1943 26 : if( sJob.hCond == NULL )
1944 : {
1945 : CPLError(CE_Warning, CPLE_AppDefined,
1946 0 : "Cannot create condition. Reverting to monothread processing");
1947 0 : nThreads = 1;
1948 : }
1949 : }
1950 :
1951 27 : if( nThreads <= 1 )
1952 : {
1953 1 : sJob.pfnProgress = GDALGridProgressMonoThread;
1954 :
1955 1 : GDALGridJobProcess(&sJob);
1956 : }
1957 : else
1958 : {
1959 26 : GDALGridJob* pasJobs = (GDALGridJob*) CPLMalloc(sizeof(GDALGridJob) * nThreads);
1960 : int i;
1961 :
1962 26 : CPLDebug("GDAL_GRID", "Using %d threads", nThreads);
1963 :
1964 26 : sJob.nYStep = nThreads;
1965 26 : sJob.hCondMutex = CPLCreateMutex(); /* and take implicitely the mutex */
1966 26 : sJob.pfnProgress = GDALGridProgressMultiThread;
1967 :
1968 : /* -------------------------------------------------------------------- */
1969 : /* Start threads. */
1970 : /* -------------------------------------------------------------------- */
1971 128 : for(i = 0; i < nThreads && !bStop; i++)
1972 : {
1973 102 : memcpy(&pasJobs[i], &sJob, sizeof(GDALGridJob));
1974 102 : pasJobs[i].nYStart = i;
1975 102 : pasJobs[i].hThread = CPLCreateJoinableThread( GDALGridJobProcess,
1976 102 : (void*) &pasJobs[i] );
1977 : }
1978 :
1979 : /* -------------------------------------------------------------------- */
1980 : /* Report progress. */
1981 : /* -------------------------------------------------------------------- */
1982 565 : while(nCounter < (int)nYSize && !bStop)
1983 : {
1984 513 : CPLCondWait(sJob.hCond, sJob.hCondMutex);
1985 :
1986 513 : int nLocalCounter = nCounter;
1987 513 : CPLReleaseMutex(sJob.hCondMutex);
1988 :
1989 513 : if( !pfnProgress( nLocalCounter / (double) nYSize, "", pProgressArg ) )
1990 : {
1991 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
1992 0 : bStop = TRUE;
1993 : }
1994 :
1995 513 : CPLAcquireMutex(sJob.hCondMutex, 1.0);
1996 : }
1997 :
1998 : /* Release mutex before joining threads, otherwise they will dead-lock */
1999 : /* forever in GDALGridProgressMultiThread() */
2000 26 : CPLReleaseMutex(sJob.hCondMutex);
2001 :
2002 : /* -------------------------------------------------------------------- */
2003 : /* Wait for all threads to complete and finish. */
2004 : /* -------------------------------------------------------------------- */
2005 128 : for(i=0;i<nThreads;i++)
2006 : {
2007 102 : if( pasJobs[i].hThread )
2008 102 : CPLJoinThread(pasJobs[i].hThread);
2009 : }
2010 :
2011 26 : CPLFree(pasJobs);
2012 26 : CPLDestroyCond(sJob.hCond);
2013 26 : CPLDestroyMutex(sJob.hCondMutex);
2014 : }
2015 :
2016 27 : CPLFree( pasGridPoints );
2017 27 : if( hQuadTree != NULL )
2018 7 : CPLQuadTreeDestroy( hQuadTree );
2019 :
2020 27 : CPLFree(pabyX);
2021 27 : CPLFree(pabyY);
2022 27 : CPLFree(pabyZ);
2023 :
2024 27 : return bStop ? CE_Failure : CE_None;
2025 : }
2026 :
2027 : /************************************************************************/
2028 : /* ParseAlgorithmAndOptions() */
2029 : /* */
2030 : /* Translates mnemonic gridding algorithm names into */
2031 : /* GDALGridAlgorithm code, parse control parameters and assign */
2032 : /* defaults. */
2033 : /************************************************************************/
2034 :
2035 27 : CPLErr ParseAlgorithmAndOptions( const char *pszAlgoritm,
2036 : GDALGridAlgorithm *peAlgorithm,
2037 : void **ppOptions )
2038 : {
2039 27 : CPLAssert( pszAlgoritm );
2040 27 : CPLAssert( peAlgorithm );
2041 27 : CPLAssert( ppOptions );
2042 :
2043 27 : *ppOptions = NULL;
2044 :
2045 27 : char **papszParms = CSLTokenizeString2( pszAlgoritm, ":", FALSE );
2046 :
2047 27 : if ( CSLCount(papszParms) < 1 )
2048 0 : return CE_Failure;
2049 :
2050 27 : if ( EQUAL(papszParms[0], szAlgNameInvDist) )
2051 5 : *peAlgorithm = GGA_InverseDistanceToAPower;
2052 22 : else if ( EQUAL(papszParms[0], szAlgNameAverage) )
2053 4 : *peAlgorithm = GGA_MovingAverage;
2054 18 : else if ( EQUAL(papszParms[0], szAlgNameNearest) )
2055 7 : *peAlgorithm = GGA_NearestNeighbor;
2056 11 : else if ( EQUAL(papszParms[0], szAlgNameMinimum) )
2057 2 : *peAlgorithm = GGA_MetricMinimum;
2058 9 : else if ( EQUAL(papszParms[0], szAlgNameMaximum) )
2059 2 : *peAlgorithm = GGA_MetricMaximum;
2060 7 : else if ( EQUAL(papszParms[0], szAlgNameRange) )
2061 2 : *peAlgorithm = GGA_MetricRange;
2062 5 : else if ( EQUAL(papszParms[0], szAlgNameCount) )
2063 2 : *peAlgorithm = GGA_MetricCount;
2064 3 : else if ( EQUAL(papszParms[0], szAlgNameAverageDistance) )
2065 2 : *peAlgorithm = GGA_MetricAverageDistance;
2066 1 : else if ( EQUAL(papszParms[0], szAlgNameAverageDistancePts) )
2067 1 : *peAlgorithm = GGA_MetricAverageDistancePts;
2068 : else
2069 : {
2070 : fprintf( stderr, "Unsupported gridding method \"%s\".\n",
2071 0 : papszParms[0] );
2072 0 : CSLDestroy( papszParms );
2073 0 : return CE_Failure;
2074 : }
2075 :
2076 : /* -------------------------------------------------------------------- */
2077 : /* Parse algorithm parameters and assign defaults. */
2078 : /* -------------------------------------------------------------------- */
2079 : const char *pszValue;
2080 :
2081 27 : switch ( *peAlgorithm )
2082 : {
2083 : case GGA_InverseDistanceToAPower:
2084 : default:
2085 : *ppOptions =
2086 5 : CPLMalloc( sizeof(GDALGridInverseDistanceToAPowerOptions) );
2087 :
2088 5 : pszValue = CSLFetchNameValue( papszParms, "power" );
2089 : ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
2090 5 : dfPower = (pszValue) ? CPLAtofM(pszValue) : 2.0;
2091 :
2092 5 : pszValue = CSLFetchNameValue( papszParms, "smoothing" );
2093 : ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
2094 5 : dfSmoothing = (pszValue) ? CPLAtofM(pszValue) : 0.0;
2095 :
2096 5 : pszValue = CSLFetchNameValue( papszParms, "radius1" );
2097 : ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
2098 5 : dfRadius1 = (pszValue) ? CPLAtofM(pszValue) : 0.0;
2099 :
2100 5 : pszValue = CSLFetchNameValue( papszParms, "radius2" );
2101 : ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
2102 5 : dfRadius2 = (pszValue) ? CPLAtofM(pszValue) : 0.0;
2103 :
2104 5 : pszValue = CSLFetchNameValue( papszParms, "angle" );
2105 : ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
2106 5 : dfAngle = (pszValue) ? CPLAtofM(pszValue) : 0.0;
2107 :
2108 5 : pszValue = CSLFetchNameValue( papszParms, "max_points" );
2109 : ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
2110 5 : nMaxPoints = (GUInt32) ((pszValue) ? CPLAtofM(pszValue) : 0);
2111 :
2112 5 : pszValue = CSLFetchNameValue( papszParms, "min_points" );
2113 : ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
2114 5 : nMinPoints = (GUInt32) ((pszValue) ? CPLAtofM(pszValue) : 0);
2115 :
2116 5 : pszValue = CSLFetchNameValue( papszParms, "nodata" );
2117 : ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
2118 5 : dfNoDataValue = (pszValue) ? CPLAtofM(pszValue) : 0.0;
2119 5 : break;
2120 :
2121 : case GGA_MovingAverage:
2122 : *ppOptions =
2123 4 : CPLMalloc( sizeof(GDALGridMovingAverageOptions) );
2124 :
2125 4 : pszValue = CSLFetchNameValue( papszParms, "radius1" );
2126 : ((GDALGridMovingAverageOptions *)*ppOptions)->
2127 4 : dfRadius1 = (pszValue) ? CPLAtofM(pszValue) : 0.0;
2128 :
2129 4 : pszValue = CSLFetchNameValue( papszParms, "radius2" );
2130 : ((GDALGridMovingAverageOptions *)*ppOptions)->
2131 4 : dfRadius2 = (pszValue) ? CPLAtofM(pszValue) : 0.0;
2132 :
2133 4 : pszValue = CSLFetchNameValue( papszParms, "angle" );
2134 : ((GDALGridMovingAverageOptions *)*ppOptions)->
2135 4 : dfAngle = (pszValue) ? CPLAtofM(pszValue) : 0.0;
2136 :
2137 4 : pszValue = CSLFetchNameValue( papszParms, "min_points" );
2138 : ((GDALGridMovingAverageOptions *)*ppOptions)->
2139 4 : nMinPoints = (GUInt32) ((pszValue) ? CPLAtofM(pszValue) : 0);
2140 :
2141 4 : pszValue = CSLFetchNameValue( papszParms, "nodata" );
2142 : ((GDALGridMovingAverageOptions *)*ppOptions)->
2143 4 : dfNoDataValue = (pszValue) ? CPLAtofM(pszValue) : 0.0;
2144 4 : break;
2145 :
2146 : case GGA_NearestNeighbor:
2147 : *ppOptions =
2148 7 : CPLMalloc( sizeof(GDALGridNearestNeighborOptions) );
2149 :
2150 7 : pszValue = CSLFetchNameValue( papszParms, "radius1" );
2151 : ((GDALGridNearestNeighborOptions *)*ppOptions)->
2152 7 : dfRadius1 = (pszValue) ? CPLAtofM(pszValue) : 0.0;
2153 :
2154 7 : pszValue = CSLFetchNameValue( papszParms, "radius2" );
2155 : ((GDALGridNearestNeighborOptions *)*ppOptions)->
2156 7 : dfRadius2 = (pszValue) ? CPLAtofM(pszValue) : 0.0;
2157 :
2158 7 : pszValue = CSLFetchNameValue( papszParms, "angle" );
2159 : ((GDALGridNearestNeighborOptions *)*ppOptions)->
2160 7 : dfAngle = (pszValue) ? CPLAtofM(pszValue) : 0.0;
2161 :
2162 7 : pszValue = CSLFetchNameValue( papszParms, "nodata" );
2163 : ((GDALGridNearestNeighborOptions *)*ppOptions)->
2164 7 : dfNoDataValue = (pszValue) ? CPLAtofM(pszValue) : 0.0;
2165 7 : break;
2166 :
2167 : case GGA_MetricMinimum:
2168 : case GGA_MetricMaximum:
2169 : case GGA_MetricRange:
2170 : case GGA_MetricCount:
2171 : case GGA_MetricAverageDistance:
2172 : case GGA_MetricAverageDistancePts:
2173 : *ppOptions =
2174 11 : CPLMalloc( sizeof(GDALGridDataMetricsOptions) );
2175 :
2176 11 : pszValue = CSLFetchNameValue( papszParms, "radius1" );
2177 : ((GDALGridDataMetricsOptions *)*ppOptions)->
2178 11 : dfRadius1 = (pszValue) ? CPLAtofM(pszValue) : 0.0;
2179 :
2180 11 : pszValue = CSLFetchNameValue( papszParms, "radius2" );
2181 : ((GDALGridDataMetricsOptions *)*ppOptions)->
2182 11 : dfRadius2 = (pszValue) ? CPLAtofM(pszValue) : 0.0;
2183 :
2184 11 : pszValue = CSLFetchNameValue( papszParms, "angle" );
2185 : ((GDALGridDataMetricsOptions *)*ppOptions)->
2186 11 : dfAngle = (pszValue) ? CPLAtofM(pszValue) : 0.0;
2187 :
2188 11 : pszValue = CSLFetchNameValue( papszParms, "min_points" );
2189 : ((GDALGridDataMetricsOptions *)*ppOptions)->
2190 11 : nMinPoints = (pszValue) ? atol(pszValue) : 0;
2191 :
2192 11 : pszValue = CSLFetchNameValue( papszParms, "nodata" );
2193 : ((GDALGridDataMetricsOptions *)*ppOptions)->
2194 11 : dfNoDataValue = (pszValue) ? CPLAtofM(pszValue) : 0.0;
2195 : break;
2196 :
2197 : }
2198 :
2199 27 : CSLDestroy( papszParms );
2200 27 : return CE_None;
2201 : }
2202 :
|
