Bilinear interpolation faster than interp2?

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:

1. The grids X and Y must be equally spaced. I need an arbitrarily spaced grids (still strictly increasing of course)
2. 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

Elapsed time is 0.025106 seconds.

% 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

Elapsed time is 0.011711 seconds.

err = max(abs(Vq(:)-Vq2(:)));
if err>1e-10
    warning('interp2 and qinterp2 give different results')
end

Warning: interp2 and qinterp2 give different results

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

[John D'Errico]

You get what you pay for...

[Walter Roberson]

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]

I use rectangular grids, just not equally spaced.

Answer 1

[Bruno Luong]

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

ans = 0

timeit(@() interp2(X,Y',V,Xq,Yq','linear'))

ans = 1.6877e-04

timeit(@() BLU_interp2(X,Y,V,Xq,Yq),1)

ans = 2.6770e-05

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

[Alessandro Maria Marco]

This is great, significantly faster, thanks!

Unfortunately it does not handle the case where Xq(1)<X(1) and/or Xq(end)>X(end), and likewise for Y and Yq. I tried setting

% Finer grids
Xq = non_linspace(-3.5,3.5,nq,curv); % column vector
Yq = non_linspace(-3.5,3.5,nq,curv); % column vector

but it returns with the following error message:

Array indices must be positive integers or logical values.

Error in BLU_interp2 (line 13)
wx = (Xq - X(xi)) ./ dx(xi);

Error in test_qinterp2 (line 30)
Vq2 = BLU_interp2(X,Y,V,Xq,Yq);

The reason is that if Xq is out of bounds, i.e. if any of the point in Xq lies outside of the range [X(1),X(end)], then

discretize(Xq,X)

returns a NaN value

Note that this problem doesn't happen with

Vq = interp2(X',Y,V,Xq',Yq,'linear');

[Alessandro Maria Marco]

I think I managed to fix the problem. I copy and paste here your function with the modification I have added

Please note that the following lines of code that I have added:

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

And here is the complete function:

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

[Torsten]

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

Vq = NaN

[Alessandro Maria Marco]

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

a = 1.6875e-04

BLU=timeit(@() BLU_interp2(X,Y,V,Xq,Yq),1)

BLU = 3.1745e-05

fprintf('BLU speed-up: %f \n', a/BLU)

BLU speed-up: 5.315572 

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]

" 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]

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

[Alessandro Maria Marco]

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

time1 = 0.0026

%% Use BLU
tic
Vq_BLU = BLU_interp2_v2(X,Y,V',Xq,Yq); %without ', it does not work!!
Vq_BLU = Vq_BLU';
time2=toc

time2 = 0.0222

err = max(abs(Vq(:)-Vq_BLU(:)))

err = 0

fprintf('Speedup of BLU relative to griddedInt: %f \n',time1/time2)

Speedup of BLU relative to griddedInt: 0.118212 

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]

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