-
Notifications
You must be signed in to change notification settings - Fork 14
/
Copy pathdvhcalc.py
208 lines (174 loc) · 7.66 KB
/
dvhcalc.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# dvhcalc.py
"""Calculate dose volume histogram (DVH) from DICOM RT Structure/Dose data."""
# Copyright (c) 2016 gluce
# Copyright (c) 2011-2016 Aditya Panchal
# Copyright (c) 2010 Roy Keyes
# This file is part of dicompyler-core, released under a BSD license.
# See the file license.txt included with this distribution, also
# available at https://github.com/dicompyler/dicompyler-core/
from __future__ import division
import pylab as pl
from dicompylercore import dicomparser, dvh, dvhcalc
import sys
import numpy as np
import numpy.ma as ma
import matplotlib.path
from six import iteritems
import logging
logger = logging.getLogger('dicompylercore.dvhcalc')
def get_dvh(structure, dose, roi, limit=None, callback=None):
"""Calculate a cumulative DVH in Gy from a DICOM RT Structure Set & Dose.
Parameters
----------
structure : pydicom Dataset
DICOM RT Structure Set used to determine the structure data.
dose : pydicom Dataset
DICOM RT Dose used to determine the dose grid.
roi : int
The ROI number used to uniquely identify the structure in the structure
set.
limit : int, optional
Dose limit in cGy as a maximum bin for the histogram.
callback : function, optional
A function that will be called at every iteration of the calculation.
"""
from dicompylercore import dicomparser
rtss = dicomparser.DicomParser(structure)
rtdose = dicomparser.DicomParser(dose)
structures = rtss.GetStructures()
s = structures[roi]
s['planes'] = rtss.GetStructureCoordinates(roi)
s['thickness'] = rtss.CalculatePlaneThickness(s['planes'])
hist = calculate_dvh(s, rtdose, limit, callback)
return dvh.DVH(counts=hist,
bins=(np.arange(0, 2) if (hist.size == 1) else
np.arange(0, hist.size + 1) / 100),
dvh_type='differential',
dose_units='gy',
name=s['name']
).cumulative
def calculate_dvh(structure, dose, limit=None, callback=None):
"""Calculate the differential DVH for the given structure and dose grid.
Parameters
----------
structure : dict
A structure (ROI) from an RT Structure Set parsed using DicomParser
dose : DicomParser
A DicomParser instance of an RT Dose
limit : int, optional
Dose limit in cGy as a maximum bin for the histogram.
callback : function, optional
A function that will be called at every iteration of the calculation.
"""
planes = structure['planes']
logger.debug(
"Calculating DVH of %s %s", structure['id'], structure['name'])
# Create an empty array of bins to store the histogram in cGy
# only if the structure has contour data or the dose grid exists
if ((len(planes)) and ("PixelData" in dose.ds)):
# Get the dose and image data information
dd = dose.GetDoseData()
id = dose.GetImageData()
# Generate a 2d mesh grid to create a polygon mask in dose coordinates
# Code taken from Stack Overflow Answer from Joe Kington:
# https://stackoverflow.com/q/3654289/74123
# Create vertex coordinates for each grid cell
x, y = np.meshgrid(np.array(dd['lut'][0]), np.array(dd['lut'][1]))
x, y = x.flatten(), y.flatten()
dosegridpoints = np.vstack((x, y)).T
maxdose = int(dd['dosemax'] * dd['dosegridscaling'] * 100)
# Remove values above the limit (cGy) if specified
if isinstance(limit, int):
if (limit < maxdose):
maxdose = limit
hist = np.zeros(maxdose)
else:
return np.array([0])
n = 0
planedata = {}
# Iterate over each plane in the structure
for z, plane in iteritems(planes):
# Get the dose plane for the current structure plane
doseplane = dose.GetDoseGrid(z)
planedata[z] = calculate_plane_histogram(
plane, doseplane, dosegridpoints,
maxdose, dd, id, structure, hist)
n += 1
if callback:
callback(n, len(planes))
# Volume units are given in cm^3
volume = sum([p[1] for p in planedata.values()]) / 1000
# Rescale the histogram to reflect the total volume
hist = sum([p[0] for p in planedata.values()])
hist = hist * volume / sum(hist)
# Remove the bins above the max dose for the structure
hist = np.trim_zeros(hist, trim='b')
return hist
def calculate_plane_histogram(plane, doseplane, dosegridpoints,
maxdose, dd, id, structure, hist):
"""Calculate the DVH for the given plane in the structure."""
contours = [[x[0:2] for x in c['data']] for c in plane]
# If there is no dose for the current plane, go to the next plane
if not len(doseplane):
return (np.arange(0, maxdose), 0)
# Create a zero valued bool grid
grid = np.zeros((dd['rows'], dd['columns']), dtype=np.uint8)
# Calculate the histogram for each contour in the plane
# and boolean xor to remove holes
for i, contour in enumerate(contours):
m = get_contour_mask(dd, id, dosegridpoints, contour)
grid = np.logical_xor(m.astype(np.uint8), grid).astype(np.bool)
hist, vol = calculate_contour_dvh(
grid, doseplane, maxdose, dd, id, structure)
return (hist, vol)
def get_contour_mask(dd, id, dosegridpoints, contour):
"""Get the mask for the contour with respect to the dose plane."""
doselut = dd['lut']
c = matplotlib.path.Path(list(contour))
grid = c.contains_points(dosegridpoints)
grid = grid.reshape((len(doselut[1]), len(doselut[0])))
return grid
def calculate_contour_dvh(mask, doseplane, maxdose, dd, id, structure):
"""Calculate the differential DVH for the given contour and dose plane."""
# Multiply the structure mask by the dose plane to get the dose mask
mask = ma.array(doseplane * dd['dosegridscaling'] * 100, mask=~mask)
# Calculate the differential dvh
hist, edges = np.histogram(mask.compressed(),
bins=maxdose,
range=(0, maxdose))
# Calculate the volume for the contour for the given dose plane
vol = sum(hist) * ((id['pixelspacing'][0]) *
(id['pixelspacing'][1]) *
(structure['thickness']))
return hist, vol
# ========================== Test DVH Calculation =========================== #
def main():
# Read the example RT structure and RT dose files
# The testdata was downloaded from the dicompyler website as testdata.zip
# Obtain the structures and DVHs from the DICOM data
rtssfile = 'testdata/rtss.dcm'
rtdosefile = 'testdata/rtdose.dcm'
RTss = dicomparser.DicomParser(rtssfile)
#RTdose = dicomparser.DicomParser("testdata/rtdose.dcm")
RTstructures = RTss.GetStructures()
# Generate the calculated DVHs
calcdvhs = {}
for key, structure in RTstructures.iteritems():
calcdvhs[key] = dvhcalc.get_dvh(rtssfile, rtdosefile, key)
if (key in calcdvhs) and (len(calcdvhs[key].counts) and calcdvhs[key].counts[0]!=0):
print ('DVH found for ' + structure['name'])
pl.plot(calcdvhs[key].counts * 100/calcdvhs[key].counts[0],
color=dvhcalc.np.array(structure['color'], dtype=float) / 255,
label=structure['name'],
linestyle='dashed')
#else:
# print("%d: no DVH"%key)
pl.xlabel('Distance (cm)')
pl.ylabel('Percentage Volume')
pl.legend(loc=7, borderaxespad=-5)
pl.setp(pl.gca().get_legend().get_texts(), fontsize='x-small')
pl.savefig('testdata/dvh.png', dpi = 75)
if __name__ == "__main__":
main()