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 ?

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 ?

• python
• numpy
• scipy

• Please provide a minimal reproducible example – Julien Commented yesterday
• Hard to tell without knowing how you extract that array, and how you create a new image from volume. But it seems that there are some dicom involved. dicom image have an "interpretation" parameter. Among which some for which 0 is black, and some for which 0 is white (a bit like tiff, with minisblack and miniswhite). So, if you extract an array from an image for which minimum level is white, and then recreate from this array (or its zoomed version) another image for which min level is black, you can expect that. But we wouldn't know without a minimal reproducible example. – chrslg Commented yesterday

1 Answer 1

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)