No License

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video

Highlights from
Fast 2-dimensional interpolation

4.5 | 13 ratings Rate this file 19 Downloads (last 30 days) File Size: 6.55 KB File ID: #10772 Version: 1.0
image thumbnail

Fast 2-dimensional interpolation



17 Apr 2006 (Updated )

Provides a 5x-50x speedup over interp2

| Watch this File

File Information

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.


Fast Interpolation inspired this file.

MATLAB release MATLAB 7.1.0 (R14SP3)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (19)
09 Jan 2017 Carolina Raposo

30 Dec 2016 Elias Pirayesh


13 Aug 2016 Hannah

Hannah (view profile)

Slower than interp2 for my matrix ([1024 1280]):
Elapsed time is 0.038501 seconds.
Elapsed time is 0.115177 seconds.

20 Jul 2016 Terence Macquart

Saved me a considerable amount of time.

Comment only
20 Jul 2016 Terence Macquart

07 Sep 2015 song sun

as fast as promised

12 Mar 2015 Pete

Pete (view profile)

Excellent stuff. Though just to note, although there are substantial benefits for small matrices, it appears to be slower than interp2 for large matrices. For example, on my 2012B Windows 7 distribution:

Win-32, matrix=1920x2560
interp2: Elapsed time is 0.436810 seconds. <<
qinterp2: Elapsed time is 0.685080 seconds.

Win-32, matrix=180x180
interp2: Elapsed time is 0.018726 seconds.
qinterp2: Elapsed time is 0.006494 seconds. <<

Win-64, matrix=1920x2560
interp2: Elapsed time is 0.369951 seconds. <<
qinterp2: Elapsed time is 0.655436 seconds.

Win-64, matrix=180x180
interp2: Elapsed time is 0.017490 seconds.
qinterp2: Elapsed time is 0.004462 seconds. <<

23 Feb 2014 ran

ran (view profile)

as fast as promised

02 Jan 2014 Davood Karimi

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?

Comment only
15 Nov 2013 Marc Crapeau

On R2012a, for very small interpolations (2x2 matrices interpolated in 1 point), speed up by a factor 10!

15 Oct 2013 Eric

Eric (view profile)

My testing indicates this is slower than interp2 for Matlab R2013a.

Comment only
07 Mar 2013 Holly Bonstrom

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

10 May 2012 Thomas

Thomas (view profile)

Exactly what i searched for, thanks!

26 Aug 2010 QU RAN

QU RAN (view profile)

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

% Library array error checking (size only)
if size(X)~=size(Y)
error('%s and %s must have the same size',inputname(1),inputname(2));
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));
if ~all(all( abs(diff(diff(Y))) < DIFF_TOL*max(max(abs(Y))) ))
error('%s is not evenly spaced',inputname(2));

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

% 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
flag1 = flagIn & flagCompare;
Zi(flag1) = ...
Z(ind1(flag1)).*(1-dfxi(flag1)) +...
Z(ind2(flag1)).*(dfxi(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)) +...

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 )

% Bilinear interpolation. See
Zi(flagIn) = Z(ind1).*(1-dfxi).*(1-dfyi) + ...
Z(ind2).*dfxi.*(1-dfyi) + ...
Z(ind4).*(1-dfxi).*dfyi + ...

error('Invalid method flag');

end %switch

Comment only
15 Oct 2008 Kyle H

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.

Comment only
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.

Comment only
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.

Contact us