This function performs 2-dimensional interpolation similar to MATLAB's built-in function interp2 with a considerable speed advantage.

qinterp2 may only be used with evenly-spaced, monotonically increasing X and Y library arrays. In addition, there is little error checking on these inputs, as qinterp2's speed boost is derived from avoiding manipulation of these arrays. See the attached image for how interp2 and qinterp2 compare vs. input length. The test was performed on a 2.4GHz Intel Windows XP machine running MATLAB 7 R14SP3.

See file id 10286 (Fast interpolation) for the one-dimensional version of this file.

Thank you for sharing this code. I really need a fast 2D interpolation function for my work; but my tests show that this code is even slightly slower than Matlab's interp2. Am I making a mistake in using this?

I am using optimization and need to interpolate for each iteration. This sped the process up over 60x faster than using interp2 and gives the same results. Thanks!!!

function Zi = qinterp2(X,Y,Z,xi,yi,methodflag)
%QINTERP2 2-dimensional fast interpolation
% qinterp2 provides a speedup over interp2 in the same way that
% qinterp1 provides a speedup over interp1
%
% Usage:
% yi = qinterp2(X,Y,Z,xi,yi) - Same usage as interp2, X and Y should be
% "plaid" (e.g., generated with meshgrid).
% Defaults to bilinear interpolation
% yi = qinterp2(...,flag)
% flag = 0 - Nearest-neighbor
% flag = 1 - Triangular-mesh linear interpolation.
% flag = 2 - Bilinear (equivalent to MATLAB's 'linear')
%
% Usage restrictions
% X(:,n) and Y(m,:) must be monotonically and evenly increasing
% e.g., [X,Y] = meshgrid(-5:5,0:0.025:1);
%
% Examples:
% % Set up the library data for each example
% [X,Y] = meshgrid(-4:0.1:4,-4:0.1:4);
% Z = exp(-X.^2-Y.^2);
%
% % Interpolate a line
% xi = -4:0.03:4; yi = xi;
% Zi = qinterp2(X,Y,Z,xi,yi);
% % Plot the interpolant over the library data
% figure, mesh(X,Y,Z), hold on, plot3(xi,yi,Zi,'-r');
%
% % Interpolate a region
% [xi,yi] = meshgrid(-3:0.3:0,0:0.3:3);
% Zi = qinterp2(X,Y,Z,xi,yi);
% % Plot the interpolant
% figure, mesh(X,Y,Zi);
%
% Error checking
% WARNING: Little error checking is performed on the X or Y arrays. If these
% are not proper monotonic, evenly increasing plaid arrays, this
% function will produce incorrect output without generating an error.
% This is done because error checking of the "library" arrays takes O(mn)
% time (where the arrays are size [m,n]). This function is
% algorithmically independent of the size of the library arrays, and its
% run time is determine solely by the size of xi and yi
%
% Using with non-evenly spaced arrays:
% See qinterp1

% Search array error checking
if size(xi)~=size(yi)
error('%s and %s must be equal size',inputname(4),inputname(5));
end

% Library array error checking (size only)
if size(X)~=size(Y)
error('%s and %s must have the same size',inputname(1),inputname(2));
end
librarySize = size(X);

%{
% Library checking - makes code super slow for large X and Y
DIFF_TOL = 1e-14;
if ~all(all( abs(diff(diff(X'))) < DIFF_TOL*max(max(abs(X))) ))
error('%s is not evenly spaced',inputname(1));
end
if ~all(all( abs(diff(diff(Y))) < DIFF_TOL*max(max(abs(Y))) ))
error('%s is not evenly spaced',inputname(2));
end
%}

% Decide the interpolation method
if nargin>=6
method = methodflag;
else
method = 2; % Default to bilinear
end

% Get X and Y library array spacing
ndx = 1/(X(1,2)-X(1,1)); ndy = 1/(Y(2,1)-Y(1,1));
% Begin mapping xi and yi vectors onto index space by subtracting library
% array minima and scaling to index spacing
xi = (xi - X(1,1))*ndx; yi = (yi - Y(1,1))*ndy;

% Fill Zi with NaNs
Zi = NaN*ones(size(xi));

switch method

% Nearest-neighbor method
case 0
% Find the nearest point in index space
rxi = round(xi)+1; ryi = round(yi)+1;
% Find points that are in X,Y range
flag = rxi>0 & rxi<=librarySize(2) & ~isnan(rxi) &...
ryi>0 & ryi<=librarySize(1) & ~isnan(ryi);
% Map subscripts to indices
ind = ryi + librarySize(1)*(rxi-1);
Zi(flag) = Z(ind(flag));

% Linear method
case 1
% Split the square bounded by (x_i,y_i) & (x_i+1,y_i+1) into two
% triangles. The interpolation is given by finding the function plane
% defined by the three triangle vertices
fxi = floor(xi)+1; fyi = floor(yi)+1; % x_i and y_i
dfxi = xi-fxi+1; dfyi = yi-fyi+1; % Location in unit square

ind1 = fyi + librarySize(1)*(fxi-1); % Indices of ( x_i , y_i )
ind2 = fyi + librarySize(1)*fxi; % Indices of ( x_i+1 , y_i )
ind3 = fyi + 1 + librarySize(1)*fxi; % Indices of ( x_i+1 , y_i+1 )
ind4 = fyi + 1 + librarySize(1)*(fxi-1); % Indices of ( x_i , y_i+1 )

% flagIn determines whether the requested location is inside of the
% library arrays
flagIn = fxi>0 & fxi<librarySize(2) & ~isnan(fxi) &...
fyi>0 & fyi<librarySize(1) & ~isnan(fyi);

% flagCompare determines which triangle the requested location is in
flagCompare = dfxi>=dfyi;

% This is the interpolation value in the x>=y triangle
% Note that the equation
% A. Is linear, and
% B. Returns the correct value at the three boundary points
% Therefore it describes a plane passing through all three points!
%
% From http://osl.iu.edu/~tveldhui/papers/MAScThesis/node33.html
flag1 = flagIn & flagCompare;
Zi(flag1) = ...
Z(ind1(flag1)).*(1-dfxi(flag1)) +...
Z(ind2(flag1)).*(dfxi(flag1)-dfyi(flag1)) +...
Z(ind3(flag1)).*dfyi(flag1);

% And the y>x triangle
flag2 = flagIn & ~flagCompare;
Zi(flag2) = ...
Z(ind1(flag2)).*(1-dfyi(flag2)) +...
Z(ind4(flag2)).*(dfyi(flag2)-dfxi(flag2)) +...
Z(ind3(flag2)).*dfxi(flag2);

case 2 % Bilinear interpolation
% Code is cloned from above for speed
% Transform to unit square
fxi = floor(xi)+1; fyi = floor(yi)+1; % x_i and y_i
dfxi = xi-fxi+1; dfyi = yi-fyi+1; % Location in unit square

% flagIn determines whether the requested location is inside of the
% library arrays
flagIn = fxi>0 & fxi<librarySize(2) & ~isnan(fxi) &...
fyi>0 & fyi<librarySize(1) & ~isnan(fyi);

% Toss all out-of-bounds variables now to save time
fxi = fxi(flagIn); fyi = fyi(flagIn);
dfxi = dfxi(flagIn); dfyi = dfyi(flagIn);

% Find bounding vertices
ind1 = fyi + librarySize(1)*(fxi-1); % Indices of ( x_i , y_i )
ind2 = fyi + librarySize(1)*fxi; % Indices of ( x_i+1 , y_i )
ind3 = fyi + 1 + librarySize(1)*fxi; % Indices of ( x_i+1 , y_i+1 )
ind4 = fyi + 1 + librarySize(1)*(fxi-1); % Indices of ( x_i , y_i+1 )

One should be careful to note that quinterp2 requires X and Y created by meshgrid, NOT ndgrid. This is in contrast with interpn.

In addition, it does a sloppy job dealing with points that are on the borders of X and Y. I was required to scale points that were on the border down by a small factor (1-2*eps), or I would get NaN values.

Finally, compared with interpn, I actually saw a slowdown using this program. It probably works well when xi and yi are large, but in my experience interpn works better when xi and yi are small. Maybe the Mathworks have improved interpn in r2008a.

15 Oct 2008

Nicolas Nicolas

Works well, but for me, MATLAB's built in interp2 with '*METHOD' runs even slightly faster. Maybe the * parameter is rather new.

12 Dec 2006

Ahmed Said

23 Aug 2006

Nathaniel Brahms

Please note that bilinear (and linear) interpolation schemes are now included. All usage and error checking properties are included in the help text.

18 Apr 2006

John D'Errico

While this code should indeed be faster than interp2, there are several factors to consider. It absolutely requires a uniform grid spacing. No error checks are done at all, and certainly not for this property. If your grid is not uniform, this code will produce improper predictions.

A second caveat is that this code does not use the same interpolation method as interp2, despite both being labeled linear. Interp2 uses a tensor product linear interpolant (often known as bilinear interpolation) whereas this code breaks each cell in the grid into a pair of triangles, then interpolates linearly within the corresponding triangle. My own opinion is the triangle interpolant is really the proper way to do "linear" interpolation, and that the tensor product linear interpolant of interp2 is a slight misnomer to call it "linear". Is the bilinear interpolat more or less accurate? Much depends upon the actual surface. The two alternatives should be reasonably close. I have seen surfaces where one or the other was preferable.

What should I rate this code? I picked a number that reflects only the fact that I had to explicitly look through the code itself to make my comments above. If your grid satisfies the requirements of this code and you are accepting of the triangle based linear interpolant, then feel free to read my rating to be a 5.