LCOV - code coverage report
Current view: directory - alg - gdalgrid.cpp (source / functions) Found Hit Coverage
Test: gdal_filtered.info Lines: 747 623 83.4 %
Date: 2012-12-26 Functions: 18 18 100.0 %

       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                 : 

Generated by: LCOV version 1.7