Get Interpolation transfer relation matrix instead of interpolated values

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!

Answer 1

[Bruno Luong]

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

[Marcus Becker]

If I read and tested it correctly, the source data needs to be gridded, not regularly, but gridded and not scattered, right?

Edit: In my case I have scattered measurements / source data points and want to get a regular grid. The nice-to-have addition would be if it also works for a scattered goal point cloud.

[Bruno Luong]

Yes.

You have to dig in and do reverse engineering to know how scatteredInterpolant works (triangulation, etc...).

It's perhaps not difficult for the first 2 methods:

• 'nearest': (the matrix can be build from an interpolation M sample-points with data is 1, 2, 3, ..., M
• 'linear': Use cartesianToBarycentric of delaunayTriangulation to build your matrix
• 'natural': I have no idea since I don't know what behind the method.

[Bruno Luong]

Here we go scatteredInterpolant matrix:

x = rand(1000,1);
y = rand(1000,1);

queryx = rand(10,1);
queryy = rand(10,1);

T = delaunayTriangulation(x,y);
P = [queryx(:) queryy(:)];
m = size(P,1);
n = size(T.Points,1);
ID = pointLocation(T,P);
W = cartesianToBarycentric(T, ID, P);
I = repmat(1:m,1,3);
IMAT = sparse(I,T.ConnectivityList(ID,:),W,m,n);

z = rand(size(x));
IMAT*z(:)

F = scatteredInterpolant(x,y,z);
F(queryx,queryy)

[Marcus Becker]

works, thanks for the effort!

[Marcus Becker]

As feedback, switching to this method of preparing the matrix and using it instead of the interpolation decreased the computational time to 0.7% (regarding this step) of what it previously was!

[Joel Lynch]

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.

Answer 2

[Marcus Becker]

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