Bilinear interpolation faster than interp2?
Show older comments
I am doing bilinear intepolation and I am looking for a function that is faster than interp2 or griddedinterpolant. I have found the function qinterp2 available here:
Nathaniel Brahms (2023). Fast 2-dimensional interpolation ( https://www.mathworks.com/matlabcentral/fileexchange/10772-fast-2-dimensional-interpolation ), MATLAB Central File Exchange. Retrieved November 4, 2023.
While qinterp2 is quite good, unfortunately it suffers from two limitations:
- The grids X and Y must be equally spaced. I need an arbitrarily spaced grids (still strictly increasing of course)
- It does not allow for extrapolation
I attach a code sample below that explains the issue. I would be really grateful is someone could help me in improving qinterp2 to relax these two limitations (especially limitation 1)
Unimportant note: qinterp2 considers also two other interpolation methods besides bilinear, but I am not interested in those.
clear
clc
close all
n = 7;
nq = 30;
curv = 2.5; %qinterp2 works properly only if curv=1 (i.e. equi-spaced grids)
% Original grids
X = non_linspace(-3,3,n,curv); % column vector
Y = non_linspace(-3,3,n,curv); % column vector
% Finer grids
Xq = non_linspace(-3,3,nq,curv); % column vector
Yq = non_linspace(-3,3,nq,curv); % column vector
% Put data in meshgrid format and plot original data
V = peaks(X',Y);
figure
surf(X,Y,V)
xlabel('X')
ylabel('Y')
title('Original Sampling');
%% Use matlab built-in interp2
% meshgrid format, need X as a row and Y as a column
tic
Vq = interp2(X',Y,V,Xq',Yq,'linear');
toc
% It does not work if the grids are not equally spaced and/or xq and yq go
% beyond the limits
tic
[XX,YY] = meshgrid(X,Y);
[Xqq,Yqq] = meshgrid(Xq,Yq);
Vq2 = qinterp2(XX,YY,V,Xqq,Yqq);
toc
err = max(abs(Vq(:)-Vq2(:)));
if err>1e-10
warning('interp2 and qinterp2 give different results')
end
figure
surf(Xq,Yq,Vq)
xlabel('X')
ylabel('Y')
title('Linear Interpolation Using Finer Grid - interp2');
function xgrid = non_linspace(a,b,n,curv)
xgrid = a+(b-a)*(linspace(0,1,n).^curv)';
end
function Zi = qinterp2(X,Y,Z,xi,yi,methodflag)
% NOTE: Consider only flag=2, ignore other two cases
%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 NOT INTERESTED IN THIS
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 NOT INTERESTED IN THIS
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 THIS IS THE RELEVANT CASE
% 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
% http://en.wikipedia.org/wiki/Bilinear_interpolation
Zi(flagIn) = Z(ind1).*(1-dfxi).*(1-dfyi) + ...
Z(ind2).*dfxi.*(1-dfyi) + ...
Z(ind4).*(1-dfxi).*dfyi + ...
Z(ind3).*dfxi.*dfyi;
otherwise
error('Invalid method flag');
end %switch
end %end function qinterp2
3 Comments
John D'Errico
on 5 Nov 2023
You get what you pay for...
Walter Roberson
on 5 Nov 2023
Are you going to be working with small original grids or large original grids? Small query spaces or large query spaces?
There are techniques that can speed up locating surrounding points on irregular grids, but some of them have a cost in pre-processing the data to make it faster to query the points, and so might become cost-effective for large queries but might be too expensive for small systems.
In particular using a knnsearch on quad-trees can be a time saver if you are querying a lot of points, but does involve set-up time.
It looks like in your example code, each of your grid locations is rectangular, just not equally spaced. If that continues to hold true for your real situation, then there are potentially speed-ups in locating points.
Alessandro Maria Marco
on 5 Nov 2023
Accepted Answer
Bruno Luong
on 5 Nov 2023
Edited: Bruno Luong
on 5 Nov 2023
Try this, about 5 time faster for your example
clc
close all
n = 7;
nq = 30;
curv = 2.5; %qinterp2 works properly only if curv=1 (i.e. equi-spaced grids)
% Original grids
X = non_linspace(-3,3,n,curv); % column vector
Y = non_linspace(-3,3,n,curv); % column vector
% Finer grids
Xq = non_linspace(-3,3,nq,curv); % column vector
Yq = non_linspace(-3,3,nq,curv); % column vector
% Put data in meshgrid format and plot original data
V = peaks(X',Y);
%% Use matlab built-in interp2
% meshgrid format, need X as a row and Y as a column
Vq = interp2(X',Y,V,Xq',Yq,'linear');
%% Bruno implementation
Vq2 = BLU_interp2(X,Y,V,Xq,Yq);
norm(Vq-Vq2,'fro')
timeit(@() interp2(X,Y',V,Xq,Yq','linear'))
timeit(@() BLU_interp2(X,Y,V,Xq,Yq),1)
figure
surf(X,Y,V)
xlabel('X')
ylabel('Y')
title('Original Sampling');
figure
surf(Xq,Yq,Vq2)
xlabel('X')
ylabel('Y')
title('Linear Interpolation Using Finer Grid - BLU\_interp2');
function xgrid = non_linspace(a,b,n,curv)
xgrid = a+(b-a)*(linspace(0,1,n).^curv)';
end
%%
function Vq = BLU_interp2(X,Y,V,Xq,Yq);
nx = length(X);
ny = length(Y);
Xe = X; Xe(1) = -Inf; Xe(end) = Inf;
Ye = Y; Ye(1) = -Inf; Ye(end) = Inf;
xi = discretize(Xq,Xe); % EDIT code
yi = discretize(Yq,Ye); % EDIT code
dx = X(2:end)-X(1:end-1);
dy = Y(2:end)-Y(1:end-1);
wx = (Xq - X(xi)) ./ dx(xi);
wy = (Yq - Y(yi)) ./ dy(yi);
wx = reshape(wx, 1, []);
wy = reshape(wy, [], 1);
wxc = 1-wx;
V00 = V(yi,xi);
V10 = V(yi+1,xi);
V01 = V(yi,xi+1);
V11 = V(yi+1,xi+1);
Vq = (1-wy).*(wxc.*V00 + wx.*V01) + wy.*(wxc.*V10 + wx.*V11);
end
8 Comments
Alessandro Maria Marco
on 5 Nov 2023
Edited: Alessandro Maria Marco
on 5 Nov 2023
Alessandro Maria Marco
on 5 Nov 2023
Note that this problem doesn't happen with
Vq = interp2(X',Y,V,Xq',Yq,'linear');
It this result better than the error message ?
[X,Y] = meshgrid(-3:3);
V = peaks(7);
Xq = -3.5;
Yq = 2;
Vq = interp2(X,Y,V,Xq,Yq,'linear')
Alessandro Maria Marco
on 5 Nov 2023
Edited: Alessandro Maria Marco
on 5 Nov 2023
Thanks for pointing this out! Actually my modification to Bruno's BLU_interp2 allows the function to do linear extrapolation. In other words, if the query point is out of bound, with interp2 you get a NaN value (without error), but BLU_interp2 returns a numeric value. I think that now BLU_interp2 replicates the Matlab griddedInterpolant with options: 'linear', 'linear' for the 2-D case, but please correct me if I am wrong.
I also realize that extrapolation has to be used with care.
Complete code example:
clear
clc
close all
n = 7;
nq = 30;
curv = 2.5;
% Original grids
X = non_linspace(-3,3,n,curv); % column vector
Y = non_linspace(-3,3,n+5,curv); % column vector
% Finer grids
Xq = non_linspace(-3.5,3.5,nq,curv); % column vector
Yq = non_linspace(-3.5,3.5,nq+5,curv); % column vector
% Put data in meshgrid format and plot original data
V = peaks(X',Y);
figure
surf(X,Y,V)
xlabel('X')
ylabel('Y')
title('Original Sampling');
%% Use matlab built-in interp2
% meshgrid format, need X as a row and Y as a column
Vq = interp2(X',Y,V,Xq',Yq,'linear');
%% Bruno implementation
Vq_BLU = BLU_interp2(X,Y,V,Xq,Yq);
err = max(abs(Vq(:)-Vq_BLU(:)));
if err>1e-10
warning('interp2 and BLU_interp2 give different results')
end
a=timeit(@() interp2(X,Y',V,Xq,Yq','linear'))
BLU=timeit(@() BLU_interp2(X,Y,V,Xq,Yq),1)
fprintf('BLU speed-up: %f \n', a/BLU)
figure
surf(Xq,Yq,Vq)
xlabel('X')
ylabel('Y')
title('Linear Interpolation Using Finer Grid - interp2. The Nan values are not plotted');
figure
surf(Xq,Yq,Vq_BLU)
xlabel('X')
ylabel('Y')
title('Linear Interpolation Using Finer Grid - BLU_interp2');
function xgrid = non_linspace(a,b,n,curv)
xgrid = a+(b-a)*(linspace(0,1,n).^curv)';
end
function Vq = BLU_interp2(X,Y,V,Xq,Yq)
nx = length(X);
ny = length(Y);
Xe = X; Xe(1) = -Inf; Xe(end) = Inf;
Ye = Y; Ye(1) = -Inf; Ye(end) = Inf;
xi = discretize(Xq,X);
yi = discretize(Yq,Y);
% Problem if Xq(1)<X(1) or Xq(end)>X(end), then the corresponding xi is NaN
% Likewise for yi
xi(Xq<X(1)) = 1;
xi(Xq>X(nx)) = nx-1;
yi(Yq<Y(1)) = 1;
yi(Yq>Y(ny)) = ny-1;
dx = X(2:end)-X(1:end-1);
dy = Y(2:end)-Y(1:end-1);
wx = (Xq - X(xi)) ./ dx(xi);
wy = (Yq - Y(yi)) ./ dy(yi);
wx = reshape(wx, 1, []);
wy = reshape(wy, [], 1);
wxc = 1-wx;
V00 = V(yi,xi);
V10 = V(yi+1,xi);
V01 = V(yi,xi+1);
V11 = V(yi+1,xi+1);
Vq = (1-wy).*(wxc.*V00 + wx.*V01) + wy.*(wxc.*V10 + wx.*V11);
end
Bruno Luong
on 5 Nov 2023
" BLU_interp2 returns a numeric value. I think that now BLU_interp2 replicates the Matlab griddedInterpolant with options: 'linear', 'linear' for the 2-D case, but please correct me if I am wrong."
I also realize that extrapolation has to be used with care.
You are correct for both statements.
Bruno Luong
on 5 Nov 2023
Edited: Bruno Luong
on 5 Nov 2023
I realize I made a small mistakein ly first code,sorry.Here is a fixed version
function Vq = BLU_interp2(X,Y,V,Xq,Yq)
nx = length(X);
ny = length(Y);
Xe = X; Xe(1) = -Inf; Xe(end) = Inf;
Ye = Y; Ye(1) = -Inf; Ye(end) = Inf;
xi = discretize(Xq,Xe); % Xe and
yi = discretize(Yq,Ye); % Ye here, you don't need the clipping you have added
dx = X(2:end)-X(1:end-1);
dy = Y(2:end)-Y(1:end-1);
wx = (Xq - X(xi)) ./ dx(xi);
wy = (Yq - Y(yi)) ./ dy(yi);
wx = reshape(wx, 1, []);
wy = reshape(wy, [], 1);
wxc = 1-wx;
V00 = V(yi,xi);
V10 = V(yi+1,xi);
V01 = V(yi,xi+1);
V11 = V(yi+1,xi+1);
Vq = (1-wy).*(wxc.*V00 + wx.*V01) + wy.*(wxc.*V10 + wx.*V11);
end
Thanks a lot, it works perfectly!
By the way, I have tested your function vs griddedInterpolant and the speed gain is smaller. I think this is expected, since griddedInterpolant should be faster than interp2.
I attach my test below.
One small question: your function wants the grids in meshgrid format, not ndgrid, correct? So if I create the grids to use griddedInterpolant, then to use your BLU function I have to do a transposition. I personally find the ndgrid format more intuitive (but I realize that in the original question I asked about interp2 that wants the meshgrid, so it's on me!)
clear
clc
close all
nx = 7;
ny = 12;
nqx = 30;
nqy = 35;
curv = 2.5;
% Original grids
X = non_linspace(-5,5,nx,curv); % column vector
Y = non_linspace(-2,2,ny,curv); % column vector
% Finer grids
Xq = non_linspace(-5.5,5.5,nqx,curv); % column vector
Yq = non_linspace(-2.5,2.5,nqy,curv); % column vector
% Put data in ndgrid format and plot original data
V = mypeaks(X,Y');
[Xmat,Ymat] = ndgrid(X,Y);
[Xqmat,Yqmat] = ndgrid(Xq,Yq);
figure
surf(Xmat,Ymat,V)
xlabel('X')
ylabel('Y')
title('Original Sampling');
%% Use matlab built-in griddedInterpolant
% ndgrid format, need X as a column and Y as a row
tic
V_interpolant = griddedInterpolant({X,Y'},V,'linear','linear');
Vq = V_interpolant({Xq,Yq'});
time1=toc
%% Use BLU
tic
Vq_BLU = BLU_interp2_v2(X,Y,V',Xq,Yq); %without ', it does not work!!
Vq_BLU = Vq_BLU';
time2=toc
err = max(abs(Vq(:)-Vq_BLU(:)))
fprintf('Speedup of BLU relative to griddedInt: %f \n',time1/time2)
figure
surf(Xqmat,Yqmat,Vq)
xlabel('X')
ylabel('Y')
title('Linear Interpolation Using Finer Grid - griddedInterpolant');
figure
surf(Xqmat,Yqmat,Vq_BLU)
xlabel('X')
ylabel('Y')
title('Linear Interpolation Using Finer Grid - BLU_interp2');
function xgrid = non_linspace(a,b,n,curv)
xgrid = a+(b-a)*(linspace(0,1,n).^curv)';
end
function F = mypeaks(x,y)
F = 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);
end
function Vq = BLU_interp2_v2(X,Y,V,Xq,Yq)
nx = length(X);
ny = length(Y);
Xe = X; Xe(1) = -Inf; Xe(end) = Inf;
Ye = Y; Ye(1) = -Inf; Ye(end) = Inf;
xi = discretize(Xq,Xe); % Xe and
yi = discretize(Yq,Ye); % Ye here, you don't need the clipping you have added
dx = X(2:end)-X(1:end-1);
dy = Y(2:end)-Y(1:end-1);
wx = (Xq - X(xi)) ./ dx(xi);
wy = (Yq - Y(yi)) ./ dy(yi);
wx = reshape(wx, 1, []);
wy = reshape(wy, [], 1);
wxc = 1-wx;
V00 = V(yi,xi);
V10 = V(yi+1,xi);
V01 = V(yi,xi+1);
V11 = V(yi+1,xi+1);
Vq = (1-wy).*(wxc.*V00 + wx.*V01) + wy.*(wxc.*V10 + wx.*V11);
end
Bruno Luong
on 5 Nov 2023
Edited: Bruno Luong
on 5 Nov 2023
If you want to change to ndgrid (for both V and Vq) I believe you can do this
function Vq = BLU_interp2_t(X,Y,V,Xq,Yq)
nx = length(X);
ny = length(Y);
Xe = X; Xe(1) = -Inf; Xe(end) = Inf;
Ye = Y; Ye(1) = -Inf; Ye(end) = Inf;
xi = discretize(Xq,Xe); % Xe and
yi = discretize(Yq,Ye); % Ye here, you don't need the clipping you have added
dx = X(2:end)-X(1:end-1);
dy = Y(2:end)-Y(1:end-1);
wx = (Xq - X(xi)) ./ dx(xi);
wy = (Yq - Y(yi)) ./ dy(yi);
wx = reshape(wx, [], 1);
wy = reshape(wy, 1, []);
wxc = 1-wx;
V00 = V(xi,yi);
V10 = V(xi,yi+1);
V01 = V(xi+1,yi);
V11 = V(xi+1,yi+1);
Vq = (1-wy).*(wxc.*V00 + wx.*V01) + wy.*(wxc.*V10 + wx.*V11);
end
More Answers (0)
Categories
Find more on Interpolation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
