I have a CT scan with a shape of (350, 512, 512) and a voxel size of (2, 1.13, 1.13).
I would like to do an interpolation to get a new voxel size of (1,1,1) by using zoom from scipy.
Here is the code I am using
volume = np.stack([dcm.pixel_array for dcm in dicom_files]).astype(np.uint16) #3D volume with original data
volume = volume * slope + intercept #convert to Hounsfield Units
full_ct_spacing = [slice_thickness, pixel_spacing[0], pixel_spacing[1]] #shape of a voxel
new_spacing = np.array([1,1,1]) #target spacing
resize_factor = full_ct_spacing / new_spacing #resize factor
volume = zoom(volume.astype(np.uint16), resize_factor, order=1, mode='reflect') #interpolated volume
I get a new CT scan with a shape of (700, 580, 580) but when I dislay it the grey levels seem to be reversed. Air is white and bones are black.
Is it due to the zoom function ? What is the best way to interpolate CT scan ?
The issue with reversed grey levels in CT scans occurs due to a data type mismatch during interpolation with scipy.ndimage.zoom. Specifically, DICOM intensity values (Hounsfield Units) can be negative (e.g., air: -1000) or positive (e.g., bones). Casting these interpolated values to np.uint16 causes negative values to wrap around (e.g., -1000 becomes 64536), leading to incorrect visualization.
To prevent this, always perform interpolation on floating-point data types (np.float32 or np.float64) to avoid overflow or clipping issues. After processing, convert the data back to an appropriate integer type, such as np.int16, if necessary. This ensures accurate CT scan visualization, where air appears dark and bones bright.
If you need better control or accuracy for medical imaging interpolation, consider libraries specifically tailored for medical imaging, such as SimpleITK or ITK. These libraries handle medical imaging data and spatial resampling more effectively.
Here's an example using SimpleITK
import SimpleITK as sitk
# Convert the 3D volume to a SimpleITK image
sitk_volume = sitk.GetImageFromArray(volume)
sitk_volume.SetSpacing(full_ct_spacing[::-1]) # Original spacing (Z, Y, X)
new_spacing = [1.0, 1.0, 1.0]
resampler = sitk.ResampleImageFilter()
resampler.SetOutputSpacing(new_spacing[::-1])
resampler.SetSize([
int(dim * full_ct_spacing[i] / new_spacing[i])
for i, dim in enumerate(sitk_volume.GetSize())
])
resampler.SetInterpolator(sitk.sitkLinear)
resampled_volume = resampler.Execute(sitk_volume)
# Convert back to NumPy array
volume = sitk.GetArrayFromImage(resampled_volume)