AAPM Thoracic Auto-segmentation Challenge Forum

Go back to competition Back to thread list Post in this thread

> Python code to read contour from RTSTRUCT

This is a code snippet I used to read contour information and optionally convert to label maps. In case anyone is interested.
pydicom, skimage, numpy, matplotlib are required. Modify the training data path as needed.

import os, sys, glob
import numpy as np
import dicom
import matplotlib.pyplot as plt
%matplotlib inline
from skimage.draw import polygon

def read_structure(structure):
contours = []
for i in range(len(structure.ROIContourSequence)):
contour = {}
contour['color'] = structure.ROIContourSequence[i].ROIDisplayColor
contour['number'] = structure.ROIContourSequence[i].RefdROINumber
contour['name'] = structure.StructureSetROISequence[i].ROIName
assert contour['number'] == structure.StructureSetROISequence[i].ROINumber
contour['contours'] = [s.ContourData for s in structure.ROIContourSequence[i].ContourSequence]
return contours

def get_mask(contours, slices):
z = [s.ImagePositionPatient[2] for s in slices]
pos_r = slices[0].ImagePositionPatient[1]
spacing_r = slices[0].PixelSpacing[1]
pos_c = slices[0].ImagePositionPatient[0]
spacing_c = slices[0].PixelSpacing[0]

label = np.zeros_like(image, dtype=np.uint8)
for con in contours:
num = int(con['number'])
for c in con['contours']:
nodes = np.array(c).reshape((-1, 3))
assert np.amax(np.abs(np.diff(nodes[:, 2]))) == 0
z_index = z.index(nodes[0, 2])
r = (nodes[:, 1] - pos_r) / spacing_r
c = (nodes[:, 0] - pos_c) / spacing_c
rr, cc = polygon(r, c)
label[rr, cc, z_index] = num

colors = tuple(np.array([con['color'] for con in contours]) / 255.0)
return label, colors

train_data_path = "./DOI"
train_patients = [os.path.join(train_data_path, name)
for name in os.listdir(train_data_path) if os.path.isdir(os.path.join(train_data_path, name))]

patient = train_patients[0] # Just get the first patient for demo
for subdir, dirs, files in os.walk(patient):
dcms = glob.glob(os.path.join(subdir, "*.dcm"))
if len(dcms) == 1:
structure = dicom.read_file(os.path.join(subdir, files[0]))
contours = read_structure(structure)
elif len(dcms) > 1:
slices = [dicom.read_file(dcm) for dcm in dcms]
slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))
image = np.stack([s.pixel_array for s in slices], axis=-1)
label, colors = get_mask(contours, slices)

# Plot to check slices, for example 50 to 59
plt.figure(figsize=(15, 15))
for i in range(9):
plt.subplot(3, 3, i + 1)
plt.imshow(image[..., i + 50], cmap="gray")
plt.contour(label[..., i + 50], levels=[0.5, 1.5, 2.5, 3.5, 4.5], colors=colors)

Posted by: xuefeng @ May 23, 2017, 1:15 a.m.

Oops.. Seems posting ruins the indentation in the code snippet. Hope it will be easy to fix.

Posted by: xuefeng @ May 23, 2017, 1:17 a.m.


I know this is a little bit late but I'm using your program to get the contour label.
I found for most of the patients, this program works very well.
But for part of them it shows like "474.5 not exist" by using the get_mask function.

Have you ever met this problem during your processing?


Posted by: jh5442 @ Feb. 14, 2018, 2:23 a.m.

Nice code, indeed, but it does not work for LCTSC-Train-S3-004, LCTSC-Train-S3-005, LCTSC-Train-S3-007, LCTSC-Train-S3-008, and LCTSC-Train-S3-010. I will try with Slicer now as suggested in a previous thread of this forum http://aapmchallenges.cloudapp.net/forums/3/4/

Posted by: ghum @ Feb. 14, 2018, 2:41 p.m.


I also encountered this error and it seems to be some rounding problem. You might be able to avoid it by replacing "z_index = z.index(nodes[0, 2])" with:

zNew = [round(elem,1) for elem in z ]
z_index = z.index(nodes[0, 2])
except ValueError:
z_index = zNew.index(nodes[0, 2])

Posted by: blu @ Feb. 14, 2018, 3:11 p.m.

It works now, thanks!

Posted by: ghum @ Feb. 14, 2018, 3:16 p.m.

Thanks for the correction. Actually I found that error but forgot to update this post. Sorry for that.
For improved accuracy, I rounded the z positions to the first decimal like this:
z_index = z.index(np.around(nodes[0, 2], 1))
It doesn't matter for this dataset since the slice thickness is all larger than 1mm.

Posted by: xuefeng @ Feb. 14, 2018, 4:59 p.m.

I made a class for this which might make it a little easier, it takes a folder with dicom images and the RT structure and creates a mask for whatever ROI you request


Posted by: brianmanderson @ April 8, 2019, 7:13 p.m.

I try this code but it is not drawing contour for my slices , l also encountered with warning contour level not found .
please help me

Posted by: kushwah.rahulcse @ Oct. 22, 2019, 5:10 a.m.

To get this code working I needed to change a couple of things:

1- I also rounded z to one digit when it's first created in the first line of the function get_mask
2- Subdirs are found by os.walk in alphabetical order. This implies that if the folder with the RTSTRUCT doesn't have a name that preceeds in alphabetical order the name of the folder containing the dicom files of the CT, then the elif condition will be evaluated first, but the get_mask function inside this condition needs the variable contour that hasn't be created yet beacause it's folder will be explored later by os.walk. So I ordered these folders in alphabetical order before evaluating the conditions.

I put a minimum working example (for me) in this gist: https://gist.github.com/Feyn-Man/7191ce053286dbd1b67e986fdf9c851e

Posted by: aelius @ Dec. 12, 2019, 1:04 p.m.


I combined my previous work into an installable, please find it here

you can install it with
pip install DicomRTTool

from DicomRTTool import DicomReaderWriter

Contour_Names = ['Liver']
path = 'C:\users\brianmanderson\Patients\Patient_1\CT_1\'
Dicom_reader = DicomReaderWriter(get_images_mask=True, Contour_Names=Contour_Names)

image = DicomImage.ArrayDicom # [Images, rows, cols]
mask = DicomImage.mask # [Images, rows, cols, # classes + 1]

You can also write predictions out with

Dicom_reader.with_annotations(mask ,output_path,ROI_Names=['written_liver'])

Posted by: brianmanderson @ Sept. 22, 2020, 8:52 p.m.
Post in this thread