Speeding up linear gridding of irregular points with multiple values (GRIDDATA)

Gridding or interpolating large amounts of data is a common task for some IDL and ENVI users. Here, I am showing a trick that can speed up bi-linear interpolation using a triangulated collection of irregularly gridded points in 2-D. The assumption here is that there are multiple values for each distinct point (x,y), and instead of using GRIDDATA repeatedly for several hundred values at the same locations, the code is pre-computing weights for the triangle corners. This saves computations in the final step and thus achieves a nice speed improvement.

The speed improvement on my computer was going from 31.4 seconds to 8.4 seconds. Here is the output produced by the example code:

 

IDL> grid_speed

% Compiled module: GRID_SPEED.

% Time elapsed: 31.422000 seconds.

% Time elapsed: 8.4380002 seconds.

Mean, Variance, Skewness, Kurtosis

    0.426667    0.0376666      1.26416     0.492909

    0.426667    0.0376666      1.26416     0.492909

 

min, mean, max difference

-1.19209e-007 1.24474e-011 1.19209e-007

 

 

Here is the example code:

 

 

pro grid_speed

 compile_opt idl2,logical_predicate

 

 ;Set up random data points

 ;Let's say 200,000 spatial (X,Y) points with 400 measurements each

 npts = 200000

 nbands = 400

 im = randomu(seed, npts, nbands)

 x = randomu(seed, npts)

 y = randomu(seed, npts)

 

 ;Set up an output gridded space for desired locations

 nx = 768

 ny = 768

 start = [0,0]

 delta = 1d / [nx, ny]

 dim = [nx, ny]

 gridIm1 = fltarr(nx, ny, nbands)

 gridIm2 = fltarr(nx, ny, nbands)

 

 ;traditional approach for bilinear gridding

 tic

 triangulate, x, y, tr

 for i=0, nbands-1 do begin

   gridIm1[0,0,i] = griddata(x, y, im[*,i], triangles=tr, /linear, $

     start=start, delta=delta, dimension=dim)

 endfor

 toc

 

 tic

 triangulate, x, y, tr

 

 ;compute triangle numbers for each input point

 ;multiply by 3 so that triangles are numbered

 ;by the starting index 0, 3, 6, 9, ...

 index = lindgen(n_elements(tr))/3*3

 xt = x[tr[*]]

 yt = y[tr[*]]

 linTr = lindgen(size(tr, /dimensions))

 tr_num = round( $

   griddata(xt, yt, float(index),triangles=linTr, /linear, $

     start=start, delta=delta, dimension=dim))

 

 ;Compute weights for each of the 3 points in the enclosing triangle

  wts = ptrarr(3)

 for i=0, 2 do begin 

   w = griddata(xt, yt, lindgen(n_elements(xt)) mod 3 eq i, $

     triangles=linTr, /linear, $

     start=start, delta=delta, dimension=dim)

   wts[i] = ptr_new(w, /no_copy)

 endfor

 

 ;Compute interpolation for all bands using weights

 ;instead of GRIDDATA

 for i=0, nbands-1 do begin

   gridIm2[0,0,i] = im[tr[tr_num] + i*n_elements(x)] * (*wts[0])

   gridIm2[*,*,i] += im[tr[tr_num+1] + i*n_elements(x)] * (*wts[1])

   gridIm2[*,*,i] += im[tr[tr_num+2] + i*n_elements(x)] * (*wts[2])

 endfor

 toc

 

 ;Verify that the results are the same for both

 ;methods.

 print, 'Mean, Variance,Skewness, Kurtosis'

 print, moment(gridIm1)

 print, moment(gridIm2)

 print

 diff = gridIm2 - gridIm1

 print, 'min, mean, maxdifference'

 print, min(diff, max=maxDiff), mean(diff), maxDiff

end