Notes on Interpolation in IDL

IDL provides routines to do 1-d and 2-d interpolation of data. This page discusses some useful techniques.

It is assumed that the JHU/APL/S1R library is on your IDL path.

1-D Interpolation

1-d interpolation will be illustrated by showing an example.

Using the user's library routine interpol

First make some 1-d data to interpolate:

    x = maken(10,20,10)    ; 10 x values from 10 to 20.
    y = randomu(k,10)      ; 10 random y values.
    plot,x,y,psym=-2       ; Take a look at the data.

Next make an array of desired interpolation locations:

    xx = maken(10,20,40)   ; Make 40 points.
    yy = interpol(y,x,xx)  ; Interpolate using interpol.
    oplot,xx,yy,psym=4     ; Look at result.

Note that, as shown in this case, interpolated values can miss extreme values. If interpolation is done to make a non-uniformly spaced data set uniform then enough new points will be needed to minimize this effect.

interpol is a user's library routine. The calling syntax is:
new_y_array = interpol(old_y_array, old_x_array, new_x_array)

Using the built-in routine interpolate

One problem with the routine interpol is that it is not very fast. This becomes important for repeated interpolations. IDL has a relatively recent routine called interpolate which is much faster.

interpolate needs as input the array indices at which to interpolate. These indices will in general be fractional. Often the indices are not directly known, what is known is the x coordinates of the desired interpolated points. These may be converted to actual indices and used in interpolateas follows.

   yi = interpol(indgen(n_elements(y)),x,xx) ; Use interpol to find indices.
   yy = interpolate(y,yi)                    ; Do interpolation.

Note that the advantage occurs for repeated interpolation at the same indices.

SUMMARY
Let old_x, old_y be a function to be interpolated. Perhaps old_x is non-uniform and you need to make it uniform. new_x is the desired x coordinates (perhaps uniform). x_index = interpol( indgen(n_elements(old_x)), old_x, new_x) new_y = interpolate( old_y, x_index) For interpolate on 2-D arrays use the keyword /GRID when possible. This applies if all the lines have the same interpolation and all the columns have the same interpolation. Need both x_index and y_index for 2-D arrays. new_img = interpolate(old_img, x_index, y_index, /grid)

   Benchmark comparing the routines interpol and interpolate:

   On our HP 7/35 to interpolate 100,000 points:
   interpol took 8.8 seconds.
   interpolate took 0.064 seconds, over 130 times faster.

Interpolating to uniformly space data points

Sometimes data is sampled non-uniformly along a curve. An example can be constructed as follows:

    x = makenlog(10,500,600)
    y = sin(x/10)

Plotting just y gives:

Plotting x and y gives the correct view:

The data can be interpolated to a uniform spacing as follows:

    xx = maken(10,500,1000)                    ; Desired new X positions.
    yi = interpol(indgen(n_elements(y)),x,xx)  ; Find their indices.
    yy = interpolate(y,yi)                     ; Interpolate corresponding Ys.

Plotting just the new y, yy, now gives:

This has importance in correcting images for non-uniformly spaced rows or columns.

2-D Interpolation

Here is a 200 x 200 pixel image, img, based on the non-uniform samples in the last example:

The rows are evenly sampled, the columns are not. The non-uniformly spaced image columns may be remapped as follows. Assume the actual data coordinates of each are given by:

    x = maken(10,20,200)
    y = makenlog(10,500,200)

Let nx and ny be the X and Y size of the interpolated image. In this example nx=200 and ny=200, this keeps the image size the same. Uniformly spaced coordinates may be found as:

    ny = 200
    yu = maken(10,500,ny)

The array indices in the Y dimension which correspond to these uniformly spaced y coordinates are given by:

    yi = interpol(indgen(n_elements(y)),y,yu)

The X dimension is already uniformly sampled so it's indices are given by:

    nx = 200
    xi = indgen(nx)

Here again the number could be changed (like xi = maken(0,199,400)).

Now 2 index arrays are set up, one for X and one for Y:

    xx = rebin(xi,nx,ny)
    yy = rebin(transpose(yi),nx,ny)

Interpolation is very simple:

    out = interpolate(img,xx,yy)

The result is:

The same indices, xx and yy, could be used to interpolate many images.

2-D interpolatin using SFITS

The IDL routine SFITS can do a 2-D polynomial fit to random data points (also regularly gridded points). Most interesting, SFITS can return the coeffecients of the 2-D polynomial fit and they can be used to compute the value for any (x,y) in the valid range. Here is an example:

Generate some random data:

	x = randomu(k,100)-.5
	y = randomu(k,100)-.5
	z = sin(x*5)+cos(y*5)

Pack it up for SFITS:

	d = [transpose(x),transpose(y),transpose(z)]

Do fit. The /irregular keywords means the data is given as random points.

	deg = 5
	r = sfit(d,deg,/irregular,kx=kx)

Don't care about the returned result which is just the fit at the given points. Really want to polynomial coeffecients returned in kx.

Now generate and plot a surface based on the returned coeffecients:

	makenxy,-.5,.5,100,-.5,.5,100,xx,yy			; 2-d X and Y arrays.
	t = fltarr(100,100)					; Surface array.
	for i=0,deg do for j=0,deg do t=t+kx(j,i)*xx^i*yy^j		; Compute surface.
	surface,t,maken(-.5,.5,100),maken(-.5,.5,100),/save	; Plot surface.
	plots,/t3d,x,y,z,psym=2,col=255				; Original points.
	tvscl,t							; Display fitted array.

The fit breaks down when the degree gets very high even if there are enough input points. Look at the results before assuming the fit is good. For 300 points the above fit was pretty good up to maybe degree 12.