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

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]