In [33]:
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.interp_fluences = list()
#TODO more rigorous checking of my_dir
if not os.path.exists(my_dir):
self.my_logs.load_dir_UI()
else:
self.my_logs.load_dir(my_dir)
def do_calcs(self, res_mm=1.0):
if len(self.my_logs) == 0:
print("No valid log files loaded ...")
return
for log in self.my_logs:
log.fluence.actual.calc_map(resolution=res_mm)
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(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
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 [34]:
my_dyn2dose = Dyn_to_Dose('/home/mpre/tmp_dynalogs_la4/OBRIEN_Michael')
In [35]:
my_dyn2dose.do_calcs(1.0)
In [36]:
my_dyn2dose.do_flip()
In [37]:
pin_x = np.linspace(-12.7, 12.7, 128)
pin_y = np.linspace(-9.1, 9.1, 92)
my_dyn2dose.make_interp(pin_y, pin_x)
In [38]:
import matplotlib.pyplot as plt
%matplotlib inline
fig = plt.figure(figsize=(6, 3.2))
ax1 = fig.add_subplot(111)
plt.imshow(my_dyn2dose.interp_fluences[2])
Out[38]:
In [32]:
pin_y
Out[32]:
In [45]:
import PinDoseMap
In [47]:
my_dose_map = PinDoseMap.PinDoseMap('/home/mpre/Documents/10')
In [48]:
my_dyn2dose.make_interp(my_dose_map.y_pos, my_dose_map.x_pos)
In [ ]:
my_dose_map.do_plot()
In [50]:
%matplotlib inline
fig = plt.figure(figsize=(6, 3.2))
ax1 = fig.add_subplot(111)
plt.imshow(my_dyn2dose.interp_fluences[2])
Out[50]:
In [51]:
my_flu_sample = my_dyn2dose.interp_fluences[2][25,31]
In [52]:
my_flu_sample
Out[52]:
In [55]:
my_dose_sample = my_dose_map.itrp_dose_array[25,31]
my_dose_sample
Out[55]:
In [56]:
my_norm = 0.81899999999999973 / 0.66608057743943438
In [57]:
my_norm
Out[57]:
In [62]:
new_flu = my_norm * my_dyn2dose.interp_fluences[2]
In [63]:
my_diff = new_flu - my_dose_map.itrp_dose_array
In [64]:
fig = plt.figure(figsize=(6, 3.2))
ax2 = fig.add_subplot(111)
In [65]:
plt.imshow(my_diff)plt.
Out[65]:
In [73]:
my_dyn2dose.my_logs[2]
Out[73]:
In [75]:
def print_attributes(obj):
for attr in obj.__dict__:
print (attr, getattr(obj, attr))
print_attributes(my_dyn2dose.my_logs[2])
In [78]:
print_attributes(my_dyn2dose.my_logs[2].axis_data.jaws)
In [91]:
print( my_dyn2dose.my_logs[2].axis_data.jaws.y2.actual)
In [95]:
print_attributes(my_dyn2dose.my_logs[2].fluence.actual._jaws.x2)
In [109]:
rows = np.linspace(-2.5, 3.5, 61)
In [110]:
rows
Out[110]:
In [96]:
cols = np.linspace(-7., 1., 81)
In [97]:
cols
Out[97]:
In [111]:
my_dyn2dose.make_interp(rows, cols)
In [112]:
plt.imshow(my_dyn2dose.interp_fluences[2])
Out[112]:
In [102]:
interpolate
In [105]:
my_tmp = interpl.RectBivariateSpline(my_dose_map.y_pos, my_dose_map.x_pos, my_dose_map.raw_dose_array)
In [113]:
my_pin_zoom = my_tmp(rows, cols)
In [114]:
plt.imshow(my_pin_zoom)
Out[114]:
In [115]:
def gauss2D(shape=(3,3),sigma=0.5):
"""
2D gaussian mask - should give the same result as MATLAB's
fspecial('gaussian',[shape],[sigma])
"""
m,n = [(ss-1.)/2. for ss in shape]
y,x = np.ogrid[-m:m+1,-n:n+1]
h = np.exp( -(x*x + y*y) / (2.*sigma*sigma) )
h[ h < np.finfo(h.dtype).eps*h.max() ] = 0
sumh = h.sum()
if sumh != 0:
h /= sumh
return h
In [116]:
from scipy import ndimage as ndi
In [117]:
kern1 = gauss2D((3,3), 0.5)
In [118]:
blur1 = ndi.convolve(my_dyn2dose.interp_fluences[2], kern1, mode='constant', cval=0.0)
In [119]:
plt.imshow(blur1)
Out[119]:
In [120]:
kern1
Out[120]:
In [183]:
kern2 = gauss2D((5,5),1)
blur2 = ndi.convolve(my_dyn2dose.interp_fluences[2], kern2, mode='reflect')
blur3 = ndi.filters.gaussian_filter(my_dyn2dose.interp_fluences[2], 3.0, mode='reflect')
plt.imshow(blur3)
Out[183]:
In [184]:
my_diff1 = blur3 - my_pin_zoom
In [185]:
plt.imshow(my_diff1)
Out[185]:
In [ ]: