This example shows how to create a `griddedInterpolant`

and how to use it effectively to perform grid-based interpolation.

A grid is a common and useful way to organize data. This data format represents values or intensities at discrete grid point locations. Grid-based data also fits naturally within the array-based environment MATLAB® provides.

When we work with gridded data we often need to know the values at locations other than at the grid points. Typically, we may need to refine the grid to improve resolution or de-refine the grid if it contains more detail than we practically need. Grid-based interpolation provides the functionality needed to carry out these tasks.

MATLAB provides the INTERP family of functions to support interpolation on grids that are in `NDGRID`

or `MESHGRID`

format. The ** griddedInterpolant** class provides similar capabilities. It is designed to support the interpolation of grids in NDGRID format and leverage memory and performance advantages where possible. These improvements are readily apparent when interpolating the same grid in a repeated manner. In this scenario the overhead is accumulative. The

`griddedInterpolant`

can generally outperform the INTERP1/2/3/N functions as it is able to cache and reuse the same interpolating function.This example shows you how to create a `griddedInterpolant`

for a gridded data set and then interpolate over a finer grid. We will begin by defining a function that generates values for X and Y input:

generatevalues = @(X,Y)(3*(1-X).^2.*exp(-(X.^2) - (Y+1).^2) ... - 10*(X/5 - X.^3 - Y.^5).*exp(-X.^2-Y.^2) ... - 1/3*exp(-(X+1).^2 - Y.^2));

We can create a 2D grid and then pass it to the `generatevalues`

function to produce values at the grid points. The grid is created from a pair of grid vectors as follows:

xgv = -1.5:0.25:1.5; ygv = -3:0.5:3; [X,Y] = ndgrid(xgv,ygv);

Now generate the value data:

V = generatevalues(X,Y);

We can create an interpolant for this data set that supports interpolation within the grid. Since the interpolant behaves like a function we will give it the variable name F. The 'cubic' option specifies cubic interpolation.

`F = griddedInterpolant(X, Y, V, 'cubic')`

F = griddedInterpolant with properties: GridVectors: {[-1.5000 -1.2500 -1 -0.7500 -0.5000 -0.2500 0 0.2500 0.5000 0.7500 1 1.2500 1.5000] [-3 -2.5000 -2 -1.5000 -1 -0.5000 0 0.5000 1 1.5000 2 2.5000 3]} Values: [13×13 double] Method: 'cubic' ExtrapolationMethod: 'cubic'

The interpolant F has 3 properties: The GridVectors are actually the vectors xgv and ygv we used to create the grid. The interpolant stores the grid in the compact form of GridVectors. This can save memory if the grid is large. GridVectors is a cell array so we can query the contents as follows:

gridvectorprop = F.GridVectors

gridvectorprop = [1×13 double] [1×13 double]

firstgridvector = F.GridVectors{1}

firstgridvector = -1.5000 -1.2500 -1.0000 -0.7500 -0.5000 -0.2500 0 0.2500 0.5000 0.7500 1.0000 1.2500 1.5000

secondgridvector = F.GridVectors{2}

secondgridvector = -3.0000 -2.5000 -2.0000 -1.5000 -1.0000 -0.5000 0 0.5000 1.0000 1.5000 2.0000 2.5000 3.0000

The values at the grid points are stored in the Values array. You can access the values using standard MATLAB syntax to index into the data. For example to inspect a 4-by-5 interval:

first4x5values = F.Values(1:4, 1:5)

first4x5values = 0.0042 0.0028 0.0452 0.3265 0.3007 -0.0050 -0.0671 -0.1285 0.3923 0.9838 -0.0299 -0.2346 -0.5921 0.1483 1.8559 -0.0752 -0.5260 -1.4478 -0.6798 2.4537

The interpolation technique is represented by the Method property. We selected cubic interpolation and our choice is reflected as follows:

theinterpolationmethod = F.Method

theinterpolationmethod = cubic

We can now create a finer grid and use the interpolant to compute the values at these points. We will call these points the query points (Xq, Yq) to distinguish them from our original sample points.

xqgv = -1.5:0.1:1.5; yqgv = -3:0.1:3; [Xq,Yq] = ndgrid(xqgv,yqgv);

We can now evaluate over the refined grid to compute the corresponding values Vq at (Xq, Yq). Since we named our interpolant F, the calling syntax is

Vq = F(Xq, Yq);

We can now generate a plot for comparison with our initial coarse plot.

figure surf(X,Y,V) title('Gridded Data Set', 'fontweight','b');

figure surf(Xq, Yq, Vq); title('Gridded Data Set Refined using Cubic Interpolation', 'fontweight','b');

We can query the interpolant at any location within the domain of the grid.

F(0.9,2.0)

ans = 2.6536

You can compare this interpolated value with the value generated by the analytical expression.

generatevalues(0.9,2.0)

ans = 2.6519

You can query the interpolant using an array of query points as opposed to arrays of query coordinates. We can show this by querying at random locations within the grid.

Xq = -1.5 + 3.*rand(5,2); Vq = F(Xq)

Vq = -1.7032 2.0861 -2.5200 2.1331 6.3799

Or using alternative syntax:

Vq = F(Xq(:,1), Xq(:,2))

Vq = -1.7032 2.0861 -2.5200 2.1331 6.3799

The interpolant supports queries within the domain of the grid. If your query point lies outside the domain of the grid the interpolant will return a NaN.

F(4,5)

ans = 29.0622

You can change the interpolation method on-the-fly. For example, if you wish to use a spline method as opposed to cubic you can change it as follows:

`F.Method = 'spline'`

F = griddedInterpolant with properties: GridVectors: {[-1.5000 -1.2500 -1 -0.7500 -0.5000 -0.2500 0 0.2500 0.5000 0.7500 1 1.2500 1.5000] [-3 -2.5000 -2 -1.5000 -1 -0.5000 0 0.5000 1 1.5000 2 2.5000 3]} Values: [13×13 double] Method: 'spline' ExtrapolationMethod: 'cubic'

We can reevaluate and plot using the spline interpolation method.

xqgv = -1.5:0.1:1.5; yqgv = -3:0.1:3; [Xq,Yq] = ndgrid(xqgv,yqgv); Vq = F(Xq, Yq);

We can now generate a plot for comparison with our initial coarse plot.

figure surf(X,Y,V) title('Gridded Data Set', 'fontweight','b');

figure surf(Xq, Yq, Vq); title('Gridded Data Set Refined using Spline Interpolation', 'fontweight','b');

`MESHGRID`

FormatThe `griddedInterpolant`

class is designed to work with gridded data that conforms to the `NDGRID`

format. This provides support for grids in general N-dimensions, including 1-D which can be regarded as a degenerate grid. In contrast, the `MESHGRID`

format can only support grids in 2D and 3D. Both grid types have identical grid point coordinates; the difference is the format of the coordinate arrays.

If you wish to create a `griddedInterpolant`

using `MESHGRID`

data, you will need to convert the data to `NDGRID`

format. In 2D this involves transposing the arrays as the following example shows.

xgv = -1.5:0.25:1.5; ygv = -3:0.5:3; [X,Y] = meshgrid(xgv,ygv); V = generatevalues(X,Y);

To convert the data to NDGRID format apply a transpose

X = X'; Y = Y'; V = V';

We can now create the interpolant

F = griddedInterpolant(X, Y, V)

F = griddedInterpolant with properties: GridVectors: {[-1.5000 -1.2500 -1 -0.7500 -0.5000 -0.2500 0 0.2500 0.5000 0.7500 1 1.2500 1.5000] [-3 -2.5000 -2 -1.5000 -1 -0.5000 0 0.5000 1 1.5000 2 2.5000 3]} Values: [13×13 double] Method: 'linear' ExtrapolationMethod: 'linear'

Converting 3D `MESHGRID`

data to `NDGRID`

format involves transposing each page of the 3D arrays. This is achieved using the `PERMUTE`

function to interchange the rows (dimension 1) and columns (dimension 2). Here's an example that shows you how:

gv = -3:3; [X,Y,Z] = meshgrid(gv); V = X.^2 + Y.^2 + Z.^2; P = [2 1 3]; X = permute(X,P); Y = permute(Y,P); Z = permute(Z,P); V = permute(V,P);

We can now create the interpolant

F = griddedInterpolant(X, Y, Z, V)

F = griddedInterpolant with properties: GridVectors: {[-3 -2 -1 0 1 2 3] [-3 -2 -1 0 1 2 3] [-3 -2 -1 0 1 2 3]} Values: [7×7×7 double] Method: 'linear' ExtrapolationMethod: 'linear'

Likewise, when querying the interpolant using a `MESHGRID`

, improved performance can be achieved by converting to `NDGRID`

format. For example, if we wish to query our interpolant F using a `MESHGRID`

composed of query points (Xq, Yq, Zq), we could convert the data to `NDGRID`

format as follows:

[Xq, Yq, Zq] = meshgrid(0:0.5:2); Xq = permute(Xq,P); Yq = permute(Yq,P); Zq = permute(Zq,P);

(Xq, Yq, Zq) is now in `NDGRID`

format and can be queried efficiently.

Vq = F(Xq,Yq,Zq);

The `griddedInterpolant`

class is not restricted to 2 and 3 dimensions. You can create an interpolant for 1D, 4D or higher. In practice, the memory required to represent the data may be the limiting factor in higher dimensions. This restriction can impact use in relatively low dimensions, less than ten, depending on the number of grid points and available computing power. The following example illustrates 1D interpolation using the PCHIP interpolation method.

```
X = 1:6;
V = [16 18 21 17 15 12];
F = griddedInterpolant(X,V,'pchip')
```

F = griddedInterpolant with properties: GridVectors: {[1 2 3 4 5 6]} Values: [16 18 21 17 15 12] Method: 'pchip' ExtrapolationMethod: 'pchip'

We can now evaluate the interpolant over a finer interval.

Xq = 1:0.05:6; Vq = F(Xq);

Plotting the query points in blue and the interpolated result in red we get:

plot(X,V,'ob',Xq,Vq,'-r') title('1D Interpolation of a Data Set using the PCHIP Method', 'fontweight','b');

We can create and query a 4D interpolant as follows:

[X1, X2, X3, X4] = ndgrid(1:6); V = X1.^2 + X2.^2 + X3.^2 + X4.^2; F = griddedInterpolant(X1,X2,X3,X4,V)

F = griddedInterpolant with properties: GridVectors: {[1 2 3 4 5 6] [1 2 3 4 5 6] [1 2 3 4 5 6] [1 2 3 4 5 6]} Values: [6×6×6×6 double] Method: 'linear' ExtrapolationMethod: 'linear'

Evaluation at a single 4D point

F(1.1,2.1,3.1,4.1)

ans = 32.4000

Compare

(1.1)^2 + (2.1)^2 + (3.1)^2 + (4.1)^2

ans = 32.0400

Evaluate at an array of 4D points

Xq = 1 + 5*rand(5,4); Vq = F(Xq)

Vq = 48.0991 68.1209 101.5174 87.5036 81.8984

In some applications there may be more than one value associated with each grid point and we may wish to interpolate each value set in turn. For example, if we have a grid representing the pixels in an image we may have three color intensities (RGB) associated with each grid point. There are two ways to interpolate this data. One approach is to create a separate interpolant for each of the three data sets. The other approach is to create a single interpolant and replace the values. The following example illustrates the replacement of values using a single interpolant.

xgv = -1.5:0.25:1.5; ygv = -3:0.5:3; [X,Y] = ndgrid(xgv,ygv); % Create two distinct value sets for this grid V1 = X.^3 - 3*(Y.^2); V2 = 0.5*(X.^2) - 0.5*(Y.^2); % Now create an interpolant for the first value set F = griddedInterpolant(X,Y,V1, 'cubic')

F = griddedInterpolant with properties: GridVectors: {[-1.5000 -1.2500 -1 -0.7500 -0.5000 -0.2500 0 0.2500 0.5000 0.7500 1 1.2500 1.5000] [-3 -2.5000 -2 -1.5000 -1 -0.5000 0 0.5000 1 1.5000 2 2.5000 3]} Values: [13×13 double] Method: 'cubic' ExtrapolationMethod: 'cubic'

We can evaluate the V1 data set on a refined grid and plot the result

xqgv = -1.5:0.1:1.5; yqgv = -3:0.1:3; [Xq,Yq] = ndgrid(xqgv,yqgv); Vq1 = F(Xq,Yq); figure surf(Xq,Yq,Vq1); title('Cubic Interpolation of V1 Dataset', 'fontweight','b');

We can reuse the interpolant to interpolate the second dataset by replacing the Values data.

F.Values = V2

F = griddedInterpolant with properties: GridVectors: {[-1.5000 -1.2500 -1 -0.7500 -0.5000 -0.2500 0 0.2500 0.5000 0.7500 1 1.2500 1.5000] [-3 -2.5000 -2 -1.5000 -1 -0.5000 0 0.5000 1 1.5000 2 2.5000 3]} Values: [13×13 double] Method: 'cubic' ExtrapolationMethod: 'cubic'

Vq2 = F(Xq,Yq); figure surf(Xq,Yq,Vq2); title('Cubic Interpolation of V2 Dataset', 'fontweight','b');

The `griddedInterpolant`

class handles large data sets relatively efficiently. These data sets may consist of a grid of values generated externally and imported into MATLAB. For example, large 2D or 3D images scanned by an external source. In addition, such data sets may not have an explicitly defined grid of coordinate arrays. If the dataset is a large 3D image, the introduction of grid coordinate arrays would quadruple the memory.

The `griddedInterpolant`

class allows you to create an interpolant from the grid of values and a **default grid** is then deduced from the size of the array. This default grid is defined in terms of **grid vectors** - a compact representation of the grid that uses very little memory.

To show this, we can use the PEAKS function to generate an array of values and then create an interpolant for this data set as follows:

```
V = peaks(10);
F = griddedInterpolant(V,'cubic')
```

F = griddedInterpolant with properties: GridVectors: {[1 2 3 4 5 6 7 8 9 10] [1 2 3 4 5 6 7 8 9 10]} Values: [10×10 double] Method: 'cubic' ExtrapolationMethod: 'cubic'

Looking at the **GridVectors** property we can observe the vectors are deduced from the size of the V array. V is a 10x10 array and the corresponding grid vectors are 1:10 and 1:10

firstgridvector = F.GridVectors{1}

firstgridvector = 1 2 3 4 5 6 7 8 9 10

secondgridvector = F.GridVectors{2}

secondgridvector = 1 2 3 4 5 6 7 8 9 10

We can interpolate over a refined grid to improve resolution and fortunately we do not have to create a full grid to do so. We can evaluate using a pair of grid vectors and we package them within curly braces { } to communicate this intent. The default grid has a scaling of 1:10 so we will refine to pick up half intervals using 1:0.5:10. The corresponding query value Vq is as follows:

Vq = F({1:0.5:10, 1:0.5:10});

We can now plot the results side by side.

figure surf(V); title('Sample Values', 'fontweight','b');

figure surf(Vq); title('Cubic Interpolation using a Compact Grid', 'fontweight','b');

**Note:** When we plot the surface the SURF function also uses a default grid to produce the plot. In the first plot the values are represented by a 10-by-10 array and the second a 20-by-20. Hence the 0-to-10 and the 0-to-20 scales on the axes.

The Default Grid can be overridden by specifying grid vectors when you create the interpolant. For example, we could have constructed the interpolant as follows:

`F = griddedInterpolant({10:19, 20:29}, V,'cubic')`

F = griddedInterpolant with properties: GridVectors: {[10 11 12 13 14 15 16 17 18 19] [20 21 22 23 24 25 26 27 28 29]} Values: [10×10 double] Method: 'cubic' ExtrapolationMethod: 'cubic'

Evaluation would follow the same scaling

F(15,25)

ans = -0.0531

The default grid vectors could also have been replaced as follows:

F.GridVectors = {10:19, 20:29}

F = griddedInterpolant with properties: GridVectors: {[10 11 12 13 14 15 16 17 18 19] [20 21 22 23 24 25 26 27 28 29]} Values: [10×10 double] Method: 'cubic' ExtrapolationMethod: 'cubic'

In some applications it may be necessary to interpolate the same dataset in a repeated manner. The `griddedInterpolant`

class can generally handle this scenario more efficiently than the INTERP functions. The `griddedInterpolant`

class is able to reuse data computed during a previous query to speed up the computation of subsequent queries. The following example shows this advantage:

Sample data set

[X, Y, Z] = ndgrid(1:100); V = X.^2 + Y.^2 + Z.^2;

Performance data for `INTERPN`

tic; for i = 1:1000 Xq = 100*rand(); Yq = 100*rand(); Zq = 100*rand(); Vq = interpn(X,Y,Z,V,Xq,Yq,Zq,'cubic'); end interpnTiming = toc

interpnTiming = 22.5707

Performance data for `griddedInterpolant`

tic; F = griddedInterpolant(X,Y,Z,V, 'cubic'); for i = 1:1000 Xq = 100*rand(); Yq = 100*rand(); Zq = 100*rand(); Vq = F(Xq,Yq,Zq); end griddedInterpolantTiming = toc

griddedInterpolantTiming = 2.5020

Was this topic helpful?