%% Grid-based Interpolation
% 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.
%
%% Refining a Coarse Grid
% 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')

%
% 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
firstgridvector = F.GridVectors{1}
secondgridvector = F.GridVectors{2}

%
% 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)

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

theinterpolationmethod = F.Method

%
% 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');

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

F(0.9,2.0)

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

generatevalues(0.9,2.0)

%
% 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)

%
% Or using alternative syntax:

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

%
% 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)

%
%% Selecting a Different Interpolation Method
% 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'

%
% 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');

%
%% Interpolating Data in MESHGRID Format
% The 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)

%
% 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)

%
% 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);

%
%% Interpolating Grids in General Dimensions
% 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')

%
% 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)

%
% Evaluation at a single 4D point

F(1.1,2.1,3.1,4.1)

%
% Compare

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

%
% Evaluate at an array of 4D points

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

%
%% Interpolating Grids that have Multiple Values at Each Gridpoint
% 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')

%
% 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

%
%

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

%
%% Interpolating Large Data Sets
% 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')

%
% 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}
secondgridvector = F.GridVectors{2}

%
% 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')

%
% Evaluation would follow the same scaling

F(15,25)

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

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

%
%% Interpolating a Dataset in a Repeated Manner
% 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

%
% 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

%
% Copyright 2011-2014 The MathWorks, Inc.