Fast 2-dimensional interpolation

Provides a 5x-50x speedup over interp2

Description

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.

Acknowledgements

The author wishes to acknowledge the following in the creation of this submission:
Fast interpolation

• MATLAB release: MATLAB 7.1.0 (R14SP3)

[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.

[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.

[Ahmed Said]

 

[Nicolas Nicolas]

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

[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.

[QU RAN]

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