Get Interpolation transfer relation matrix instead of interpolated values
8 views (last 30 days)
Show older comments
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!
0 Comments
Accepted Answer
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
6 Comments
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.
More Answers (1)
See Also
Categories
Find more on Map Display in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!