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]
contours.append(contour)
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)
plt.axis('off')

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.

Hi,

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?

Thanks,
Jin

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.

Hi,

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 ]
try:
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.
Post in this thread