Grid operations for the Julia language

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, Grid.jl is deprecated in favor of Interpolations. Julia 0.5 is the latest release on which Grid.jl will be installable.

Installation

Within Julia, use the package manager:

Pkg.add("Grid")

Usage

To use the Grid module, begin your code with

using Grid

Interpolation

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]
5.569000000000003

julia> yi[3]
5.569000000000003

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

julia> yi[1.9]
2.8959999999999995

julia> qfunc(1.9)
2.8959999999999995

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

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

julia> 2*a*(4.25-c)
31.59

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

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

julia> 2a
16.2

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> zi[-1.0]
-0.8414709848078965

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]
-0.9156696045493824

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
end

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)
(1920,1080,3)

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

julia> size(imgr)
(961,541,3)

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)
(1920,1080,3)

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.

Credits

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