Artefacts in 3D volume visualization with non-isotropic voxels

Hi everyone,

I am trying to visualize a stack of CT images as a 3D volume. The voxels in the stack are non-isotropic: the pixel size in the x (dxu variable in the code below) and y (dyu variable) directions are equal, but the pixel size in the z direction (dzu variable) is different.

I applied a 3D affine transformation using affinetform3d and then displayed the volume. However, depending on the viewing angle, I notice some artefacts in the visualization. Specifically, certain parts of the volume appear to be duplicated or distorted.

I have attached a figure to illustrate the issue. Below is the code I used:

threshold = 0.2;
Vbone(Vbone < threshold) = 0;
alphaMap = linspace(0, 1, 256).^2;

% Define the affine transformation matrix
A = [dxu 0   0   0; 0   dyu 0   0; 0   0   dzu 0; 0   0   0   1];

% Create the affine transformation object
tform = affinetform3d(A);

% Create the 3D viewer
viewer = viewer3d('BackgroundColor', [0 0 0], 'CameraZoom', 1.3);

% Display the volume with the transformation
vol = volshow(Vbone, ...
    'Parent', viewer, ...
    'Colormap', jet, ...
    'Alphamap', alphaMap, ...
    'RenderingStyle', 'CinematicRendering', ...
    'Transformation', tform);

The figure is attached.

Does anyone have an idea of what might be causing these artefacts? Is there a way to correct the visualization to avoid these distortions?

Thank you in advance for your help!

[Matt J]

What happens if you use imresize3 instead?

[GMabilleau]

Thanks, Matt, for your suggestion to use imresize3 instead of performing an affine transformation. It works well.

Answer 1

[GMabilleau]

The following script has been adapted to use imresize3 instead of an affine transformation, following Matt J’s suggestion:

threshold = 0.2;

Vbone(Vbone < threshold) = 0;

alphaMap = linspace(0, 1, 256).^2;

% Resize the volume to obtain isotropic voxels if not already isotropic

% dxu, dyu, dzu = original voxel sizes (mm/voxel)

voxelSizeIso = min([dxu, dyu, dzu]);

% Compute the scaling factor to apply along each axis

scaleX = dxu / voxelSizeIso;

scaleY = dyu / voxelSizeIso;

scaleZ = dzu / voxelSizeIso;

% 3D resizing with trilinear interpolation

Vbone_iso = imresize3(Vbone, [ ...

round(size(Vbone_visu,1) * scaleY), ...

round(size(Vbone_visu,2) * scaleX), ...

round(size(Vbone_visu,3) * scaleZ)], 'linear');

% Create the 3D viewer and the volshow properties

viewer = viewer3d('BackgroundColor', [0 0 0], 'CameraZoom', 1.3);

vol = volshow(Vbone_iso, ...

'Parent', viewer, ...

'Colormap', jet, ...

'Alphamap', alphaMap, ...

'RenderingStyle', 'CinematicRendering');