2D slice from 3D volume (interpolation)

Question

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.

You're using `ndgrid`, but from the documentation for `griddata` it looks like it may expect values with the `meshgrid` convention (see [here](https://stackoverflow.com/questions/19563272/practical-differences-between-meshgrid-and-ndgrid) or [here](https://stackoverflow.com/questions/75465349/how-to-use-ndgrid-with-interp2). Have you accounted for that? Beyond that, it's hard to say much more with a [minimal reproducible example](https://stackoverflow.com/help/minimal-reproducible-example). — horchler, May 30 '23 at 20:48
`griddata` is for when the input data is not sampled on a grid. If you have a 3D image, then I assume it is sampled on a grid. You can then use the much simpler `interp3`. The computation is a lot cheaper. — Cris Luengo, May 31 '23 at 00:48
@CrisLuengo My first thought was interp3, but then when I found it, it seems only available after MATLAB 2023.....whereas mine is MATLAB 2021..... — Nick X Tsui, May 31 '23 at 12:25
I’ve been using MATLAB for a very long time, and `interp3` has always been there. It is a very basic function. — Cris Luengo, May 31 '23 at 12:55
@CrisLuengo Gotcha. Looks like I have it too. But if the plane I am trying to sample is tilted as it is shown above, how can I visualize it as a 2D slice? Because the meshgrid() will generate a 3D volume of coordinates. — Nick X Tsui, May 31 '23 at 13:55

Cris Luengo · Accepted Answer · 2023-05-31T15:50:47.217

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.

2D slice from 3D volume (interpolation)

1 Answers1