2D slice from 3D volume (interpolation)

I have a 3D volume (black) with each voxel having a known 3D coordinates (x,y,z), and a corresponding intensity value f(x,y,z).

Now I like to get a slice (blue) from the 3D volume. All (xq, yq, zq) coordinates of each voxel of the slice is know, and I'd like to get the value f(xq, yq, zq) for each voxel on the slice.

The psuedo code I tried is something like this

[xxq,yyq,zzq]= ndgrid(xq, yq, zq);   % coordinates of sampling points
[xx,yy,zz] = ndgrid(xx,yy,zz);       % coordinates of data points
fq = griddata(x,y,z,f,xq,yq,zq);     % 3D interpolation 

but the results look completely ridiculous, and I don't know why.

Any pointer of how to do it is very appreciated.

Because your image data is sampled on a regular grid, you should use interp3 for interpolation. griddata is for the case where the input data is sampled at arbitrary locations, it does a Delaunay triangulation of the data coordinates to be able to interpolate. This is obviously much more expensive than needed on the simple case of the data sampled on a regular grid.

Here's an example for how to apply interp3 to a 3D image to sample in a new 2D grid along some plane within the 3D volume:

% Generate example input image
x = -50:50;
y = permute(x,[2,1]);
z = permute(x,[1,3,2]);
img = x.^2 + y.^2 + z.^2;

% We need the input grid coordinates as full 3D arrays, not sure why...
[y,x,z] = ndgrid(y,x,z);

% Create 2D arrays with output coordinates on a plane
% TODO: generate the xq, yq, zq arrays containing the 3D coordinates
%       for each output pixel
[xq,yq] = meshgrid(-50:50, -50:50);
zq = zeros(size(xq)) + xq/10 + 6;

% Interpolate
out = interp3(x,y,z,img,xq,yq,zq);

If the image can be assumed to have coordinates 1:N in each dimension, then we don't need to generate the x,y,z grids:

[xq,yq] = meshgrid(1:100, 1:100);
zq = zeros(size(xq)) + xq/10 + 56;
out = interp3(img,xq,yq,zq);

Please note the order of the x and y arguments to the various functions. MATLAB makes a mess of dimension orders. The first array dimension is vertical, which makes sense for matrices when doing linear algebra, but doesn't match the expected direction of the x axis in other fields (i.e. the first dimension is y, the second is x). Some functions take inputs in the array axis order (row, column, ...), and some functions take the inputs in the x, y, z order. The documentation is always clear about this by naming arguments "row" and "column", or "x1", "x2" and "x3", for the first case, and naming them "x", "y" and "z" for the second case.

meshgrid and interp3 in the x, y, z order, but ndgrid takes arguments in the array axis order, and hence the ugly syntax [y,x,z] = ndgrid(y,x,z) in the code above.