DICOM slice thickness in Python

Question

I'm working on a project that requires slicing a 3D .stl file into a stack of DICOM slices, and I need to get the distance between the slices. I have figured out how to include the real .stl dimensions (in mm) in each DICOM file, using Autodesk Meshmixer; however, I can't seem to get the correct distance/thickness between slices.

I use ImageJ to read the DICOM stacks, and the voxel depth (distance between slices, if I'm being correct) gets taken the same as the pixel dimensions in the X axis. I want to make the voxel depth to be equal to slice height in the Z axis

As you can see in the last screenshot, the z values are arbitrary.

This is a summary of the code I've done so far:

import numpy as np
import pyvista as pv
import os
import cv2
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from PIL import Image
import pydicom
import datetime
import time

# Load .stl
stl_file = pv.read(stl_path)

# Dimensions of the 3d file
xmin, xmax, ymin, ymax, zmin, zmax = stl_file.bounds

# Number of slices
num_slices = 128
slice_height = (zmax - zmin) / num_slices

# Dimensions of the images in mm 
# (reescaled because the slice gets plotted inside a 396x392 rectangle in a 512x512 image)
# (multiplied by 10 because .stl units are in cm)
pixel_dist_x = (xmax - xmin)/396*10
pixel_dist_y = (ymax - ymin)/392*10

# Date and time
datatime = datetime.datetime.now()

#Loop 
for i in range(num_slices):
    # Here I process each slice as an image (img), and obtain an array
    arr = pil_to_ndarray(img)     
    
    # DICOM conversion
    arr = arr.astype(np.float16)
    ds = pydicom.dataset.FileDataset(os.path.join(slice_folder, f'slice_{i}.dcm'), \
                                     {}, file_meta=pydicom.dataset.FileMetaDataset())
    
    # DICOM Metadata
    ds.PixelData = arr.tobytes()
    ds.SamplesPerPixel = 1
    ds.PhotometricInterpretation = 'MONOCHROME2'
    ds.Rows, ds.Columns = arr.shape[:2]
    ds.BitsAllocated = 16
    ds.BitsStored = 16
    ds.HighBit = 15
    ds.PixelRepresentation = 1
    ds.RescaleIntercept = 1 
    # This slope is set in a way so that Hounsfield units are correct
    ds.RescaleSlope = -1000/15360
    ds.WindowCenter = np.max(arr) / 2
    ds.WindowWidth = np.max(arr)
    
    # Conversion from pixel size to real units
    ds.PixelSpacing = [pixel_dist_x, pixel_dist_y]

    # This is, I believe, the part I should adjust to get the proper
    # spacing between slices
    ds.ImagePositionPatient = [0.0, 0.0, 0.0]
    ds.ImageOrientationPatient = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0]

    # Other data
    ds.PatientName = 'Name'
    ds.PatientID = 'ID'    
    ds.Modality = 'CT'
    ds.InstanceCreationDate = datatime.strftime('%Y%m%d')
    ds.InstanceCreationTime = datatime.strftime('%H%M%S.%f')


    # Save the DICOM file
    ds.save_as(os.path.join(slice_folder, f'slice_{i}.dcm'))

Answer 1

You wrote that:

I need to get the distance between the slices

so I think you need the distance between the DICOM slices and that's easier than it seems. I hope that I have properly understood you question.

First of all, we need to clarify the difference between slice spacing and slice thickness and it seems to me that what you need is the spacing which is defined as

the spacing between adjacent slices, in mm. The spacing is measured from the center-to-center of each slice.

while the slice thickness refers to the (often axial) resolution of the scan. The image below by Materialise should be clear:

Using pydicom you can access the Image Position tag (0020,0032) which

specifies the x, y, and z coordinates of the upper left hand corner of the image; it is the center of the first voxel transmitted.

So the z coordinate allows you to compute the spacing. Given the z of the i slice and the z of the i+1 slice you can get the spacing between this two slices as the difference between the z coordinates. There is no guarantee that the spacing is constant. Moreover not all DICOM tags are always available, some of them are mandatory (Type 1 and Type 2), others are optional (Type 3). The Image Position tag is a Type 1 attribute, i.e. is required and must have a valid value.

The following code allows you to access the z coordinate of the slices:

import os
from pathlib import Path

import pydicom

dicom_dir = r"path\to\dir"


def main() -> None:

    for filename in os.listdir(dicom_dir):
        dicom = pydicom.dcmread(Path(dicom_dir, filename), force=True)                
        z = dicom.ImagePositionPatient[2]


if __name__ == "__main__":
    main()

Other attributes exist, such as Spacing Between Slices (0018, 0088) and Slice Location (0020,1041) but they are optional (Type 3).

Answer 2

Considering what @blunova said, using the 'SliceThickness' attribute also worked for what I need, since (for this case) it's the same to consider thick slices with a null spacing between them, than infinitesimal-thickness slices with non-zero spacing between them.

This way, the thickness is equal to the slice height, and the 'patient' position (initial slice position) is equal to the minimum z position:

import numpy as np
import pyvista as pv
import os
import cv2
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from PIL import Image
import pydicom
import datetime
import time

# Load .stl
stl_file = pv.read(stl_path)

# Dimensions of the 3d file
xmin, xmax, ymin, ymax, zmin, zmax = stl_file.bounds
# Number of slices
num_slices = 128

# Distance between slices
slice_height = (zmax - zmin) / num_slices

#Loop 
for i in range(num_slices):
    # Image and DICOM processing...

    # Slice spacing/thickness
    ds.SliceThickness = slice_height
    ds.ImagePositionPatient = [0.0, 0.0, zmin]