# Get Interpolation transfer relation matrix instead of interpolated values

6 views (last 30 days)
Marcus Becker on 1 Aug 2020
Commented: Joel Lynch on 19 Jan 2021
Hey,
I am working on a simulation where I have scattered measurements which I would like to inter- and extrapolate, ideally to any coordinates but to gridded data would also do. Currently I am working with F = scatteredInterpolant().
The point is, that I am always interpolating to the same coordinates. You can set F.Values = newValues; but the computational time is still unreasonably long in my eyes since it should be a simple matrix multiplication.
So what I am looking for is a function which gives me a transformation matrix A from a set of changing variables at fixed coordinates to another set of fixed coordinates.
So the code should look like this:
% Measurement coordinates
meas_x = data.x;
meas_y = data.y;
% Coordinates to transform to
goal_x = linspace(a,b);
goal_y = linspace(a,b);
A = interpolationTransfMatrix(meas_x,meas_y,goal_x,goal_y,'linear');
for i = 1:n
% Apply the same interploation with different values (reshape would be added)
output = A*data.v(:,i);
end
I feel like that is a too-simple problem to not be solved already, but I haven't found a solution so far, maybe I have the wrong keywords.
Thanks for the help!

Bruno Luong on 1 Aug 2020
Edited: Bruno Luong on 1 Aug 2020
The matrix can be derived from linear interpolation method only, by linear I meant the output is linear to the input. SPLINE method is linear in this sense, PCHIP not.
Here is the code for 2D bilnear (which is linear method).
% Generate source image and 2D grid
z = 10+peaks(50);
z = z(1:30,:);
xgrid = linspace(0,1,size(z,2));
ygrid = linspace(0,1,size(z,1));
% Generate goalx goaly, it can be scattered points, rotate coordinates, morphology etc...
xi = linspace(0,1.1,500);
yi = linspace(0,1.1,500);
[goalx,goaly] = meshgrid(xi, yi);
% Build bilinear interpolation matrix
if isequal(size(goalx),size(goaly))
szout = size(goalx);
else
szout = [];
end
goalx = goalx(:);
goaly = goaly(:);
xgrid = xgrid(:);
ygrid = ygrid(:);
nx = size(xgrid,1);
ny = size(ygrid,1);
dx = diff(xgrid);
dy = diff(ygrid);
if any(dx==0) || any(dy==0)
error('grid not distinct')
end
if ~issorted(xgrid) || ~issorted(ygrid)
error('grid is not sorted')
end
[~,~,ix] = histcounts(goalx,xgrid);
[~,~,iy] = histcounts(goaly,ygrid);
valid = ix > 0 & iy > 0;
ixvalid = ix(valid);
iyvalid = iy(valid);
wx = (goalx(valid)-xgrid(ixvalid)) ./ dx(ixvalid);
wy = (goaly(valid)-ygrid(iyvalid)) ./ dy(iyvalid);
wx = reshape(wx,1,1,[]);
wy = reshape(wy,1,1,[]);
wx = cat(1, 1-wx, wx);
wy = cat(2, 1-wy, wy);
W = wx .* wy;
j = sub2ind([ny nx], iyvalid, ixvalid);
J = reshape(j,1,1,[]);
J = J + [ 0, 1;
ny, ny+1];
nvalid = sum(valid);
i = 1:nvalid;
I = reshape(i,1,1,[]);
I = I + [ 0, 0;
0, 0];
IMAT = sparse([],[],[],length(valid),nx*ny);
IMAT(valid,:) = sparse(I(:),J(:),W(:),nvalid,nx*ny);
% Interpolation using matrix
zi_matmult = IMAT*z(:);
if ~isempty(szout)
zi_matmult = reshape(zi_matmult, szout);
goalx = reshape(goalx, szout);
goaly = reshape(goaly, szout);
end
% and using MATLAB
zi_intepr = interp2(xgrid,ygrid,z,goalx,goaly,'linear',0);
% Check
if ~isvector(zi_matmult)
close all
subplot(1,2,1);
imagesc(zi_intepr)
subplot(1,2,2);
imagesc(zi_matmult)
end
Joel Lynch on 19 Jan 2021
Thanks Bruno for your answser; it helped greatly with a research project! I needed to extend your method to include cubic interpolation as well as various numerical integration techniques. By the time I put evertyhing together, I realized someone else may find it useful too. I put together a toolkit called "mat-op-ex" and realased it on the file exchange.
I'm utterly amazed MATrix-LABoratory doesn't support this approach nativley, so maybe I'm missing something obvious. The only other solution I've found is FUNC2MAT, which even the author admits is extremely expensive to run.
Feel free to provide any feedback.

Marcus Becker on 1 Aug 2020
I have solved the issue by identifying the influence of each point using a value vector which has the value 1 at index i and is otherwise zero and then iterating over all source points. It might not be the most elegant way but this way I can build a function which can use 'nearest','linear' inter-/extrapolation and 'natural' interpolation just like scatteredInterpolant() without buliding an algorithm from scratch:
% Interpolation relation matrix IR
% Measurements
n = 10;
meas_xy = rand(n,2);
% Goal points
m = 100;
goal_xy = rand(m,2);
IR = zeros(m,n);
% First value vector and iteration:
v = zeros(n,1);
v(1) = 1;
F = scatteredInterpolant(meas_xy(:,1),meas_xy(:,2),v,'linear','linear');
IR(:,1) = F(goal_xy);
% Second till n-th iteration
for i = 2:n
v(i-1:i) = [0;1];
F.Values = v;
IR(:,i) = F(goal_xy);
end
% Test and plot
testdata = rand(n,1);
F.Values = testdata;
d = IR*testdata - F(goal_xy);
figure
scatter3(goal_xy(:,1),goal_xy(:,2),IR*testdata)
hold on
scatter3(goal_xy(:,1),goal_xy(:,2),F(goal_xy),10)
hold off