Interpolation and related operations on grids

License: MIT

Language: Julia

Grid operations for the Julia language

Grid Build Status Coverage Status

Continuous objects, such as functions or images, are frequently sampled on a regularly-spaced grid of points. The Grid module provides support for common operations on such objects. Currently, the two main operations are interpolation and restriction/prolongation. Restriction and prolongation are frequently used for solving partial differential equations by multigrid methods, but can also be used simply as fast, antialiased methods for two-fold resolution changes (e.g., in computing thumbnails).

Note: for interpolation, the intended successor to Grid is Interpolations. However, as of August 2015, not all features of Grid are yet implemented in Interpolations (and vice versa, of course).


Within Julia, use the package manager:



To use the Grid module, begin your code with

using Grid


Let's define a quadratic function in one dimension, and evaluate it on an evenly-spaced grid of 5 points:

c = 2.3  # center
a = 8.1  # quadratic coefficient
o = 1.6  # vertical offset
qfunc = x -> a*(x-c).^2 + o
xg = Float64[1:5]
y = qfunc(xg)

From these evaluated points, let's define an interpolation grid:

yi = InterpGrid(y, BCnil, InterpQuadratic)

The last two arguments will be described later. Note that only y is supplied; as with Julia's arrays, the x-coordinates implicitly start at 1 and increase by 1 with each grid point. It's easy to check that yi[i] is equal to y[i], within roundoff error:

julia> y[3]

julia> yi[3]

However, unlike y, yi can also be evaluated at off-grid points:

julia> yi[1.9]

julia> qfunc(1.9)

It's also possible to evaluate the slope of the interpolation function:

julia> v,g = valgrad(yi, 4.25)

julia> 2*a*(4.25-c)

or the second-order derivative (in multidimensional cases, the Hessian matrix):

julia> v,g,h = valgradhess(yi, 4.25)

julia> 2a

Interpolation of a function on a (evenly-spaced) grid that is scaled and/or shifted can be created with CoordInterpGrid that is similar to InterpGrid but takes a range as additional argument:

x = -1.0:0.1:1.0
z = sin(x)

zi = CoordInterpGrid(x, z, BCnil, InterpQuadratic)

julia> z1[-1.0]

While these examples are one-dimensional, you can do interpolation on arrays of any dimensionality. In two or more dimensions CoordInterpGrid expects a tuple of ranges for scaling. For example:

x = -1.0:0.1:1.0
y = 2.0:0.5:10.0
z_2d = Float64[sin(i+j) for i in x, j in y]

z_2di = CoordInterpGrid((x,y), z_2d, BCnil, InterpQuadratic)

julia> z_2di[0.2, 4.1]

The last two parameters of InterpGrid and CoordInterpGrid specify the boundary conditions (what happens near, or beyond, the edges of the grid) and the interpolation order. The choices are specified below:

mode Meaning
BCnil generates an error when points beyond the grid edge are needed
BCnan generate NaN when points beyond the grid edge are needed
BCreflect reflecting boundary conditions
BCperiodic periodic boundary conditions
BCnearest when beyond the edge, use the value at the closest interior point
BCfill produce a specified value when beyond the edge

Most of these modes are activated by passing them as an argument to InterpGrid or CoordInterpGrid as done in the description above. To activate BCfill, pass the number which is to be produced outside the grid as an argument, instead of the mode "BCfill".

The interpolation order can be one of the following:

InterpNearest nearest-neighbor (one-point) interpolation
InterpLinear piecewise linear (two-point) interpolation (bilinear in two dimensions, etc.)
InterpQuadratic quadratic (three-point) interpolation
InterpCubic cubic (four-point) interpolation

Note that quadratic and cubic interpolation are implemented through B-splines which are technically "non-interpolating", meaning that the coefficients of the interpolating polynomial are not the function values at the grid points. InterpGrid solves the tridiagonal system of equations for you, so in simple cases you do not need to worry about such details. InterpQuadratic is the lowest order of interpolation that yields a continuous gradient, and hence is suitable for use in gradient-based optimization, and InterpCubic is similarly the lowest order of interpolation that yields a continuous Hessian.

Note that InterpCubic currently doesn't support all types of boundary conditions; only BCnil and BCnan are supported.

In d dimensions, interpolation references n^d grid points, where n is the number of grid points used in one dimension. InterpQuadratic corresponds to n=3, and InterpCubic corresponds to n=4. Consequently, in higher dimensions quadratic interpolation can be a significant savings relative to cubic spline interpolation.

Low-level interface

It should be noted that, in addition to the high-level InterpGrid interface, Grid also has lower-level interfaces. Users who need to extract values from multi-valued functions (e.g., an RGB image, which has three colors at each position) can achieve significant savings by using this low-level interface. The main cost of interpolation is computing the coefficients, and by using the low-level interface you can do this just once at each x location and use it for each color channel.

Here's one example using the low-level interface, starting from the one-dimensional quadratic example above:

y = qfunc(xg)
# Do the following once
interp_invert!(y, BCnan, InterpQuadratic)   # solve for generalized interp. coefficients
ic = InterpGridCoefs(y, InterpQuadratic)    # prepare for interpolation on this grid

# Do the following for each evaluation point
set_position(ic, BCnil, true, true, [1.8])  # set position to x=1.8, computes the coefs
v = interp(ic, y)                           # extract the value
# Do this if you want the slope at the same point
set_gradient_coordinate(ic, 1)              # change coefs to calc gradient along coord 1
g = interp(ic, y)                           # extract the gradient
# Do this to evaluate the Hessian at that point
set_hessian_coordinate(ic, 1, 2)            # change coefs to calc hessian element H[1,2], i.e. d2/dxdy
h = interp(ic, y)                           # extract hessian element

If this were an RGB image, you could call interp once for the red color channel, once for the green, and once for the blue, with just one call to set_position.

Here is a second example, using a multi-dimensional grid to do multi-value interpolation. Consider that instead of an image, the main data unit of interest is a one-dimensional spectrum, with length Npixels, with three parameters which describe a single spectrum: x1, x2, and x3. If you have a grid of simulations that output spectra for regularly spaced values of these parameters, the following example shows how to interpolate a spectrum with an arbitrary value of x1, x2, and x3.

# Our spectra are stored as pixels in a 4d grid
Npixels = 200
grid = rand(10, 10, 10, Npixels)

# For example, the first spectrum (x1=1, x2=1, x3=1), length Npixels, is stored in 
# the grid as `raw_spec = grid[:,:,:,1]`
grid_size = size(grid) #(10,10,10,200)
grid_strides = strides(grid) #(1,10,100,1000)

strd = stride(grid, 4) #1000

# solve for the generalized interp. coefficients
# but select only the 1st through 3rd axes for interpolation `dimlist = 1:3`
interp_invert!(grid, BCnil, InterpCubic, 1:3)   

# prepare for interpolation on this grid
# but also specify the dimensions and strides of the first three dimensions which we want 
# to interpolate over
ic = InterpGridCoefs(eltype(grid), InterpCubic, grid_size[1:3], grid_strides[1:3])    

# Set the grid position to the indices corresponding to the x1=1.5, x2=1.5, x3=1.5 
# value we wish to interpolate 
set_position(ic, BCnil, false, false, [1.5, 1.5, 1.5])

# Iterate over the pixels in the spectrum to interpolate each pixel into a new array
spec = Array(Float64, (Npixels,))
index = 1
for i = 1:Npixels
    spec[i] = interp(ic, grid, index)
    index += strd

Restriction and prolongation

Suppose you have an RGB image stored in an array img, where the third dimension is of length 3 and specifies the color. You can create a 2-fold smaller version of the image using restrict:

julia> size(img)

julia> imgr = restrict(img, [true,true,false]);

julia> size(imgr)

The second argument to restrict specifies which dimensions should be down-sampled.

Notice that the sizes are not precisely 2-fold smaller; this is because restriction is technically defined as the adjoint of prolongation, and prolongation interpolates (linearly) at intermediate points. For prolongation, you also have to supply the desired size:

julia> img2 = prolong(imgr, [size(img)...]);

julia> size(img2)

If a given dimension has size n, then the prolonged dimension can be either of size 2n-2 (if you want it to be even) or 2n-1 (if you want it to be odd). For odd-sized outputs, the interpolation is at half-grid points; for even-sized outputs, all output values are interpolated, at 1/4 and 3/4 grid points. Having both choices available makes it possible to apply restrict to arrays of any size.

If you plan multiple rounds of restriction, you can get the "schedule" of sizes from the function restrict_size:

julia> pyramid = restrict_size(size(img), [true true true true; true true true true; false false false false])
3x4 Int64 Array:
 961  481  241  121
 541  271  136   69
   3    3    3    3

Restriction and prolongation are extremely fast operations, because the coefficients can be specified in advance. For floating-point data types, this implementation makes heavy use of the outstanding performance of BLAS's axpy.


The main authors of this package are Tim Holy, Tomas Lycken, Simon Byrne, and Ron Rock.

Project Statistics

Sourcerank 7
Repository Size 188 KB
Stars 38
Forks 27
Watchers 8
Open issues 8
Dependencies 0
Contributors 13
Tags 15
Last updated
Last pushed

Top Contributors See all

Tim Holy Tomas Aschan Eric Forgy Simon Byrne Matt Bauman Spencer Lyon Ron Rock Andreas Scheidegger Sebastian Good Scott Thomas Blake Johnson Ian Czekala Andrew Gibb

Packages Referencing this Repo

Interpolation and related operations on grids
Latest release v0.4.1 - Published - 38 stars

Recent Tags See all

v0.4.1 July 29, 2016
v0.4.0 September 24, 2015
v0.3.11 September 16, 2015
v0.3.10 June 03, 2015
v0.3.9 May 20, 2015
v0.3.8 April 13, 2015
v0.3.7 January 05, 2015
v0.3.6 October 18, 2014
v0.3.5 October 09, 2014
v0.3.3 July 04, 2014
v0.3.1 May 08, 2014
v0.3.0 April 17, 2014
v0.2.12 April 02, 2014
v0.2.11 April 02, 2014
v0.2.10 March 01, 2014

Something wrong with this page? Make a suggestion

Last synced: 2016-05-27 14:34:26 UTC

Login to resync this repository