Get Interpolation transfer relation matrix instead of interpolated values
10 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
Bruno Luong
2020-8-1
编辑:Bruno Luong
2020-8-1
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 个评论
Joel Lynch
2021-1-19
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.
更多回答(1 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Map Display 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!