-
-
Notifications
You must be signed in to change notification settings - Fork 14
/
qinterp2.m
171 lines (150 loc) · 6.39 KB
/
qinterp2.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
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