In [3]:
from scipy import interpolate as interpl
import numpy as np
import os
from pylinac import log_analyzer as lga
class Dyn_to_Dose:
def __init__(self,my_dir):
"""class to import a folder full of dynalog files for a treatment"""
self.my_logs = lga.MachineLogs();
self.flip_fluences = list()
self.gant_angle = list()
self.interp_fluences = list()
self.res_mm = 1.0
self.num_logs = None
self.pin_dose_index = list()
#TODO more rigorous checking of my_dir
if not os.path.exists(my_dir):
self.my_logs.from_folder_UI()
else:
self.my_logs.load_folder(my_dir)
self.num_logs = len(self.my_logs)
self.interp_fluences = [None] * self.num_logs
self.pin_dose_index = [None] * self.num_logs
def do_calcs(self, res_mm=1.0):
if len(self.my_logs) == 0:
print("No valid log files loaded ...")
return
self.gant_angle = list()
for log in self.my_logs:
log.fluence.actual.calc_map(resolution=res_mm)
self.gant_angle.append((540.0 - log.axis_data.gantry.actual[0])%360.0)
def do_flip(self):
for log in self.my_logs:
tmp = log.fluence.actual.calc_map(resolution=1.0)
self.flip_fluences.append(np.flipud(tmp))
def make_interp_single(self, index, y_dim, x_dim):
assert index <= self.num_logs, "index out of bounds"
flu_y = [-19.5, -18.5, -17.5, -16.5, -15.5, -14.5, -13.5, -12.5, -11.5, -10.5, -9.75,
-9.25, -8.75, -8.25, -7.75, -7.25, -6.75, -6.25, -5.75, -5.25, -4.75, -4.25,
-3.75, -3.25, -2.75, -2.25, -1.75, -1.25, -0.75, -0.25, 0.25, 0.75, 1.25,
1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75, 5.25, 5.75, 6.25, 6.75, 7.25, 7.75,
8.25, 8.75, 9.25, 9.75, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5,
19.5]
flu_x = np.linspace(-19.95, 19.95, 400) #TODO get rid of hardcoding dimensions of 1mm res
my_flu = self.flip_fluences[index]
#first create interpolating object tmp
tmp = interpl.RectBivariateSpline(flu_y, flu_x, my_flu)
# now create new interpolated fluences and store
self.interp_fluences[index] = tmp(y_dim, x_dim)
def make_interp(self,y_dim, x_dim):
"""y_dim must be 1D array containing co-ords in cms
for fluence and x_dim 1D array likewise. x is parallel to leaf motion
"""
flu_y = [-19.5, -18.5, -17.5, -16.5, -15.5, -14.5, -13.5, -12.5, -11.5, -10.5, -9.75,
-9.25, -8.75, -8.25, -7.75, -7.25, -6.75, -6.25, -5.75, -5.25, -4.75, -4.25,
-3.75, -3.25, -2.75, -2.25, -1.75, -1.25, -0.75, -0.25, 0.25, 0.75, 1.25,
1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75, 5.25, 5.75, 6.25, 6.75, 7.25, 7.75,
8.25, 8.75, 9.25, 9.75, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5,
19.5]
flu_x = np.linspace(-19.95, 19.95, 400) #TODO get rid of hardcoding dimensions of 1mm res
self.interp_fluences = list() # need to empty out list first
for flu in self.flip_fluences:
#first create interpolating object tmp
tmp = interpl.RectBivariateSpline(flu_y, flu_x, flu)
# now create new interpolated fluences and store
self.interp_fluences.append(tmp(y_dim, x_dim))
In [4]:
my_dyn2dose = Dyn_to_Dose('/home/mpre/tmp_dynalogs_la4/FORD')
In [5]:
my_dyn2dose.do_calcs(1.0)
In [6]:
my_dyn2dose.do_flip()
In [7]:
import PinDoseMap
import os
os.chdir('/home/mpre/Documents/PlaDosTool')
import PinDoseMapList as pdml
In [8]:
my_dose_maps = pdml.PinDoseMapList('/home/mpre/tmp_Pinn_dose_maps/FORD_20150611')
In [9]:
for x in my_dose_maps.pdm_list:
print(x.gant_angle)
In [10]:
import matplotlib.pyplot as plt
%matplotlib inline
#ax1 = fig.add_subplot(911)
flu_y = [-19.5, -18.5, -17.5, -16.5, -15.5, -14.5, -13.5, -12.5, -11.5, -10.5, -9.75,
-9.25, -8.75, -8.25, -7.75, -7.25, -6.75, -6.25, -5.75, -5.25, -4.75, -4.25,
-3.75, -3.25, -2.75, -2.25, -1.75, -1.25, -0.75, -0.25, 0.25, 0.75, 1.25,
1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75, 5.25, 5.75, 6.25, 6.75, 7.25, 7.75,
8.25, 8.75, 9.25, 9.75, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5,
19.5]
flu_x = np.linspace(-19.95, 19.95, 400)
for i in range (0,len(my_dose_maps.pdm_list)):
j=0
for ang in my_dyn2dose.gant_angle:
#print(int(ang))
#print (int(my_dose_maps.pdm_list[i].gant_angle))
if (ang == int(my_dose_maps.pdm_list[i].gant_angle)):
print('j is {0:d}'.format(j))
my_x = my_dose_maps.pdm_list[i].x_pos
my_y = my_dose_maps.pdm_list[i].y_pos
# NB - use j, not i here
flu1 = my_dyn2dose.flip_fluences[j]
tmp = interpl.RectBivariateSpline(flu_y, flu_x, flu1)
blah = tmp(my_y, my_x)
my_dyn2dose.interp_fluences[j] = blah
my_dyn2dose.pin_dose_index[j] = i
# FIXME have to sort out negative x jaws
fig = plt.figure(i+1)
#plt.subplot(9,1,(i+1))
plt.imshow(blah)
j+=1
In [11]:
my_dose_maps.pdm_list[0].do_plot()
In [12]:
print(my_dose_maps.pdm_list[0].y_pos)
print(my_dose_maps.pdm_list[0].x_pos)
In [13]:
my_diff = my_dyn2dose.interp_fluences[2] - my_dose_maps.pdm_list[0].itrp_dose_array
In [14]:
plt.imshow(my_diff)
Out[14]:
In [15]:
my_dyn2dose.my_logs[2]
Out[15]:
In [16]:
def print_attributes(obj):
for attr in obj.__dict__:
print (attr, getattr(obj, attr))
print_attributes(my_dyn2dose.my_logs[2])
In [17]:
print_attributes(my_dyn2dose.my_logs[2].axis_data.jaws)
In [18]:
print( my_dyn2dose.my_logs[2].axis_data.jaws.y2.actual)
In [19]:
print(my_dyn2dose.my_logs[2].axis_data.jaws.x2.actual)
In [20]:
print(my_dyn2dose.my_logs[2].axis_data.jaws.y1.actual)
In [21]:
print(my_dyn2dose.my_logs[2].axis_data.jaws.x1.actual)
In [22]:
my_pin_zoom_list = [None]*my_dyn2dose.num_logs
for i in range (0,my_dyn2dose.num_logs):
print(i)
y2 = (my_dyn2dose.my_logs[i].axis_data.jaws.y2.actual)[0]
x2 = (my_dyn2dose.my_logs[i].axis_data.jaws.x2.actual)[0]
y1 = (my_dyn2dose.my_logs[i].axis_data.jaws.y1.actual)[0]
x1 = (my_dyn2dose.my_logs[i].axis_data.jaws.x1.actual)[0]
print(y2)
print(y1)
print(x1)
print(x2)
print ("=====", '\r\n')
rows = np.linspace(-1.0*(y2-1.0), (y1 - 1.0), ((y2+y1)*10 +1-20))
cols = np.linspace(-1.0*(x1), x2, (x1+x2)*10 +1)
#print(rows)
#print (cols)
dm = my_dose_maps.pdm_list[my_dyn2dose.pin_dose_index[i]]
my_tmp = interpl.RectBivariateSpline(dm.y_pos, dm.x_pos, dm.raw_dose_array)
my_pin_zoom_list[i] = my_tmp(rows, cols)
my_dyn2dose.make_interp_single(i, rows, cols)
In [23]:
for i in range (0, my_dyn2dose.num_logs):
fig = plt.figure(i*2+1)
plt.imshow(my_dyn2dose.interp_fluences[i])
fig = plt.figure(i*2+2)
plt.imshow(my_pin_zoom_list[i])
In [24]:
for log in my_dyn2dose.my_logs:
log.fluence.actual.calc_map(resolution=0.1)
log.fluence.actual.plot_map()
In [ ]:
In [25]:
my_dyn2dose.make_interp_single(2, rows, cols)
In [26]:
plt.imshow(my_dyn2dose.interp_fluences[2])
Out[26]:
In [27]:
my_dose_map = my_dose_maps.pdm_list[0]
In [28]:
my_tmp = interpl.RectBivariateSpline(my_dose_map.y_pos, my_dose_map.x_pos, my_dose_map.raw_dose_array)
In [29]:
my_pin_zoom = my_tmp(rows, cols)
In [147]:
blur3 = ndi.filters.gaussian_filter(my_dyn2dose.interp_fluences[2],(3.1)
, mode='reflect')
plt.imshow(blur3)
Out[147]:
In [120]:
pin_max = np.max(my_pin_zoom)
dyn_max = np.max(blur3)
print(pin_max, " ", dyn_max)
In [132]:
from scipy import ndimage as ndi
sigmas = list()
IQRs = list()
for l in range(20, 59, ):
blur3 = ndi.filters.gaussian_filter(my_dyn2dose.interp_fluences[2],(0.1*l)
, mode='reflect')
my_diff1 = 100.0*abs(blur3 - my_pin_zoom)
my_grad = np.gradient(my_diff1)
#result = np.percentile(my_diff1, [25,50,75])
#print('iteration {0:d}, Q75 {1}, Q50 {2}, Q25 {2}'.format(l,result[2],result[1], result[0]))
Q25 = np.percentile(my_grad, 25)
m = np.mean(my_grad)
Q75 = np.percentile(my_grad, 75)
res2 = np.mean(my_grad)
IQR = Q75 - Q25
#print("Q25 is {0:8.4f} and MEAN is {1:8.4f}, Q75 {2:8.4f}\t sigma is {3:3.1f} \n IQR is {4:6.4f}"
# .format(Q25, m, Q75, (0.1*l), IQR))
print("Sigma: {0:3.1f}\t IQR: {1:6.4f}".format((0.1*l), IQR))
sigmas.append(0.1*l)
IQRs.append(IQR)
blah = plt.scatter(np.asarray(sigmas), np.asarray(IQRs))
plt.ylabel('differential metric')
plt.xlabel('Sigma')
plt.show()
#plt.imshow(blur3)
In [133]:
blur3 = ndi.filters.gaussian_filter(my_dyn2dose.interp_fluences[2],(3.1)
, mode='reflect')
my_diff1 = abs(blur3 - my_pin_zoom)
In [134]:
from matplotlib import cm
fig, ax = plt.subplots()
data = 100*(my_diff1/np.max(my_diff1))
im = ax.imshow(data)
cax = fig.add_axes([0.9, 0.1, 0.03, 0.8])
cbar = fig.colorbar(im, cax=cax, ticks=[0.0, 50.0, 100.0])
#cbar.ax.set_yticklabels(['0.0', '0.1', '0.2'])
plt.show()
#cax = ax.imshow(data, interpolation='nearest', cmap=cm.coolwarm)
#ax.set_title('Diffs with vertical colorbar')
print (100*np.max(my_diff1))
#plt.imshow(my_diff1)
In [135]:
q75, q25 = np.percentile(my_diff1, [75 ,25])
print("quartile 25% is {0} and quartile 75% is {1}".format(q25, q75))
print("interquartile range is {0}".format(1000*(q75-q25)))
In [43]:
q75
Out[43]:
In [44]:
np.max(my_diff1)
Out[44]:
In [136]:
matplotlib qt
In [137]:
from matplotlib import cm
fig, ax = plt.subplots()
data = 100*(my_diff1/np.max(my_diff1))
im = ax.imshow(data)
cax = fig.add_axes([0.9, 0.1, 0.03, 0.8])
cbar = fig.colorbar(im, cax=cax, ticks=[0.0, 50.0, 100.0])
#cbar.ax.set_yticklabels(['0.0', '0.1', '0.2'])
plt.show()
#cax = ax.imshow(data, interpolation='nearest', cmap=cm.coolwarm)
#ax.set_title('Diffs with vertical colorbar')
print (100*np.max(my_diff1))
#plt.imshow(my_diff1)
In [142]:
from scipy import ndimage as ndi
sigmas = list()
IQRs = list()
for l in range(20, 59, ):
blur3 = ndi.filters.gaussian_filter(my_dyn2dose.interp_fluences[2],(0.1*l)
, mode='reflect')
my_diff1 = 100.0*abs(blur3 - my_pin_zoom)
my_grad = np.gradient(my_diff1)
#result = np.percentile(my_diff1, [25,50,75])
#print('iteration {0:d}, Q75 {1}, Q50 {2}, Q25 {2}'.format(l,result[2],result[1], result[0]))
Q25 = np.percentile(my_grad, 25)
m = np.mean(my_grad)
Q75 = np.percentile(my_grad, 75)
res2 = np.mean(my_grad)
IQR = Q75 - Q25
#print("Q25 is {0:8.4f} and MEAN is {1:8.4f}, Q75 {2:8.4f}\t sigma is {3:3.1f} \n IQR is {4:6.4f}"
# .format(Q25, m, Q75, (0.1*l), IQR))
print("Sigma: {0:3.1f}\t IQR: {1:6.4f}".format((0.1*l), IQR))
sigmas.append(0.1*l)
IQRs.append(IQR)
blah = plt.scatter(np.asarray(sigmas), np.asarray(IQRs))
plt.ylabel('differential metric')
plt.xlabel('Sigma')
plt.show()
#plt.imshow(blur3)
In [140]:
In [ ]: