6 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!

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.

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!