This notebook aims to independently calculate the dose based on the dwell positions within an Oncentra® RT Dicom Plan, and compare the resulting dose grid to an exported Oncentra® dose grid within an RT Dicom Dose file.
Collect a range of brachytherapy dicom files that are able to be placed within the Gihub repository that can be used for testing. Aim to support as many brachy dicom formats as possible.
Once egs_brachy (https://doi.org/10.1088/0031-9155/61/23/8214) is available for use I would like to directly implement that model based dose calculation method to be able to provide checking for both TG43 and Monte Carlo based algorithms.
It is also planned that upper and lower 95% confidence interval doses will be reported when uncertainties due to catheter movement, catheter reconstruction, and calculation uncertainties are taken into account. Using the structure dicom file these 95% confidence interval doses can be converted to comparative DVHs.
The final stage is for this code to rerun a dwell time optimisation to see if any plan improvement is possible and compare the robustness of the original plan to positioning uncertainties with the calculated comparative plan.
Copyright (C) 2016-2017 Simon Biggs
This program is free software: you can redistribute it and/or modify it under the terms of the GNU Affero General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for more details.
You should have received a copy of the GNU Affero General Public License along with this program. If not, see http://www.gnu.org/licenses/.
Brachytherapy dicom files make use of a large number of private dicom tags. What some of these tags mean needs to be reverse engineered at times. This program has as of yet only been tested with Oncentra® dicom files, not BrachyVision™ dicom files. It has only been tested with a small subset of Oncentra® dicom files, with some of these files for certain configurations this code does not yet work. Be sure when using this code to investigate the testing figures produced to confirm that they represent what is expected within the plan.
Pay particular attention to x,y,z definitions, catheter definitions, and source orientation.
import numpy as np
import pandas as pd
import xarray as xr
from scipy.interpolate import RegularGridInterpolator
from copy import copy
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.collections import LineCollection
import matplotlib.gridspec as gridspec
%matplotlib inline
import itertools
from IPython.display import display, Markdown
source_date_directory = 'source_data/estro_website/'
# I couldn't seem to get Granero's data to work well
# source_date_directory = 'source_data/granero_paper/'
# I would like to also try Taylor's data
dose_rate_constant = np.squeeze(
pd.DataFrame.from_csv(source_date_directory + 'Lambda.csv', index_col=None, header=None).values)
dose_rate_constant
array(1.10942)
# When using the Estro data it appears that a length of 0.36 agrees better than a length of 0.35.
# I am worried that during the consensus data merging that data lookups that were based on a length of 0.36
# were potentially mixed with data lookups based on a length of 0.35.
length = np.squeeze(
pd.DataFrame.from_csv(source_date_directory + 'L.csv', index_col=None, header=None).values)
# length = 0.36
length
array(0.35)
radial_function_table = pd.DataFrame.from_csv(source_date_directory + 'g_L(r).csv')
radial_function_data = xr.DataArray(
np.squeeze(radial_function_table.values),
coords=[
radial_function_table.index.values.astype(float)],
dims=['radius_cm'])
radial_function_data
<xarray.DataArray (radius_cm: 17)> array([ 1.276, 1.199, 1.11 , 1.018, 1.001, 0.995, 0.997, 0.998, 1. , 1.003, 1.005, 1.008, 1.007, 1.003, 0.996, 0.972, 0.939]) Coordinates: * radius_cm (radius_cm) float64 0.06 0.08 0.1 0.15 0.2 0.25 0.5 0.75 1.0 ...
anisotropy_function_table = pd.DataFrame.from_csv(source_date_directory + 'F(r,theta).csv')
anisotropy_function_data = xr.DataArray(
anisotropy_function_table.values,
coords=[
anisotropy_function_table.index.values.astype(float),
anisotropy_function_table.columns.values.astype(float)],
dims=['theta_deg', 'radius_cm'])
anisotropy_function_data
<xarray.DataArray (theta_deg: 47, radius_cm: 19)> array([[ 0.951, 0.934, 0.917, ..., 0.733, 0.768, 0.798], [ 0.947, 0.93 , 0.914, ..., 0.744, 0.775, 0.801], [ 0.944, 0.927, 0.91 , ..., 0.755, 0.782, 0.804], ..., [ 0.575, 0.575, 0.576, ..., 0.718, 0.749, 0.773], [ 0.536, 0.537, 0.537, ..., 0.682, 0.717, 0.746], [ 0.497, 0.498, 0.499, ..., 0.646, 0.685, 0.718]]) Coordinates: * theta_deg (theta_deg) float64 0.0 2.0 4.0 6.0 8.0 10.0 12.0 14.0 16.0 ... * radius_cm (radius_cm) float64 0.06 0.08 0.1 0.15 0.2 0.25 0.3 0.35 0.4 ...
QA_along_away_table = pd.DataFrame.from_csv(source_date_directory + 'QA_along_away.csv')
QA_along_away_data = xr.DataArray(
QA_along_away_table.values.astype(float),
coords=[
QA_along_away_table.index.values.astype(float),
QA_along_away_table.columns.values.astype(float)],
dims=['z_cm', 'y_cm'])
QA_along_away_data
<xarray.DataArray (z_cm: 19, y_cm: 12)> array([[ 0.0169 , 0.01709, 0.01724, ..., 0.01385, 0.01209, 0.01044], [ 0.0227 , 0.0231 , 0.0235 , ..., 0.01719, 0.01455, 0.01226], [ 0.032 , 0.0328 , 0.0337 , ..., 0.0213 , 0.0175 , 0.01432], ..., [ 0.0279 , 0.0304 , 0.0324 , ..., 0.0213 , 0.0175 , 0.01432], [ 0.02 , 0.0213 , 0.0225 , ..., 0.01717, 0.01455, 0.01226], [ 0.015 , 0.01576, 0.01647, ..., 0.01384, 0.01207, 0.01044]]) Coordinates: * z_cm (z_cm) float64 7.0 6.0 5.0 4.0 3.0 2.0 1.5 1.0 0.5 0.0 -0.5 ... * y_cm (y_cm) float64 0.0 0.25 0.5 0.75 1.0 1.5 2.0 3.0 4.0 5.0 6.0 7.0
TG43U1S1 states:
...a log-linear function for $g(r)$ interpolation was recommended in the 2004 AAPM TG-43U1 report. ... The log-linear interpolation should be performed using data points immediately adjacent to the radius of interest.
...
Appendix C 3 of the 2004 AAPM TG-43U1 report clearly
specified an extrapolation method (nearest neighbor or zeroth order) for 1D dose rate distributions when $r<r_{min}$ for $g(r_{min})$. Due to the great variability in $g(r)$ based on choice of L and features of source construction, use of nearestneighbor or zeroth-order data is still recommended for extrapolation of $g_L(r)$ for $r<r_{min}$.
...
The AAPM now recommends adoption of a single exponential
function based on fitting $g_L(r)$ data points for the largest two consensus $r$ values ... Specifically, a log-linear extrapolation ...
plt.scatter(
radial_function_data.radius_cm,
radial_function_data)
plt.xlabel("Radius (cm)")
plt.ylabel("Radial dose function")
<matplotlib.text.Text at 0xd078390>
plt.scatter(
np.log(radial_function_data.radius_cm[1::]),
radial_function_data[1::])
plt.xlabel("$ln$( Radius (cm) )")
plt.ylabel("Radial dose function")
<matplotlib.text.Text at 0xd45ebe0>
radial_function_interpolation = RegularGridInterpolator(
(np.log(radial_function_data.radius_cm),),
radial_function_data,
bounds_error=False,
fill_value=None
)
def radial_function(radius):
interpolation = radial_function_interpolation(np.log(radius))
too_close_ref = radius < np.array(radial_function_data.radius_cm[0])
interpolation[too_close_ref] = radial_function_data[0]
return interpolation
r_test = np.linspace(0.03, 30, 10000)
y_test = radial_function(r_test)
plt.figure()
plt.scatter(
np.log(radial_function_data.radius_cm),
radial_function_data)
plt.plot(
np.log(r_test),
y_test)
plt.xlabel("$ln$( Radius (cm) )")
plt.ylabel("Radial dose function")
plt.figure()
plt.scatter(
radial_function_data.radius_cm,
radial_function_data)
plt.plot(
r_test,
y_test)
plt.xlabel("Radius (cm)")
plt.ylabel("Radial dose function")
<matplotlib.text.Text at 0xd588d68>
TG43U1S1 states:
"...linear-linear interpolation method for $F(r,\theta)$, as a function of $r$ and $\theta$ is appropriate, and should be based on the two data points for each variable located immediately adjacent to the interpolated point of interest."
"...the nearestneighbor or zeroth-order approach presented in Appendix C of the 2004 AAPM TG-43U1 report is still recommended for $F(r,\theta)$ extrapolation for $r < r_{min}$ and also for $r > r_{max}$."
radius_cm_mesh, theta_deg_mesh = np.meshgrid(
anisotropy_function_data.radius_cm,
anisotropy_function_data.theta_deg)
plt.scatter(
radius_cm_mesh,
theta_deg_mesh,
c=anisotropy_function_data,
s=2)
c = plt.colorbar()
plt.xlabel("Radius (cm)")
plt.ylabel("Theta (degrees)")
c.set_label("Anisotropy function")
find_nan = np.where(np.isnan(anisotropy_function_data))
row = find_nan[0]
column = find_nan[1]
unique_row = np.unique(row)
for row_val in unique_row:
ref = row == row_val
last_nan = np.max(column[ref])
given_row = anisotropy_function_data[row_val, :]
nan_ref = np.isnan(given_row)
anisotropy_function_data[row_val, nan_ref] = given_row[last_nan + 1]
plt.scatter(
radius_cm_mesh,
theta_deg_mesh,
c=anisotropy_function_data,
s=2)
c = plt.colorbar()
plt.xlabel("Radius (cm)")
plt.ylabel("Theta (degrees)")
c.set_label("Anisotropy function")
marker_cycle = itertools.cycle((',', 'o', '*', '^'))
colour_list = [
item['color']
for item in mpl.rcParams['axes.prop_cycle']
]
colour_cycle = itertools.cycle(colour_list)
for data in anisotropy_function_data.T:
plt.scatter(
data.theta_deg, data, label="radius = {} cm".format(np.array(data.radius_cm)),
marker=next(marker_cycle), c=next(colour_cycle))
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.xlabel('Theta (degrees)')
plt.ylabel('Anisotropy function')
<matplotlib.text.Text at 0xe7a1ef0>
anisotropy_theta_deg = np.array(anisotropy_function_data.theta_deg)
anisotropy_radius_cm = np.array(anisotropy_function_data.radius_cm)
anisotropy_function_data_array = np.array(anisotropy_function_data)
anisotropy_function_linear_interpolation = RegularGridInterpolator(
(anisotropy_theta_deg, anisotropy_radius_cm), anisotropy_function_data_array,
bounds_error=False,
fill_value=np.nan)
anisotropy_function_nearest_interpolation = RegularGridInterpolator(
(anisotropy_theta_deg, anisotropy_radius_cm), anisotropy_function_data_array,
bounds_error=False,
fill_value=None,
method='nearest')
def anisotropy_function(radius, theta):
points = np.vstack([theta, radius]).T
linear_interpolation = anisotropy_function_linear_interpolation(points)
ref = np.isnan(linear_interpolation)
nearest_interpolation = anisotropy_function_nearest_interpolation(
points[ref,:])
linear_interpolation[ref] = nearest_interpolation
return linear_interpolation
r_test = np.linspace(0, 20, 100)
theta_test = np.linspace(0, 180, 100)
r_test_mesh, theta_test_mesh = np.meshgrid(
r_test, theta_test)
r_test_mesh_flat = r_test_mesh.ravel()
theta_test_mesh_flat = theta_test_mesh.ravel()
anisotropy_test = anisotropy_function(r_test_mesh_flat, theta_test_mesh_flat)
anisotropy_test = np.reshape(anisotropy_test, np.shape(r_test_mesh))
plt.contourf(r_test, theta_test, anisotropy_test, 1000)
c = plt.colorbar()
plt.xlabel("Radius (cm)")
plt.ylabel("Theta (degrees)")
c.set_label("Anisotropy function")
marker_cycle = itertools.cycle((',', 'o', '*', '^'))
colour_list = [
item['color']
for item in mpl.rcParams['axes.prop_cycle']
]
colour_cycle = itertools.cycle(colour_list)
theta_linspace = np.linspace(
np.array(np.min(anisotropy_function_data.theta_deg)),
np.array(np.max(anisotropy_function_data.theta_deg)), 1000)
for data in anisotropy_function_data.T:
colour = next(colour_cycle)
plt.scatter(
data.theta_deg, data, label="radius = {} cm".format(np.array(data.radius_cm)),
marker=next(marker_cycle), c=colour)
interpolation = anisotropy_function(np.array(data.radius_cm)*np.ones_like(theta_linspace), theta_linspace)
plt.plot(theta_linspace, interpolation, c=colour)
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.xlabel('Theta (degrees)')
plt.ylabel('Anisotropy function')
<matplotlib.text.Text at 0xfabe710>
r_test_mesh, theta_test_mesh = np.meshgrid(
anisotropy_radius_cm, anisotropy_theta_deg)
r_test_mesh_flat = r_test_mesh.ravel()
theta_test_mesh_flat = theta_test_mesh.ravel()
anisotropy_test = anisotropy_function(r_test_mesh_flat, theta_test_mesh_flat)
anisotropy_test = np.reshape(anisotropy_test, np.shape(r_test_mesh))
assert np.all(anisotropy_test - anisotropy_function_data_array == 0)
To better understand the geometry function see the following interactive demonstration:
The method here for calculating the geometry function is based on the following paper:
King, R. P., Anderson, R. S. and Mills, M. D. (2001), Geometry function of a linear brachytherapy source. Journal of Applied Clinical Medical Physics, 2: 69–72. doi:10.1120/jacmp.v2i2.2615
def calculate_beta(length, radius, theta):
# https://dx.doi.org/10.1120/jacmp.v2i2.2615
theta2 = np.arctan(
radius * np.sin(theta) /
(radius * np.cos(theta) - length/2))
AP = radius * np.sin(theta)
SA = radius * np.cos(theta) + length/2
SP = np.sqrt(AP**2 + SA**2)
beta = np.arcsin(
length * np.sin(theta2) /
SP)
return beta
def geometry_function(length, radius, theta):
theta = theta / 180 * np.pi
result = np.nan * np.ones_like(radius)
ref0 = (theta == 0) | (theta == np.pi)
ref1 = (theta != 0) & (theta != np.pi)
result[ref0] = (
1 / (radius[ref0]**2 - length**2 / 4))
beta = calculate_beta(
length, radius[ref1], theta[ref1])
result[ref1] = (
beta / (length * radius[ref1] * np.sin(theta[ref1])))
return np.abs(result)
# def geometry_function(length, radius, theta):
# return radius**-2
def normalised_geometry_function(length, radius, theta):
geometry = geometry_function(length, radius, theta)
ref = radius == 0
geometry_0 = geometry_function(length, np.array([0.000001]*np.sum(ref)), theta[ref])
geometry[ref] = geometry_0
return geometry / geometry_function(length, np.array([1]), np.array([90]))
def tg43(radius, theta):
initial_shape = np.shape(radius)
radius_flat = np.ravel(radius)
theta_flat = np.ravel(theta)
geometry = normalised_geometry_function(length, radius_flat, theta_flat)
radial = radial_function(radius_flat)
anisotropy = anisotropy_function(radius_flat, theta_flat)
result = dose_rate_constant * geometry * radial * anisotropy
result = np.reshape(result, initial_shape)
return result
def calc_on_grid(calc_grid_positions, dwell_positions, dwell_directions):
calc_grid_vector = calc_grid_positions[:,None,:] - dwell_positions[None,:,:]
radius = np.sqrt(
calc_grid_vector[:,:,0]**2 +
calc_grid_vector[:,:,1]**2 +
calc_grid_vector[:,:,2]**2)
calc_unit_vector = calc_grid_vector / radius[:,:,None]
dot_product = (
dwell_directions[None,:,0] * calc_unit_vector[:,:,0] +
dwell_directions[None,:,1] * calc_unit_vector[:,:,1] +
dwell_directions[None,:,2] * calc_unit_vector[:,:,2])
theta = np.arccos(dot_product) * 180 / np.pi
nan_ref = np.isnan(theta)
theta[nan_ref] = 90
tg43_dose = tg43(radius, theta)
return tg43_dose, radius
def determine_calc_grid_positions(calc_x, calc_y, calc_z):
calc_x_mesh, calc_y_mesh, calc_z_mesh = np.meshgrid(
calc_x, calc_y, calc_z)
calc_x_mesh_flat = calc_x_mesh.ravel()
calc_y_mesh_flat = calc_y_mesh.ravel()
calc_z_mesh_flat = calc_z_mesh.ravel()
calc_grid_positions = np.array([
[calc_x_mesh_flat[i], calc_y_mesh_flat[i], calc_z_mesh_flat[i]]
for i in range(len(calc_x_mesh_flat))
])
return calc_grid_positions, np.shape(calc_x_mesh)
def tg43_on_grid(calc_x, calc_y, calc_z,
reference_air_kerma_rate, dwell_times, dwell_positions, dwell_directions,
return_radius=False):
calc_grid_positions, shape = determine_calc_grid_positions(
calc_x, calc_y, calc_z)
tg43_dose, radius = calc_on_grid(
calc_grid_positions, dwell_positions, dwell_directions)
total_dose = np.sum(
tg43_dose * reference_air_kerma_rate * dwell_times[None, :], axis=1)
# too_close = np.any(radius < length/2, axis=1)
# total_dose[too_close] = np.nan
total_dose = np.reshape(total_dose, shape)
if return_radius:
result = (total_dose, np.reshape(radius, shape))
else:
result = total_dose
return result
testing_calc_x = np.array(0)
testing_calc_y = np.array(QA_along_away_data.y_cm)
testing_calc_z = np.array(QA_along_away_data.z_cm)
test_reference_air_kerma_rate = 1
test_dwell_times = np.array([1])
testing_dwell_positions = np.array([[0, 0, 0]])
testing_dwell_directions = np.array([[0, 0, 1]])
testing_tg43_dose, testing_radius = tg43_on_grid(
testing_calc_x, testing_calc_y, testing_calc_z,
test_reference_air_kerma_rate, test_dwell_times,
testing_dwell_positions, testing_dwell_directions, return_radius=True)
C:\Users\sbiggs\AppData\Local\Continuum\Anaconda3\lib\site-packages\ipykernel\__main__.py:10: RuntimeWarning: invalid value encountered in true_divide C:\Users\sbiggs\AppData\Local\Continuum\Anaconda3\lib\site-packages\ipykernel\__main__.py:31: RuntimeWarning: invalid value encountered in true_divide C:\Users\sbiggs\AppData\Local\Continuum\Anaconda3\lib\site-packages\ipykernel\__main__.py:9: RuntimeWarning: divide by zero encountered in log C:\Users\sbiggs\AppData\Local\Continuum\Anaconda3\lib\site-packages\scipy\interpolate\interpolate.py:2444: RuntimeWarning: invalid value encountered in add values += np.asarray(self.values[edge_indices]) * weight[vslice]
test_data = np.array(QA_along_away_data.T)
test_data[test_data > 1000] = np.nan
mpl.rcParams['contour.negative_linestyle'] = 'solid'
dose_difference = 100 * (testing_tg43_dose[:,0,:] - test_data) / test_data
levels = np.arange(-0.8,0.8,0.1)
CS = plt.contour(testing_calc_z, testing_calc_y, dose_difference, levels, colors='k', linewidths=0.2)
plt.clabel(CS,
inline=1,
fmt='%1.1f',
fontsize=9)
plt.contourf(testing_calc_z, testing_calc_y, dose_difference, np.arange(-0.8,0.8,0.01), extend='both', cmap='Spectral_r')
c = plt.colorbar()
c.set_label('Percent relative dose difference')
plt.xlabel('z (cm)')
plt.ylabel('y (cm)')
<matplotlib.text.Text at 0xfa8e320>
np.nanmax(np.abs(dose_difference))
0.92238097135892283
reference = np.squeeze(testing_radius > 0.6)
np.nanmax(np.abs(dose_difference[reference]))
0.40067723212877082
import dicom
from glob import glob
from scipy.interpolate import make_interp_spline
dose_filename = glob("private/RD1*")[0]
plan_filename = glob("private/RP1*")[0]
# dose_filename = "test_dicom/RD_y5cm.dcm"
# plan_filename = "test_dicom/RP_y5cm.dcm"
dcm_dose = dicom.read_file(dose_filename, force=True)
dcm_plan = dicom.read_file(plan_filename, force=True)
dcm_plan.SourceSequence[0]
(300a, 0212) Source Number IS: '0' (300a, 0214) Source Type CS: 'LINE' (300a, 0226) Source Isotope Name LO: 'Ir-192' (300a, 0228) Source Isotope Half Life DS: '73.8300000000000' (300a, 0229) Source Strength Units CS: 'AIR_KERMA_RATE' (300a, 022a) Reference Air Kerma Rate DS: '43939.3000000000' (300a, 022c) Source Strength Reference Date DA: '20151201' (300a, 022e) Source Strength Reference Time TM: '000000' (300b, 0010) Private Creator LO: 'NUCLETRON' (300b, 1006) Private tag data DS: '-100557' (300b, 1007) Private tag data DS: '-100631' (300b, 1008) Private tag data DT: '20151130140000' (300b, 100c) Private tag data DS: ''
reference_air_kerma_rate = np.float(dcm_plan.SourceSequence[0].ReferenceAirKermaRate) / 360000
reference_air_kerma_rate
0.12205361111111111
dcm_plan.ApplicationSetupSequence[0].ChannelSequence[0]
(3006, 0084) Referenced ROI Number IS: '' (300a, 0110) Number of Control Points IS: '26' (300a, 0282) Channel Number IS: '1' (300a, 0284) Channel Length DS: '1239.00000000000' (300a, 0286) Channel Total Time DS: '22.0807651173517' (300a, 0288) Source Movement Type CS: 'STEPWISE' (300a, 0290) Source Applicator Number IS: '1' (300a, 0291) Source Applicator ID SH: '1' (300a, 0292) Source Applicator Type CS: 'FLEXIBLE' (300a, 0294) Source Applicator Name LO: '' (300a, 0296) Source Applicator Length DS: '158.898261747671' (300a, 02a0) Source Applicator Step Size DS: '2.50000000000000' (300a, 02a2) Transfer Tube Number IS: '1' (300a, 02a4) Transfer Tube Length DS: '' (300a, 02c8) Final Cumulative Time Weight DS: '5.07299481332302' (300a, 02d0) Brachy Control Point Sequence 26 item(s) ---- (300a, 0112) Control Point Index IS: '0' (300a, 02d2) Control Point Relative Position DS: '30.0000000000000' (300a, 02d4) Control Point 3D Position DS: ['-21.196375', '8.158049', '36.385686'] (300a, 02d6) Cumulative Time Weight DS: '0.00000000000000' --------- (300a, 0112) Control Point Index IS: '1' (300a, 02d2) Control Point Relative Position DS: '30.0000000000000' (300a, 02d4) Control Point 3D Position DS: ['-21.196375', '8.158049', '36.385686'] (300a, 02d6) Cumulative Time Weight DS: '1.02102994918823' --------- (300a, 0112) Control Point Index IS: '2' (300a, 02d2) Control Point Relative Position DS: '32.5000000000000' (300a, 02d4) Control Point 3D Position DS: ['-20.892436', '8.876451', '34.010497'] (300a, 02d6) Cumulative Time Weight DS: '1.02102994918823' --------- (300a, 0112) Control Point Index IS: '3' (300a, 02d2) Control Point Relative Position DS: '32.5000000000000' (300a, 02d4) Control Point 3D Position DS: ['-20.892436', '8.876451', '34.010497'] (300a, 02d6) Cumulative Time Weight DS: '1.83374685049057' --------- (300a, 0112) Control Point Index IS: '4' (300a, 02d2) Control Point Relative Position DS: '35.0000000000000' (300a, 02d4) Control Point 3D Position DS: ['-20.588497', '9.594852', '31.635308'] (300a, 02d6) Cumulative Time Weight DS: '1.83374685049057' --------- (300a, 0112) Control Point Index IS: '5' (300a, 02d2) Control Point Relative Position DS: '35.0000000000000' (300a, 02d4) Control Point 3D Position DS: ['-20.588497', '9.594852', '31.635308'] (300a, 02d6) Cumulative Time Weight DS: '1.95902469754219' --------- (300a, 0112) Control Point Index IS: '6' (300a, 02d2) Control Point Relative Position DS: '37.5000000000000' (300a, 02d4) Control Point 3D Position DS: ['-20.293912', '10.336907', '29.266542'] (300a, 02d6) Cumulative Time Weight DS: '1.95902469754219' --------- (300a, 0112) Control Point Index IS: '7' (300a, 02d2) Control Point Relative Position DS: '37.5000000000000' (300a, 02d4) Control Point 3D Position DS: ['-20.293912', '10.336907', '29.266542'] (300a, 02d6) Cumulative Time Weight DS: '2.15143331885338' --------- (300a, 0112) Control Point Index IS: '8' (300a, 02d2) Control Point Relative Position DS: '40.0000000000000' (300a, 02d4) Control Point 3D Position DS: ['-20.020003', '11.131245', '26.911972'] (300a, 02d6) Cumulative Time Weight DS: '2.15143331885338' --------- (300a, 0112) Control Point Index IS: '9' (300a, 02d2) Control Point Relative Position DS: '40.0000000000000' (300a, 02d4) Control Point 3D Position DS: ['-20.020003', '11.131245', '26.911972'] (300a, 02d6) Cumulative Time Weight DS: '2.44192209839821' --------- (300a, 0112) Control Point Index IS: '10' (300a, 02d2) Control Point Relative Position DS: '42.5000000000000' (300a, 02d4) Control Point 3D Position DS: ['-19.746093', '11.925583', '24.557402'] (300a, 02d6) Cumulative Time Weight DS: '2.44192209839821' --------- (300a, 0112) Control Point Index IS: '11' (300a, 02d2) Control Point Relative Position DS: '42.5000000000000' (300a, 02d4) Control Point 3D Position DS: ['-19.746093', '11.925583', '24.557402'] (300a, 02d6) Cumulative Time Weight DS: '2.64625501632690' --------- (300a, 0112) Control Point Index IS: '12' (300a, 02d2) Control Point Relative Position DS: '45.0000000000000' (300a, 02d4) Control Point 3D Position DS: ['-19.465747', '12.727309', '22.206254'] (300a, 02d6) Cumulative Time Weight DS: '2.64625501632690' --------- (300a, 0112) Control Point Index IS: '13' (300a, 02d2) Control Point Relative Position DS: '45.0000000000000' (300a, 02d4) Control Point 3D Position DS: ['-19.465747', '12.727309', '22.206254'] (300a, 02d6) Cumulative Time Weight DS: '2.79906022548676' --------- (300a, 0112) Control Point Index IS: '14' (300a, 02d2) Control Point Relative Position DS: '47.5000000000000' (300a, 02d4) Control Point 3D Position DS: ['-19.140841', '13.580187', '19.878802'] (300a, 02d6) Cumulative Time Weight DS: '2.79906022548676' --------- (300a, 0112) Control Point Index IS: '15' (300a, 02d2) Control Point Relative Position DS: '47.5000000000000' (300a, 02d4) Control Point 3D Position DS: ['-19.140841', '13.580187', '19.878802'] (300a, 02d6) Cumulative Time Weight DS: '3.51364403963089' --------- (300a, 0112) Control Point Index IS: '16' (300a, 02d2) Control Point Relative Position DS: '50.0000000000000' (300a, 02d4) Control Point 3D Position DS: ['-18.815936', '14.433065', '17.551349'] (300a, 02d6) Cumulative Time Weight DS: '3.51364403963089' --------- (300a, 0112) Control Point Index IS: '17' (300a, 02d2) Control Point Relative Position DS: '50.0000000000000' (300a, 02d4) Control Point 3D Position DS: ['-18.815936', '14.433065', '17.551349'] (300a, 02d6) Cumulative Time Weight DS: '3.54612701758742' --------- (300a, 0112) Control Point Index IS: '18' (300a, 02d2) Control Point Relative Position DS: '52.5000000000000' (300a, 02d4) Control Point 3D Position DS: ['-18.545826', '15.345881', '15.239715'] (300a, 02d6) Cumulative Time Weight DS: '3.54612701758742' --------- (300a, 0112) Control Point Index IS: '19' (300a, 02d2) Control Point Relative Position DS: '52.5000000000000' (300a, 02d4) Control Point 3D Position DS: ['-18.545826', '15.345881', '15.239715'] (300a, 02d6) Cumulative Time Weight DS: '3.90212383493781' --------- (300a, 0112) Control Point Index IS: '20' (300a, 02d2) Control Point Relative Position DS: '55.0000000000000' (300a, 02d4) Control Point 3D Position DS: ['-18.276953', '16.260050', '12.928438'] (300a, 02d6) Cumulative Time Weight DS: '3.90212383493781' --------- (300a, 0112) Control Point Index IS: '21' (300a, 02d2) Control Point Relative Position DS: '55.0000000000000' (300a, 02d4) Control Point 3D Position DS: ['-18.276953', '16.260050', '12.928438'] (300a, 02d6) Cumulative Time Weight DS: '4.02017613127828' --------- (300a, 0112) Control Point Index IS: '22' (300a, 02d2) Control Point Relative Position DS: '57.5000000000000' (300a, 02d4) Control Point 3D Position DS: ['-18.008080', '17.174219', '10.617161'] (300a, 02d6) Cumulative Time Weight DS: '4.02017613127828' --------- (300a, 0112) Control Point Index IS: '23' (300a, 02d2) Control Point Relative Position DS: '57.5000000000000' (300a, 02d4) Control Point 3D Position DS: ['-18.008080', '17.174219', '10.617161'] (300a, 02d6) Cumulative Time Weight DS: '4.07139252126217' --------- (300a, 0112) Control Point Index IS: '24' (300a, 02d2) Control Point Relative Position DS: '62.5000000000000' (300a, 02d4) Control Point 3D Position DS: ['-17.523745', '19.171616', '6.059713'] (300a, 02d6) Cumulative Time Weight DS: '4.07139252126217' --------- (300a, 0112) Control Point Index IS: '25' (300a, 02d2) Control Point Relative Position DS: '62.5000000000000' (300a, 02d4) Control Point 3D Position DS: ['-17.523745', '19.171616', '6.059713'] (300a, 02d6) Cumulative Time Weight DS: '5.07299481332302' --------- (300b, 0010) Private Creator LO: 'NUCLETRON' (300b, 1000) Private tag data IS: '0' (300b, 100e) Private tag data DS: '1' (300b, 1013) Private tag data LO: 'NO_MODEL_ID' (300b, 1027) Private tag data DS: '6.00000000000000' (300b, 1032) Private tag data DS: '0.00000000000000' (300c, 000e) Referenced Source Number IS: '0'
def pull_dwells(dcm):
dwell_channels = []
all_dwell_positions = []
dwell_times = []
number_of_channels = len(dcm.ApplicationSetupSequence[0].ChannelSequence)
for i in range(number_of_channels):
ChannelSequence = dcm.ApplicationSetupSequence[0].ChannelSequence[i]
if len(ChannelSequence.dir("FinalCumulativeTimeWeight")) != 0:
BrachyControlPointSequence = (
ChannelSequence.BrachyControlPointSequence)
number_of_dwells_in_channel = len(BrachyControlPointSequence)//2
channel_time_weight = (
ChannelSequence.ChannelTotalTime / ChannelSequence.FinalCumulativeTimeWeight)
channel_dwells_positions = []
for j in range(0,number_of_dwells_in_channel*2,2):
if len(BrachyControlPointSequence[j].dir("ControlPoint3DPosition")) != 0:
channel_dwells_positions.append(
BrachyControlPointSequence[j].ControlPoint3DPosition)
dwell_channels.append(i+1)
assert (
BrachyControlPointSequence[j].ControlPoint3DPosition ==
BrachyControlPointSequence[j+1].ControlPoint3DPosition
)
uncorrected_dwell_time = (
BrachyControlPointSequence[j+1].CumulativeTimeWeight -
BrachyControlPointSequence[j].CumulativeTimeWeight
)
dwell_times.append(
uncorrected_dwell_time * channel_time_weight)
all_dwell_positions += channel_dwells_positions
dwell_channels = np.array(dwell_channels)
dwell_positions = np.array(all_dwell_positions).astype(float) / 10 # convert from mm to cm
dwell_times = np.array(dwell_times)
return dwell_positions, dwell_channels, dwell_times
dwell_positions, dwell_channels, dwell_times = pull_dwells(dcm_plan)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
scatt = ax.scatter(
dwell_positions[:,0],
dwell_positions[:,1],
dwell_positions[:,2], c=dwell_times, s=dwell_times*15)
c = plt.colorbar(scatt)
c.set_label('Dwell time (s)')
ax.set(xlabel='x (cm)', ylabel='y (cm)', zlabel='z (cm)')
# plt.figure()
# plt.scatter(dwell_positions[:,0], dwell_positions[:,1])
[<matplotlib.text.Text at 0x104f15c0>, <matplotlib.text.Text at 0x104e9ba8>, <matplotlib.text.Text at 0x104e5080>]
def pull_catheter_coords(dcm, channel):
index = channel - 1
# This line isn't correct for every dicom file. Specifically the "index" as given here
# doesn't uniquely define that given catheter channel. There is a specific tag within
# the contour sequence which should be used instead.
# The risk is that there is a single contour sequence that doesn't refer to a catheter
# instead I suspect that is where info such as catheter storage direction and the like
# is stored. It just so happens in this file that this private sequence is at the end
# in this example file. That is not always the case.
catheter_coords = dcm[0x300f,0x1000][0].ROIContourSequence[index].ContourSequence[0].ContourData
# mm to cm
x = np.array(catheter_coords[0::3]) / 10
y = np.array(catheter_coords[1::3]) / 10
z = np.array(catheter_coords[2::3]) / 10
return x, y, z
catheter_x, catheter_y, catheter_z = pull_catheter_coords(dcm_plan, 1)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(
catheter_x,
catheter_y,
catheter_z)
ax.set(xlabel='x', ylabel='y', zlabel='z')
[<matplotlib.text.Text at 0x119dbd68>, <matplotlib.text.Text at 0x119e35c0>, <matplotlib.text.Text at 0xd6dd358>]
print(dcm_plan[0x300f,0x1000][0].ROIContourSequence[0])
(0021, 0010) Private Creator LO: 'NUCLETRON' (0021, 1000) Private tag data UI: 1.2.840.113619.2.278.3.380434001.123.1424034575.743 (3006, 002a) ROI Display Color IS: ['255', '0', '0'] (3006, 0040) Contour Sequence 1 item(s) ---- (3006, 0042) Contour Geometric Type CS: 'OPEN_NONPLANAR' (3006, 0046) Number of Contour Points IS: '16' (3006, 0048) Contour Number IS: '1' (3006, 0050) Contour Data DS: ['-14.505572', '70.248421', '-75.000000', '-14.271858', '68.371307', '-70.000000', '-13.748368', '61.304198', '-57.500000', '-13.486623', '52.666620', '-45.000000', '-14.097361', '42.371325', '-32.500000', '-14.969844', '33.123009', '-20.000000', '-16.016823', '27.015631', '-10.000000', '-16.889306', '21.867983', '0.000000', '-17.936285', '17.418322', '10.000000', '-18.808767', '14.451881', '17.500000', '-19.506753', '12.619668', '22.500000', '-20.379236', '10.089468', '30.000000', '-21.338967', '7.821013', '37.500000', '-22.298698', '5.639807', '45.000000', '-23.258429', '3.284104', '52.500000', '-24.916146', '-0.031330', '65.000000'] --------- (3006, 0084) Referenced ROI Number IS: '0'
def dwell_direction_from_catheter_coords(catheter_x, catheter_y, catheter_z,
dwell_x, dwell_y, dwell_z):
t = np.linspace(0, 1, len(catheter_x))
spl = make_interp_spline(t, np.c_[
catheter_x, catheter_y, catheter_z
], k=1)
t_find_closest = np.linspace(0,1,10000)
x_find_closest, y_find_closest, z_find_closest = spl(t_find_closest).T
distance = np.sqrt(
(dwell_x[None,:] - x_find_closest[:,None])**2 +
(dwell_y[None,:] - y_find_closest[:,None])**2 +
(dwell_z[None,:] - z_find_closest[:,None])**2
)
dwell_params = t_find_closest[np.argmin(distance, axis=0)]
dwell_derivs = spl.derivative()(dwell_params).T
# dwell_directions = -dwell_derivs / np.linalg.norm(dwell_derivs, axis=0)
dwell_directions = dwell_derivs / np.linalg.norm(dwell_derivs, axis=0)
return dwell_directions.T
def determine_dwell_directions(dcm, dwell_positions, dwell_channels):
catheter_numbers = np.unique(dwell_channels)
dwell_directions = np.nan * np.ones_like(dwell_positions)
for catheter_number in catheter_numbers:
catheter_x, catheter_y, catheter_z = pull_catheter_coords(dcm, catheter_number)
current_ref = dwell_channels==catheter_number
current_dwell_positions = dwell_positions[current_ref,:]
current_dwell_directions = dwell_direction_from_catheter_coords(
catheter_x, catheter_y, catheter_z,
current_dwell_positions[:,0], current_dwell_positions[:,1],
current_dwell_positions[:,2])
dwell_directions[current_ref, :] = current_dwell_directions
return dwell_directions
dwell_directions = determine_dwell_directions(dcm_plan, dwell_positions, dwell_channels)
# Need to make a dicom file specific correction based off end/catheter or catheter/end
# private dicom tag
# dwell_directions = -dwell_directions
plt.plot(dwell_directions)
[<matplotlib.lines.Line2D at 0x11ab0278>, <matplotlib.lines.Line2D at 0x11ab0518>, <matplotlib.lines.Line2D at 0x11ab0710>]
def show_dwell_directions(catheter_x, catheter_y, catheter_z,
dwell_x, dwell_y, dwell_z, dwell_directions):
fig = plt.figure(figsize=(10, 10))
gs = gridspec.GridSpec(2, 2)
ax3 = fig.add_subplot(gs[1,0])
ax2 = fig.add_subplot(gs[0,1], sharex=ax3)
ax1 = fig.add_subplot(gs[0,0], sharey=ax2, sharex=ax3)
# ax1 = fig.add_subplot(gs[0,0])
# ax2 = fig.add_subplot(gs[0,1], sharey=ax1)
# ax3 = fig.add_subplot(gs[1,0], sharex=ax1)
ax4 = fig.add_subplot(gs[1,1], projection='3d')
ax1.plot(catheter_x, catheter_z, '--', alpha=0.2)
ax2.plot(catheter_y, catheter_z, '--', alpha=0.2)
ax3.plot(catheter_x, catheter_y, '--', alpha=0.2)
ax1.scatter(dwell_x, dwell_z)
ax2.scatter(dwell_y, dwell_z)
ax3.scatter(dwell_x, dwell_y)
ax1.quiver(
dwell_x,
dwell_z,
dwell_directions[:,0], dwell_directions[:,2],
pivot='mid', scale=10
)
ax2.quiver(
dwell_y,
dwell_z,
dwell_directions[:,1], dwell_directions[:,2],
pivot='mid', scale=10
)
ax3.quiver(
dwell_x,
dwell_y,
dwell_directions[:,0], dwell_directions[:,1],
pivot='mid', scale=10
)
ax1.set(aspect=1, xlabel='x (cm)', ylabel='z (cm)')
ax2.set(aspect=1, xlabel='y (cm)', ylabel='z (cm)')
ax3.set(aspect=1, xlabel='x (cm)', ylabel='y (cm)')
return ax4
test_catheter_x = np.arange(0, 3, 0.1)
test_catheter_y = test_catheter_x**2
test_catheter_z = 3*test_catheter_x
test_dwell_x = np.array([0.5, 1, 2, 2.5])
test_dwell_y = test_dwell_x**2
test_dwell_z = 3*test_dwell_x
test_dwell_directions = dwell_direction_from_catheter_coords(
test_catheter_x, test_catheter_y, test_catheter_z,
test_dwell_x, test_dwell_y, test_dwell_z)
show_dwell_directions(
test_catheter_x, test_catheter_y, test_catheter_z,
test_dwell_x, test_dwell_y, test_dwell_z, test_dwell_directions)
<matplotlib.axes._subplots.Axes3DSubplot at 0x11bad7b8>
def test_catheter(catheter_number):
catheter_x, catheter_y, catheter_z = pull_catheter_coords(dcm_plan, catheter_number)
min_z = np.min(dwell_positions[:,2]) - 1
max_z = np.max(dwell_positions[:,2]) + 1
catheter_ref = (catheter_z >= min_z) & (catheter_z <= max_z)
ax = show_dwell_directions(
catheter_x[catheter_ref], catheter_y[catheter_ref], catheter_z[catheter_ref],
dwell_positions[dwell_channels == catheter_number][:,0],
dwell_positions[dwell_channels == catheter_number][:,1],
dwell_positions[dwell_channels == catheter_number][:,2],
dwell_directions[dwell_channels==catheter_number,:])
ax.scatter(
dwell_positions[dwell_channels == catheter_number][:,0],
dwell_positions[dwell_channels == catheter_number][:,1],
dwell_positions[dwell_channels == catheter_number][:,2])
ax.scatter(
dwell_positions[dwell_channels != catheter_number][:,0],
dwell_positions[dwell_channels != catheter_number][:,1],
dwell_positions[dwell_channels != catheter_number][:,2])
ax.plot(
catheter_x[catheter_ref],
catheter_y[catheter_ref],
catheter_z[catheter_ref])
ax.set(xlabel='x', ylabel='y', zlabel='z')
catheter_numbers = np.unique(dwell_channels)
for catheter_number in catheter_numbers:
display(Markdown("#### Catheter Number {}".format(catheter_number)))
test_catheter(catheter_number)
plt.show()
# The x, y, and z defined here have not been sufficiently verified
# They do not necessarily match either what is within Dicom nor what is within your
# TPS. Please verify these and see if they are what you expect them to be.
# If these functions are incorrect or there is a better choice of dimension definitions
# please contact me by creating an issue within the github repository:
# https://github.com/SimonBiggs/npgamma/issues
# If you are able to validate these functions please contact me in the same way.
def load_dose_from_dicom(dcm):
"""Imports the dose in matplotlib format, with the following index mapping:
i = y
j = x
k = z
Therefore when using this function to have the coords match the same order,
ie. coords_reference = (y, x, z)
"""
pixels = np.transpose(
dcm.pixel_array, (1, 2, 0))
dose = pixels * dcm.DoseGridScaling
return dose
def load_xyz_from_dicom(dcm):
"""Although this coordinate pull from Dicom works in the scenarios tested
this is not an official x, y, z pull. It needs further confirmation.
"""
resolution = np.array(
dcm.PixelSpacing).astype(float)
# Does the first index match x?
# Haven't tested with differing grid sizes in x and y directions.
dx = resolution[0]
# The use of dcm.Columns here is under question
x = (
dcm.ImagePositionPatient[0] +
np.arange(0, dcm.Columns * dx, dx))
# Does the second index match y?
# Haven't tested with differing grid sizes in x and y directions.
dy = resolution[1]
# The use of dcm.Rows here is under question
y = (
dcm.ImagePositionPatient[1] +
np.arange(0, dcm.Rows * dy, dy))
# Is this correct?
z = (
np.array(dcm.GridFrameOffsetVector) +
dcm.ImagePositionPatient[2])
return x, y, z
tps_dose = load_dose_from_dicom(dcm_dose)
max_dose = np.max(tps_dose)
# Oncentra
estimated_perscription = max_dose / 10
relevant_dose_deviation = estimated_perscription * 0.03
max_dose_per_slice = np.max(np.max(tps_dose, axis=0), axis=0)
relevant_slice = max_dose_per_slice > estimated_perscription
x_raw, y_raw, z_raw = load_xyz_from_dicom(dcm_dose)
calc_grid_slice_skip = 5
calc_x = x_raw / 10
calc_y = y_raw / 10
calc_z = z_raw[relevant_slice] / 10
calc_z = calc_z[::calc_grid_slice_skip]
calc_grid_positions = determine_calc_grid_positions(calc_x, calc_y, calc_z)
tg43_dose = tg43_on_grid(
calc_x, calc_y, calc_z, reference_air_kerma_rate,
dwell_times, dwell_positions, dwell_directions)
tg43_dose[tg43_dose>max_dose] = max_dose
# Could be axis definition issues here
for z in calc_z:
tps_z_ref = np.where(z_raw/10 == z)[0][0]
calc_z_ref = np.where(calc_z == z)[0][0]
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, sharey=True, figsize=(16,4))
c1 = ax1.contourf(
calc_x, calc_y, tg43_dose[:,:,calc_z_ref], np.arange(0,max_dose/2,0.1), extend='max')
ax1.set_title("Calculated dose @ z = {0:.2f}cm".format(z))
ax1.set_xlabel("x (cm)")
ax1.set_ylabel("y (cm)")
cbar = plt.colorbar(c1, ax=ax1)
cbar.ax.set_ylabel('Dose (Gy)')
c2 = ax2.contourf(
calc_x, calc_y, tps_dose[:,:,tps_z_ref], np.arange(0,max_dose/2,0.01), extend='max')
ax2.set_title("TPS dose @ z = {0:.2f}cm".format(z))
ax2.set_xlabel("x (cm)")
ax2.set_ylabel("y (cm)")
cbar = plt.colorbar(c2, ax=ax2)
cbar.ax.set_ylabel('Dose (Gy)')
dose_diff = tg43_dose[:,:,calc_z_ref] - tps_dose[:,:,tps_z_ref]
c3 = ax3.contourf(
calc_x, calc_y, dose_diff, np.arange(
-relevant_dose_deviation,relevant_dose_deviation,0.001), extend='both', cmap='Spectral_r')
ax3.set_title("Dose difference @ z = {0:.2f}cm".format(z))
ax3.set_xlabel("x (cm)")
ax3.set_ylabel("y (cm)")
cbar = plt.colorbar(c3, ax=ax3)
cbar.ax.set_ylabel('Dose Difference (Gy)')
plt.show()
# z = calc_z[len(calc_z)//2]
# tps_z_ref = np.where(z_raw/10 == z)[0][0]
# calc_z_ref = np.where(calc_z == z)[0][0]
# fig, (ax1, ax2, ax3) = plt.subplots(1, 3, sharey=True, figsize=(16,4))
# c1 = ax1.contourf(
# calc_x, calc_y, tg43_dose[:,:,calc_z_ref], np.arange(0,max_dose/2,0.001), extend='max')
# ax1.set_title("Calculated dose @ z = {0:.2f}cm".format(z))
# ax1.set_xlabel("x (cm)")
# ax1.set_ylabel("y (cm)")
# cbar = plt.colorbar(c1, ax=ax1)
# cbar.ax.set_ylabel('Dose (Gy)')
# dose_diff = tg43_dose[:,:,calc_z_ref] - tps_dose[:,:,tps_z_ref]
# c2 = ax2.contourf(
# calc_x, calc_y, dose_diff, np.arange(-relevant_dose_deviation,relevant_dose_deviation,0.0001), extend='both')
# ax2.set_title("Dose difference @ z = {0:.2f}cm".format(z))
# ax2.set_xlabel("x (cm)")
# ax2.set_ylabel("y (cm)")
# cbar = plt.colorbar(c2, ax=ax2)
# cbar.ax.set_ylabel('Dose Difference (Gy)')
# rel_dose_diff = 100 * dose_diff / tg43_dose[:,:,calc_z_ref]
# c3 = ax3.contourf(
# calc_x, calc_y, rel_dose_diff, np.arange(-3,3,0.0001), extend='both')
# ax3.set_title("Percent relative dose difference @ z = {0:.2f}cm".format(z))
# ax3.set_xlabel("x (cm)")
# ax3.set_ylabel("y (cm)")
# cbar = plt.colorbar(c3, ax=ax3)
# cbar.ax.set_ylabel('Percent Dose Difference')
# plt.show()