Multiple griddata calls into a single one (same grid)

Hi,

I'm using griddata to regrid points (xpoint, ypoint, pointval) into a regular grid (xgrid, ygrid). I used to do this in 2D as follows:

x = linspace(-1,1,200);
y = linspace(200,600,200);
[xgrid,ygrid] = meshgrid(x,y);
gridded = griddata(xpoint,ypoint,pointval,xgrid,ygrid,'v4'); 

But now I have expanded my analysis and I have different sets of points. All have the same (xpoint,ypoint) coordinates but different pointval. The straight way of doing this would be:

gridded = arrayfun(@(nn)griddata(xpoint,ypoint,pointval(:,nn),xgrid,ygrid,'v4'),1:npoints,'uniformoutput',false)

However this is effectively similar to a for loop and takes some time. Do you think there is a smarter way of calling griddata only once?

Thanks!

Answer 1

[Star Strider]

The arrayfun function is significantly slower than an explicit loop, at least in my experience.

I would just do something like this:

for nn = 1:npoints
    gridded{k} = griddata(xpoint,ypoint,pointval(:,nn),xgrid,ygrid,'v4'); 
end

I am sure there are applications where arrayfun is useful (perhaps with small arrays in anonymous functions), however llikely not here.

Also, the original arguments were (xpoint,xpoint). I changed them here to (xpoint,ypoint) since that also appears, so check that to be certain it is correct.

[Albert]

Thank you, I changed the original question it was (xpoint,ypoint) indeed. In my case I checked arrayfun and a for loop and they both last the same...

[Star Strider]

Perhaps arrayfun has been optimised since I last tried it.

I will delete my answer in a few hours.

[Albert]

No problem, maybe it depends on the type of operation inside

[Bruno Luong]

IMO arrayfun is always slower than for-loop. It is just more visible when the unitary loop operation is fast. Here the interpolation is not fast so using arrayfun or for-loop doesn't matter speedwise.

Answer 2

[Bruno Luong]

If you are ready to trade 'v4' method for something else, you can use scatteredInterpolant

% example of fake data
x = -3 + 6*rand(50,1);
y = -3 + 6*rand(50,1);
v = sin(x).^4 .* cos(y);
pointval = v + (0:10);
% fake grid
[xgrid,ygrid] = meshgrid(-3:0.1:3);

F = scatteredInterpolant(x,y,pointval(:,1));
F.Method = 'natural';
nv = size(pointval,2);
gridded = zeros([size(xgrid),nv]);
for k=1:size(pointval,2)
    F.Values = pointval(:,k);
    gridded(:,:,k)  = F(xgrid,ygrid);
end 

[Albert]

I always had the impression that v4 performed best but takes more time. I can also use the former griddata function with the 'cubic' method and that's superfast too. Thanks!

Answer 3

[Bruno Luong]

For nearest/linear/cubic method you can build the matrix https://www.mathworks.com/matlabcentral/fileexchange/85939-mat-op-ex, followed up from this thread

The output is simply matrix x vector of values, so you can multiply by many input vectors in one shot.