Here we provide code and analysis of the kernel shots for all cups done using Google Computing Engine. Total of 12,398 shots were computed, processed and analyzed.
In [1]:
from __future__ import print_function
import sys
import math
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from functools import partial
%matplotlib inline
In [2]:
print('Python version ' + sys.version)
print('Pandas version ' + pd.__version__)
print('Matplotlib version ' + matplotlib.__version__)
First, data shall be prepared by reading all 12,398 shots. Smaller number of shots could be processed as well. We advice to create list of shots to be processed and run script dmax_list.py. All listed shots will be processed and dmax values and (20,50,80) percent of dmax values would be found. For dmax values we average 8 voxels around the shot point, so dose max value is by default, in 2.4x2.4x2.4 mm^3 voxel. X, Y, Z curves to compute (20,50,80) percent values are averaged over 4 voxels. Values are checked for consistency, so that value on curve around the max averaged over two neighbours is exactly equal to dmax.
In [3]:
dt = pd.read_csv("Shots", header=None, delim_whitespace=True, names=["Cup", "Num", "Col", "ShY", "ShZ", "Dm", "X20L", "X20R", "X50L", "X50R", "X80L", "X80R", "Y20L", "Y20R", "Y50L", "Y50R", "Y80L", "Y80R", "Z20L", "Z20R", "Z50L", "Z50R", "Z80L", "Z80R"])
dt = dt.sort(["Cup", "Num", "Col", "ShY", "ShZ"], ascending=[0, 1, 0, 1, 1])
dt.index = range(0, len(dt)) # make it reindex
Print first five rows of the data
In [4]:
dt.head() # first five rows
Out[4]:
Now last five rows of the data
In [5]:
dt.tail() # last five rows
Out[5]:
In [6]:
dt.dtypes # types of the columns
Out[6]:
In [7]:
dt.info() # some summary
So far looks right, we have first three classification columns (with Cup name, Cup number and collimator), then two columns with shot Y,Z position and Dmax. Last 18 columns are six per axis, with two values (L marking Left and R marking Right) per one percentage value out of (20, 50, 80).
We will plot all cups with dmax placed at particular shot location. Per page we make two plots, left one being actual cup with its size and number. Right plot will be for a largest cup of the series. We expect dose ponts first, to fill the cup shape, and second, to be consistent with general pattern of being maximum at close to inner cup are and a bit smaller when going to central cup area.
In [8]:
def find_idx_nearby_shot(df, sh_y, sh_z):
"""
Given shot Y and shot Z,
find row index for this shot and
"""
cY = df["ShY"]
cZ = df["ShZ"]
for idx in df.index:
if math.fabs( cY.loc[idx] - sh_y ) < 0.5 and math.fabs( cZ.loc[idx] - sh_z ) < 0.5:
return idx
return -1
def make_an_image(df, imfun):
"""
Make an image out of data frame
"""
ymin = df["ShY"].min()
ymax = df["ShY"].max()
zmin = df["ShZ"].min()
zmax = df["ShZ"].max()
#print(ymin, ymax, zmin, zmax)
step = 5.0
ny = int( np.around((ymax - ymin)/step) ) + 1
nz = int( np.around((zmax - zmin)/step) ) + 1
#print(ny, nz)
sh_dm = np.empty((nz, ny))
for iz in range(0, nz):
z = zmin + float(iz) * step
for iy in range(0, ny):
y = ymin + float(iy) * step
dm = imfun(df, y, z)
sh_dm[iz, iy] = dm
return sh_dm
def max_image(dt, cup, col, imfun):
q = (dt["Cup"] == cup) & (dt["Col"] == col)
cups = dt[q]
nums = cups["Num"]
cmax = nums.max()
dc = cups.loc[cups["Num"] == cmax]
return (cmax, make_an_image(dc, imfun))
In [9]:
def dmax_at_nearby_shot(df, sh_y, sh_z):
"""
Given shot Y and shot Z,
find row index for this shot and
return dmax
"""
idx = find_idx_nearby_shot(df, sh_y, sh_z)
if idx < 0:
return 0.0
return df["Dm"].loc[idx]
max_imgs = {}
max_imgs["L25"] = max_image(dt, "L", 25, dmax_at_nearby_shot)
max_imgs["M25"] = max_image(dt, "M", 25, dmax_at_nearby_shot)
max_imgs["S25"] = max_image(dt, "S", 25, dmax_at_nearby_shot)
max_imgs["L15"] = max_image(dt, "L", 15, dmax_at_nearby_shot)
max_imgs["M15"] = max_image(dt, "M", 15, dmax_at_nearby_shot)
max_imgs["S15"] = max_image(dt, "S", 15, dmax_at_nearby_shot)
kk = 0
for cn, dc in dt.groupby(['Cup', 'Num', 'Col']):
cup, num, col = cn
sh_dm = make_an_image(dc, dmax_at_nearby_shot)
img = None
fig, axes = plt.subplots(1, 2, figsize=(12, 7), subplot_kw={'xticks': [], 'yticks': []})
left = True
ncup = num
for ax in axes.flat:
if left:
img = ax.imshow(sh_dm, interpolation="none")
left = False
else:
img = ax.imshow(max_imgs[cup+str(col)][1], interpolation="none")
ncup = max_imgs[cup+str(col)][0]
title = str(cup) + "_"
if ncup < 10:
title += "0"
title += str(ncup) + " " + str(col) + "mm"
ax.set_title(title)
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.25, 0.05, 0.5])
fig.colorbar(img, cax=cbar_ax)
plt.show()
kk += 1
We define gates as left and right values at the same dose percentage level. Code below is checking if gates are always consistent, i.e. left gate value is always smaller than right gate value. We expect all gates to be consistent, i.e. all left values to be smaller than right values.
In [10]:
for cn, dc in dt.groupby(['Cup', 'Num', 'Col']):
cup, num, col = cn
q = dc["X20L"] < dc["X20R"]
good = True
if len(q) > q.sum():
print("Inconsistent X-20% gates for cup {0}, num {1}".format(cup, num))
good = False
if good:
print("Consistent X-20% gates for all cups".format(cup, num))
In [11]:
for cn, dc in dt.groupby(['Cup', 'Num', 'Col']):
cup, num, col = cn
q = dc["X50L"] < dc["X50R"]
good = True
if len(q) > q.sum():
print("Inconsistent X-50% gates for cup {0}, num {1}".format(cup, num))
good = False
if good:
print("Consistent X-50% gates for all cups".format(cup, num))
In [12]:
for cn, dc in dt.groupby(['Cup', 'Num', 'Col']):
cup, num, col = cn
q = dc["X80L"] < dc["X80R"]
good = True
if len(q) > q.sum():
print("Inconsistent X-80% gates for cup {0}, num {1}".format(cup, num))
good = False
if good:
print("Consistent X-80% gates for all cups".format(cup, num))
In [13]:
for cn, dc in dt.groupby(['Cup', 'Num', 'Col']):
cup, num, col = cn
q = dc["Y20L"] < dc["Y20R"]
good = True
if len(q) > q.sum():
print("Inconsistent Y-20% gates for cup {0}, num {1}".format(cup, num))
good = False
if good:
print("Consistent Y-20% gates for all cups".format(cup, num))
In [14]:
for cn, dc in dt.groupby(['Cup', 'Num', 'Col']):
cup, num, col = cn
q = dc["Y50L"] < dc["Y50R"]
good = True
if len(q) > q.sum():
print("Inconsistent Y-50% gates for cup {0}, num {1}".format(cup, num))
good = False
if good:
print("Consistent Y-50% gates for all cups".format(cup, num))
In [15]:
for cn, dc in dt.groupby(['Cup', 'Num', 'Col']):
cup, num, col = cn
q = dc["Y80L"] < dc["Y80R"]
good = True
if len(q) > q.sum():
print("Inconsistent Y-80% gates for cup {0}, num {1}".format(cup, num))
good = False
if good:
print("Consistent Y-80% gates for all cups".format(cup, num))
In [16]:
for cn, dc in dt.groupby(['Cup', 'Num', 'Col']):
cup, num, col = cn
q = dc["Z20L"] < dc["Z20R"]
good = True
if len(q) > q.sum():
print("Inconsistent Z-20% gates for cup {0}, num {1}".format(cup, num))
good = False
if good:
print("Consistent Z-20% gates for all cups".format(cup, num))
In [17]:
for cn, dc in dt.groupby(['Cup', 'Num', 'Col']):
cup, num, col = cn
q = dc["Z50L"] < dc["Z50R"]
good = True
if len(q) > q.sum():
print("Inconsistent Z-50% gates for cup {0}, num {1}".format(cup, num))
good = False
if good:
print("Consistent Z-50% gates for all cups".format(cup, num))
In [18]:
for cn, dc in dt.groupby(['Cup', 'Num', 'Col']):
cup, num, col = cn
q = dc["Z80L"] < dc["Z80R"]
good = True
if len(q) > q.sum():
print("Inconsistent Z-80% gates for cup {0}, num {1}".format(cup, num))
good = False
if good:
print("Consistent Z-80% gates for all cups".format(cup, num))
This check is for data to be monotonic on each side, so that on the left, for example, 20% gate will be larger than and contains 50% gate, which, in turn, will be larger than and contains 80% gate. We expect such test will pass for all directions.
In [19]:
for cn, dc in dt.groupby(['Cup', 'Num', 'Col']):
cup, num, col = cn
q = (dc["X20L"] < dc["X50L"]) & (dc["X20R"] > dc["X50R"])
good = True
if len(q) > q.sum():
print("Inconsistent X-20-50% gates for cup {0}, num {1}".format(cup, num))
good = False
if good:
print("Consistent X-20-50% gates for all cups".format(cup, num))
In [20]:
for cn, dc in dt.groupby(['Cup', 'Num', 'Col']):
cup, num, col = cn
q = (dc["X50L"] < dc["X80L"]) & (dc["X50R"] > dc["X80R"])
good = True
if len(q) > q.sum():
print("Inconsistent X-50-80% gates for cup {0}, num {1}".format(cup, num))
good = False
if good:
print("Consistent X-50-80% gates for all cups".format(cup, num))
In [21]:
for cn, dc in dt.groupby(['Cup', 'Num', 'Col']):
cup, num, col = cn
q = (dc["Y20L"] < dc["Y50L"]) & (dc["Y20R"] > dc["Y50R"])
good = True
if len(q) > q.sum():
print("Inconsistent Y-20-50% gates for cup {0}, num {1}".format(cup, num))
good = False
if good:
print("Consistent Y-20-50% gates for all cups".format(cup, num))
In [22]:
for cn, dc in dt.groupby(['Cup', 'Num', 'Col']):
cup, num, col = cn
q = (dc["Y50L"] < dc["Y80L"]) & (dc["Y50R"] > dc["Y80R"])
good = True
if len(q) > q.sum():
print("Inconsistent Y-50-80% gates for cup {0}, num {1}".format(cup, num))
good = False
if good:
print("Consistent Y-50-80% gates for all cups".format(cup, num))
In [23]:
for cn, dc in dt.groupby(['Cup', 'Num', 'Col']):
cup, num, col = cn
q = (dc["Z20L"] < dc["Z50L"]) & (dc["Z20R"] > dc["Z50R"])
good = True
if len(q) > q.sum():
print("Inconsistent Z-20-50% gates for cup {0}, num {1}".format(cup, num))
good = False
if good:
print("Consistent Z-20-50% gates for all cups".format(cup, num))
In [24]:
for cn, dc in dt.groupby(['Cup', 'Num', 'Col']):
cup, num, col = cn
q = (dc["Z50L"] < dc["Z80L"]) & (dc["Z50R"] > dc["Z80R"])
good = True
if len(q) > q.sum():
print("Inconsistent Z-50-80% gates for cup {0}, num {1}".format(cup, num))
good = False
if good:
print("Consistent Z-50-80% gates for all cups".format(cup, num))
We are plotting centers of the the gates for each dose level value of (20, 50, 80)% of dmax. We expect computed centers to form good pattern similar to the shot centers set for a compuation (steps of 5mm in Y and Z directions). Shots on the borders are expected to be in less regular pattern and have some anomalies.
In [25]:
# compute centers
dt["X20C"] = 0.5*(dt.X20L + dt.X20R)
dt["X50C"] = 0.5*(dt.X50L + dt.X50R)
dt["X80C"] = 0.5*(dt.X80L + dt.X80R)
dt["Y20C"] = 0.5*(dt.Y20L + dt.Y20R)
dt["Y50C"] = 0.5*(dt.Y50L + dt.Y50R)
dt["Y80C"] = 0.5*(dt.Y80L + dt.Y80R)
dt["Z20C"] = 0.5*(dt.Z20L + dt.Z20R)
dt["Z50C"] = 0.5*(dt.Z50L + dt.Z50R)
dt["Z80C"] = 0.5*(dt.Z80L + dt.Z80R)
dt.info()
In [26]:
for cn, dc in dt.groupby(['Cup', 'Num', 'Col']):
cup, num, col = cn
fig, axes = plt.subplots(1, 3, figsize=(15, 8))
k = 0
ymax = max(dc.Z20C.max(), dc.Z50C.max(), dc.Z80C.max())
ymax = float(int(ymax)/10 + 1 ) * 10.0
for ax in axes.flat:
if k == 0:
x = dc.Y20C
y = dc.Z20C
title = "20% "
if k == 1:
x = dc.Y50C
y = dc.Z50C
title = "50% "
if k == 2:
x = dc.Y80C
y = dc.Z80C
title = "80% "
title += str(cup) + "_"
if num < 10:
title += "0"
title += str(num) + " " + str(col) + "mm"
ax.set_title(title)
ax.set_ylim([-5,ymax])
ax.set_xlabel('Y, mm')
ax.set_ylabel('Z, mm')
ax.scatter(x, y)
k += 1
#fig.subplots_adjust(right=0.8)
plt.show()
We plot X centers here, first Z vs X, then Y vs X. We expect plots to be centered about 0mm exactly.
In [27]:
for cn, dc in dt.groupby(['Cup', 'Num', 'Col']):
cup, num, col = cn
fig, axes = plt.subplots(1, 3, figsize=(15, 8))
k = 0
ymax = max(dc.Z20C.max(), dc.Z50C.max(), dc.Z80C.max())
ymax = float(int(ymax)/10 + 1 ) * 10.0
for ax in axes.flat:
if k == 0:
x = dc.X20C
y = dc.Z20C
title = "20% "
if k == 1:
x = dc.X50C
y = dc.Z50C
title = "50% "
if k == 2:
x = dc.X80C
y = dc.Z80C
title = "80% "
title += str(cup) + "_"
if num < 10:
title += "0"
title += str(num) + " " + str(col) + "mm"
ax.set_title(title)
ax.set_ylim([-5,ymax])
ax.set_xlabel('X, mm')
ax.set_ylabel('Z, mm')
ax.scatter(x, y)
k += 1
plt.show()
We plot X vs Y. We expect plots to be centered about 0 exactly.
In [28]:
for cn, dc in dt.groupby(['Cup', 'Num', 'Col']):
cup, num, col = cn
fig, axes = plt.subplots(1, 3, figsize=(15, 8))
k = 0
ymax = max(dc.Y20C.max(), dc.Y50C.max(), dc.Y80C.max())
ymax = float(int(ymax)/10 + 1 ) * 10.0
for ax in axes.flat:
if k == 0:
x = dc.X20C
y = dc.Y20C
title = "20% "
if k == 1:
x = dc.X50C
y = dc.Y50C
title = "50% "
if k == 2:
x = dc.X80C
y = dc.Y80C
title = "80% "
title += str(cup) + "_"
if num < 10:
title += "0"
title += str(num) + " " + str(col) + "mm"
ax.set_title(title)
ax.set_ylim([-5,ymax])
ax.set_xlabel('X, mm')
ax.set_ylabel('Y, mm')
ax.scatter(x, y)
k += 1
plt.show()
We're plotting the width of the gates for each cup for each value of (20,50,80) percent of Dmax. Gate value is plotted at exact shot position. We expect gates to have somewhat constant values in the bulk of the cup, but dropping out and down to zero when on the border and outside the cup.
In [29]:
def delta_at_nearby_shot(df, sh_y, sh_z, R, L):
"""
Given shot Y and shot Z,
find row index for this shot and
return dmax
"""
idx = find_idx_nearby_shot(df, sh_y, sh_z)
if idx < 0:
return 0.0
return df[R].loc[idx] - df[L].loc[idx]
def draw_delta_images(what, deltaXX_at_nearby_shot):
"""
"""
max_imgs = {}
max_imgs["L25"] = max_image(dt, "L", 25, deltaXX_at_nearby_shot)
max_imgs["M25"] = max_image(dt, "M", 25, deltaXX_at_nearby_shot)
max_imgs["S25"] = max_image(dt, "S", 25, deltaXX_at_nearby_shot)
max_imgs["L15"] = max_image(dt, "L", 15, deltaXX_at_nearby_shot)
max_imgs["M15"] = max_image(dt, "M", 15, deltaXX_at_nearby_shot)
max_imgs["S15"] = max_image(dt, "S", 15, deltaXX_at_nearby_shot)
for cn, dc in dt.groupby(['Cup', 'Num', 'Col']):
cup, num, col = cn
sh_dm = make_an_image(dc, deltaXX_at_nearby_shot)
img = None
fig, axes = plt.subplots(1, 2, figsize=(12, 7), subplot_kw={'xticks': [], 'yticks': []})
left = True
ncup = num
for ax in axes.flat:
if left:
img = ax.imshow(sh_dm, interpolation="none")
left = False
else:
img = ax.imshow(max_imgs[cup+str(col)][1], interpolation="none")
ncup = max_imgs[cup+str(col)][0]
title = what + ", " + str(cup) + "_"
if ncup < 10:
title += "0"
title += str(ncup) + " " + str(col) + "mm"
ax.set_title(title)
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.25, 0.05, 0.5])
fig.colorbar(img, cax=cbar_ax)
plt.show()
draw_delta_images("Y 20% delta", partial(delta_at_nearby_shot, R = "Y20R", L = "Y20L"))
In [30]:
draw_delta_images("Y 50% delta", partial(delta_at_nearby_shot, R = "Y50R", L = "Y50L"))
In [31]:
draw_delta_images("Y 80% delta", partial(delta_at_nearby_shot, R = "Y80R", L = "Y80L"))
In [32]:
draw_delta_images("Z 20% delta", partial(delta_at_nearby_shot, R = "Z20R", L = "Z20L"))
In [33]:
draw_delta_images("Z 50% delta", partial(delta_at_nearby_shot, R = "Z50R", L = "Z50L"))
In [34]:
draw_delta_images("Z 80% delta", partial(delta_at_nearby_shot, R = "Z80R", L = "Z80L"))
In [35]:
draw_delta_images("X 20% delta", partial(delta_at_nearby_shot, R = "X20R", L = "X20L"))
In [36]:
draw_delta_images("X 50% delta", partial(delta_at_nearby_shot, R = "X50R", L = "X50L"))
In [37]:
draw_delta_images("X 80% delta", partial(delta_at_nearby_shot, R = "X80R", L = "X80L"))
Here we will plot value of Dmax at (0,0) versus cup size. We expect such values to be quite constant with small decrease toward larger cups. When we compare S cups vs M cups vs L cups, we expect S cups to have larger doses at (0,0) than M, and, in turn, M cups to have larger doses than L cups.
In [38]:
for cc, dc in dt.groupby(['Cup', 'Col']):
cup, col = cc
fig = plt.figure()
ax = fig.add_subplot(111)
## the data
q = (dc.ShY == 0) & (dc.ShZ == 0)
t = dc[q]
x = t.Num
y = t.Dm
title = "Plotting Dmax vs Cup size for " + str(cup) + " " + str(col) + "mm"
ax.set_title(title)
ax.set_ylim([0.0, 1.2*dc.Dm.max()])
ax.set_xlabel('Cup number')
ax.set_ylabel('Dose, Gy/photon')
ax.scatter(x, y)
plt.show()
We examined doses and profiles/gates using Pandas and Matplotlib. A lot of plots were drawn, and checks where done. After examining all available data, I believe dose computation for all 12.398 was consistent with input and expectations. Profiles/gates were consistent with expected collimator spot size, and dose values were quite reasonable. No outliers were discovered so far. This conclude non-measurement related part of kernels validation and verification procedure.
In [ ]: