Proof of concept

This is a proof of concept for the inclusion of positional uncertainty within the dwell time optimisation algorthim. The focus is on the optimisation method, and as such it is assumed that each seed radially deposits its energy purely based upon the inverse square law.

Two regions are defined, a central target region and a off centre avoid region.

Initialisation



In [1]:

    
import numpy as np

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.colors import BoundaryNorm
from matplotlib.colorbar import ColorbarBase
%matplotlib inline

from utilities import BasinhoppingWrapper, create_green_cm



In [2]:

    
green_cm = create_green_cm()

Create the calculation grid



In [3]:

    
grid_spacing = 0.1 # Was 0.1, using 0.2 to speed up for testing

x_ = np.arange(-0.5, 0.5 + grid_spacing, grid_spacing)
y_ = np.arange(-0.5, 0.5 + grid_spacing, grid_spacing)
z_ = np.arange(-0.5, 0.5 + grid_spacing, grid_spacing)

x, y, z = np.meshgrid(x_, y_, z_)
x = np.ravel(x)
y = np.ravel(y)
z = np.ravel(z)

Define the target and avoid cubes



In [4]:

    
def target(x, y, z):
    return (
        (x < 0.45) & (x > -0.45) & 
        (y < 0.45) & (y > -0.45) & 
        (z < 0.45) & (z > -0.45))

def avoid(x, y, z):
    return (
        (x < 0.25) & (x > 0.05) & 
        (y < 0.15) & (y > -0.15) & 
        (z < 2) & (z > -2))

target_cube = target(x, y, z)

avoid_cube = avoid(x, y, z)

target_cube = target_cube & ~avoid_cube

Display the target and avoid cubes



In [5]:

    
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.scatter(
    x[target_cube], y[target_cube], z[target_cube], 
    alpha=0.2, s=2, color='blue')
ax.scatter(
    x[avoid_cube], y[avoid_cube], z[avoid_cube], 
    alpha=0.2, s=2, color='red')

ax.set_xlim([np.min(x_), np.max(x_)])
ax.set_ylim([np.min(y_), np.max(y_)])
ax.set_zlim([np.min(z_), np.max(z_)])









    Out[5]:





(-0.5, 0.49999999999999978)



In [6]:

    
ax.view_init(azim=-90, elev=90)
fig









    Out[6]:

Create initial equidistant parrallel lines to represent catheters. Inbuild a slight random skew to emulate what occurs physically



In [7]:

    
number_of_lines = 16 # Was 16, reduced to 9 for testing

line_start = np.meshgrid(
    [-0.39, -0.13, 0.13, 0.39],
    [-0.39, -0.13, 0.13, 0.39],
    [1])

line_finish = np.array([
    line_start[0] + np.random.normal(scale=0.02, size=[4, 4, 1]),
    line_start[1] + np.random.normal(scale=0.02, size=[4, 4, 1]),
    -line_start[2]])



In [8]:

    
line_start = np.array([np.ravel(mesh) for mesh in line_start])
line_finish = np.array([np.ravel(mesh) for mesh in line_finish])

Display the lines overlayed



In [9]:

    
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.scatter(
    x[target_cube], y[target_cube], z[target_cube], 
    alpha=0.2, s=2, color='blue')
ax.scatter(
    x[avoid_cube], y[avoid_cube], z[avoid_cube], 
    alpha=0.2, s=2, color='red')
ax.scatter(*line_start)
ax.scatter(*line_finish)

for i in range(len(line_start[0])):
    plt_coords = [
        [line_start[j][i], line_finish[j][i]]
        for j in range(len(line_start))]
    ax.plot(*plt_coords, color='black', alpha=0.5)

ax.set_xlim([-1, 1])
ax.set_ylim([-1, 1])
ax.set_zlim([-1, 1])









    Out[9]:





(-1, 1)

Create a function to return x, y, z coords when a distance along a line is requested



In [10]:

    
diff = (line_finish - line_start)
line_length = np.sqrt(diff[0]**2 + diff[1]**2 + diff[2]**2)

def find_distance_coords(line_num=None, distance=None):
    relative_dist = distance / line_length[line_num]
    
    if (relative_dist > 1) | (relative_dist < 0):
        return np.array([np.nan]*3)
    
    x = (
        line_start[0][line_num] * (1 - relative_dist) + 
        line_finish[0][line_num] * relative_dist)
    
    y = (
        line_start[1][line_num] * (1 - relative_dist) + 
        line_finish[1][line_num] * relative_dist)
        
    z = (
        line_start[2][line_num] * (1 - relative_dist) + 
        line_finish[2][line_num] * relative_dist)
    
    coords = np.array([x, y, z])
    
    return coords

Pick dwell positons starting at a random position along the line



In [11]:

    
dwell_spacing = 0.1 # Was 0.1, increased to 0.3 for testing
dwell_distances_from_initial = np.arange(0, 2, dwell_spacing)
number_of_dwells = len(dwell_distances_from_initial)



In [12]:

    
inital_dwell_position = np.random.uniform(
    low=0, high=dwell_spacing, size=number_of_lines)
inital_dwell_position









    Out[12]:





array([ 0.00623742,  0.02435552,  0.09851191,  0.04731143,  0.04131299,
        0.09314668,  0.04570622,  0.02615803,  0.07655093,  0.0195772 ,
        0.07244818,  0.0767939 ,  0.03721719,  0.0875253 ,  0.04447357,
        0.07471193])



In [13]:

    
dwell_distances = np.reshape(inital_dwell_position, (-1, 1)) + np.reshape(dwell_distances_from_initial, (1, -1))



In [14]:

    
def find_dwell_coords(line_num=None, dwell_num=None):
    distance = dwell_distances[line_num, dwell_num]
    
    coords = find_distance_coords(
        line_num=line_num, distance=distance)
    
    return coords

Find all the dwell positions that are on the grid



In [15]:

    
dwell_positions = np.array([
    [
        find_dwell_coords(
            line_num=line_num, dwell_num=dwell_num)
        for dwell_num in range(number_of_dwells)]
 for line_num in range(number_of_lines)])

Plot the dwell positions



In [16]:

    
line_colours = np.random.uniform(size=(number_of_lines,3))



In [17]:

    
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.scatter(
    x[target_cube], y[target_cube], z[target_cube], 
    alpha=0.2, s=2, color='blue')
ax.scatter(
    x[avoid_cube], y[avoid_cube], z[avoid_cube], 
    alpha=0.2, s=2, color='red')

for line_num in range(number_of_lines):
    ax.scatter(*np.transpose(dwell_positions[line_num]), 
               c=line_colours[line_num], alpha=0.7)


ax.set_xlim([-1, 1])
ax.set_ylim([-1, 1])
ax.set_zlim([-1, 1])









    Out[17]:





(-1, 1)

Select the dwell positions that fall withing the target region. Only use these dwell positions.



In [18]:

    
dwell_positions_to_be_filtered = np.reshape(dwell_positions, (-1, 3))

distance_to_dwell_pos_initial = np.array([
    np.sqrt(
        (x[i] - dwell_positions_to_be_filtered[:,0])**2 + 
        (y[i] - dwell_positions_to_be_filtered[:,1])**2 + 
        (z[i] - dwell_positions_to_be_filtered[:,2])**2
    )
    for i in range(len(x))
])



In [19]:

    
closest_voxel_to_dwell = np.reshape(np.argmin(distance_to_dwell_pos_initial, axis=0), [1, -1])
target_cube_voxels = np.reshape(np.where(target_cube)[0], [-1, 1])
is_dwell_closest_to_target = np.any(closest_voxel_to_dwell == target_cube_voxels, axis=0)

relevant_dwell_positions = dwell_positions_to_be_filtered[is_dwell_closest_to_target]



In [20]:

    
line_number_index = np.reshape(np.arange(number_of_lines), (number_of_lines, 1))
line_number_index = np.ones([number_of_lines, number_of_dwells]) * line_number_index
line_number_index = np.ravel(line_number_index)[is_dwell_closest_to_target]
line_number_index = line_number_index.astype(int)



In [21]:

    
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.scatter(
    x[target_cube], y[target_cube], z[target_cube], 
    alpha=0.2, s=2, color='blue')

ax.scatter(
    x[avoid_cube], y[avoid_cube], z[avoid_cube], 
    alpha=0.2, s=2, color='red')

for line_num in range(number_of_lines):
    ref = line_number_index == line_num
    ax.scatter(*np.transpose(relevant_dwell_positions[ref]), 
               c=line_colours[line_num], alpha=0.7)

ax.set_xlim([np.min(x_), np.max(x_)])
ax.set_ylim([np.min(y_), np.max(y_)])
ax.set_zlim([np.min(z_), np.max(z_)])









    Out[21]:





(-0.5, 0.49999999999999978)

Calculate distance between each dwell and each grid point

Create an array containing the distance to each dwell position for each voxel and translate this to exposure at that voxel per unit dwell time at each dwell position. This is defined in this way so that a reasonable portion of the calculation of exposure can be done outside of the optimisation.



In [22]:

    
distance_to_dwell_pos = np.array([
    np.sqrt(
        (x[i] - relevant_dwell_positions[:,0])**2 + 
        (y[i] - relevant_dwell_positions[:,1])**2 + 
        (z[i] - relevant_dwell_positions[:,2])**2
    )
    for i in range(len(x))
])

exposure_per_unit_time = 1 / distance_to_dwell_pos**2

Estimate uncertainty in distance

Estimate the uncertainty in distance given a shift along the catheter



In [23]:

    
relevant_dwell_distances = np.ravel(dwell_distances)[is_dwell_closest_to_target]

shift_uncertainty = 0.05

shifted_up = np.array([
    find_distance_coords(
        distance=relevant_dwell_distances[i] + shift_uncertainty, 
        line_num=line_number_index[i])
    for i in range(len(relevant_dwell_distances))
])

shifted_down = np.array([
    find_distance_coords(
        distance=relevant_dwell_distances[i] - shift_uncertainty, 
        line_num=line_number_index[i])
    for i in range(len(relevant_dwell_distances))
])



In [24]:

    
distance_to_dwell_shift_up = np.array([
    np.sqrt(
        (x[i] - shifted_up[:,0])**2 + 
        (y[i] - shifted_up[:,1])**2 + 
        (z[i] - shifted_up[:,2])**2
    )
    for i in range(len(x))
])

distance_to_dwell_shift_down = np.array([
    np.sqrt(
        (x[i] - shifted_down[:,0])**2 + 
        (y[i] - shifted_down[:,1])**2 + 
        (z[i] - shifted_down[:,2])**2
    )
    for i in range(len(x))
])



In [25]:

    
dwell_distance_uncertainty = (
    np.abs(distance_to_dwell_shift_up - distance_to_dwell_pos) + 
    np.abs(distance_to_dwell_shift_down - distance_to_dwell_pos)
) / 2

Create the PDF for exposure given a distance and uncertainty



In [45]:

    
num_relevant_dwells = np.sum(is_dwell_closest_to_target)

def create_pdf_given_time(exposure, distance, uncertainty):   
    exposure = np.reshape(exposure, [1, 1, -1])
    distance = np.expand_dims(distance, axis=2)
    uncertainty = np.expand_dims(uncertainty, axis=2)    
    
    bot = 2*np.sqrt(2 * np.pi) * uncertainty * exposure**1.5
    
    def pdf(time):
        time = np.reshape(time, [1, -1, 1])
        top = (
            np.sqrt(time) * 
            np.exp(
                -(np.sqrt(time/exposure) - distance)**2 / 
                (2*uncertainty**2)))
    
        return top/bot
    
    return pdf



In [62]:

    
n = 2**9
exposure_spacing = np.linspace(0.01, 100, n)
pdf = create_pdf_given_time(exposure_spacing, distance_to_dwell_pos, dwell_distance_uncertainty)



In [70]:

    
time = np.ones(num_relevant_dwells)

check = pdf(time)









    Out[70]:





[<matplotlib.lines.Line2D at 0xc6062b0>]



In [83]:

    
for i in range(num_relevant_dwells):
    plt.plot(exposure_spacing, check[6,i,:])
    plt.ylim([0, 0.4])



In [64]:

    
np.shape(check[0,:,:])









    Out[64]:





(126, 512)



In [69]:

    
fft_of_rows = np.fft.rfft(check[2,:,:])
fft_of_convolution = fft_of_rows.prod(axis=0)
convolution_not_normalised = np.fft.irfft(fft_of_convolution)
# convolution = convolution_not_normalised[0:n] / np.trapz(convolution_not_normalised[0:n], x=exposure_spacing)

# plt.plot(exposure_spacing, convolution)

plt.plot(exposure_spacing, convolution_not_normalised)









    Out[69]:





[<matplotlib.lines.Line2D at 0xc5a84a8>]



In [ ]:



In [ ]:



In [ ]:



In [ ]:



In [ ]:



In [ ]:

Create exposure calculation function



In [23]:

    
def calculate_exposure(dwell_times):
    exposure = np.sum(dwell_times * exposure_per_unit_time, axis=1)
    
    return exposure

Run a test of arbitrary dwell times



In [24]:

    
num_relevant_dwells = len(relevant_dwell_positions)

random_pick = np.random.uniform(
    size=2, high=num_relevant_dwells, low=0).astype(int)

dwell_times = np.zeros([1, num_relevant_dwells])
dwell_times[0, random_pick] = 10

exposure = calculate_exposure(dwell_times)



In [25]:

    
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

reference = exposure > 80
colour = exposure[reference]
colour[colour > 200] = 200

ax.scatter(
    x[reference], y[reference], z[reference], 
    alpha=0.2, s=20, c=colour, cmap=green_cm)


ax.set_xlim([np.min(x_), np.max(x_)])
ax.set_ylim([np.min(y_), np.max(y_)])
ax.set_zlim([np.min(z_), np.max(z_)])









    Out[25]:





(-0.5, 0.49999999999999978)

Create function to display the results of the optimisation as it is being calculated.



In [26]:

    
def display_results(dwell_times):
    dwell_times = np.reshape(dwell_times, (1, num_relevant_dwells))
    exposure = calculate_exposure(dwell_times)
    
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    reference = exposure > 25
    colour = exposure[reference]
    colour[colour > 100] = 100
    
    small = exposure[reference] < 50
    large = ~small

    ax.scatter(
        x[reference][small], y[reference][small], z[reference][small], 
        alpha=0.2, s=3, c=colour[small], cmap=green_cm)
    
    ax.scatter(
        x[reference][large], y[reference][large], z[reference][large], 
        alpha=0.4, s=20, c=colour[large], cmap=green_cm)


    ax.set_xlim([np.min(x_), np.max(x_)])
    ax.set_ylim([np.min(y_), np.max(y_)])
    ax.set_zlim([np.min(z_), np.max(z_)])
    
    cost_function(dwell_times, debug=True)
    
    plt.show()

Create the optimisation cost function

The maximum exposure to a grid position not within the target is aimed to be less than 100. This cost function aims to achieve this.



In [27]:

    
def hot_exterior_cost_function(max_target_exterior):
    return ((max_target_exterior)/50)**4

testx = np.linspace(0, 120)
testy = hot_exterior_cost_function(testx)
plt.plot(testx, testy)

plt.ylim([0, 20])









    Out[27]:





(0, 20)

No grid point within the target is to be less than about 50. This cost function aims to achieve this.



In [28]:

    
def cold_target_cost_function(min_target):
    return ((min_target-80)/15)**4

testx = np.linspace(0, 120)
testy = cold_target_cost_function(testx)
plt.plot(testx, testy)

plt.ylim([0, 20])









    Out[28]:





(0, 20)

The avoid cost function aims to make no point within the avoid structure more than 45.



In [29]:

    
def hot_avoid_cost_function(max_avoid):
    return ((max_avoid)/25)**4

testx = np.linspace(0, 120)
testy = hot_avoid_cost_function(testx)
plt.plot(testx, testy)

plt.ylim([0, 20])









    Out[29]:





(0, 20)

Create the cost function to be used by the optimiser



In [30]:

    
def cost_function(dwell_times, debug=False):
    dwell_times = np.reshape(dwell_times, (1, num_relevant_dwells))
    exposure = calculate_exposure(dwell_times)
    
    min_target = np.min(exposure[target_cube])
    max_avoid = np.max(exposure[avoid_cube])
    max_target_exterior = np.max(exposure[~target_cube])
    
    cold_target_cost = cold_target_cost_function(min_target)
    hot_exterior_cost = hot_exterior_cost_function(max_target_exterior)
    hot_avoid_cost = hot_avoid_cost_function(max_avoid)
    
    total_cost = hot_exterior_cost + cold_target_cost + hot_avoid_cost
    
    if debug:
        print("Minimum target = %0.4f, resulting cost = %0.4f" %
              (min_target, cold_target_cost))
        print("Maximum exterior = %0.4f, resulting cost = %0.4f" %
              (max_target_exterior, hot_exterior_cost))
        print("Maximum avoid = %0.4f, resulting cost = %0.4f" %
              (max_avoid, hot_avoid_cost))
        print("Total cost = %0.4f" % (total_cost))
        
    
    return total_cost

Create initial conditions



In [31]:

    
num_relevant_dwells









    Out[31]:





126



In [32]:

    
initial_conditions = np.ones(num_relevant_dwells)*0.1

Step noise



In [33]:

    
step_noise = np.ones(num_relevant_dwells) * 0.3

Bounds



In [34]:

    
bounds = ((0, None),)*num_relevant_dwells

Run the optimiser



In [35]:

    
optimisation = BasinhoppingWrapper(
    to_minimise=cost_function,
    initial=initial_conditions,
    step_noise=step_noise,
    basinhopping_confidence=5,
    optimiser_confidence=0.0001,
    n=5,
    debug=display_results,
    bounds=bounds
)









    



---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-35-c46968f8b13d> in <module>()
      7     n=5,
      8     debug=display_results,
----> 9     bounds=bounds
     10 )

C:\Users\sbiggs\Documents\GitHub\poc-brachyoptimisation\utilities.py in __init__(self, n, optimiser_confidence, basinhopping_confidence, debug, bounds, **kwargs)
     42             )
     43 
---> 44         self.result = self.run_basinhopping()
     45 
     46     def step_function(self, optimiser_input):

C:\Users\sbiggs\Documents\GitHub\poc-brachyoptimisation\utilities.py in run_basinhopping(self)
     96             minimizer_kwargs=minimizer_config,
     97             take_step=self.step_function,
---> 98             callback=self.callback_function
     99         )
    100 

C:\Users\sbiggs\AppData\Local\Continuum\Anaconda3\lib\site-packages\scipy\optimize\_basinhopping.py in basinhopping(func, x0, niter, T, stepsize, minimizer_kwargs, take_step, accept_test, callback, interval, disp, niter_success)
    605 
    606     bh = BasinHoppingRunner(x0, wrapped_minimizer, take_step_wrapped,
--> 607                             accept_tests, disp=disp)
    608 
    609     # start main iteration loop

C:\Users\sbiggs\AppData\Local\Continuum\Anaconda3\lib\site-packages\scipy\optimize\_basinhopping.py in __init__(self, x0, minimizer, step_taking, accept_tests, disp)
     70 
     71         # do initial minimization
---> 72         minres = minimizer(self.x)
     73         if not minres.success:
     74             self.res.minimization_failures += 1

C:\Users\sbiggs\AppData\Local\Continuum\Anaconda3\lib\site-packages\scipy\optimize\_basinhopping.py in __call__(self, x0)
    277             return self.minimizer(x0, **self.kwargs)
    278         else:
--> 279             return self.minimizer(self.func, x0, **self.kwargs)
    280 
    281 

C:\Users\sbiggs\AppData\Local\Continuum\Anaconda3\lib\site-packages\scipy\optimize\_minimize.py in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
    442     elif meth == 'l-bfgs-b':
    443         return _minimize_lbfgsb(fun, x0, args, jac, bounds,
--> 444                                 callback=callback, **options)
    445     elif meth == 'tnc':
    446         return _minimize_tnc(fun, x0, args, jac, bounds, callback=callback,

C:\Users\sbiggs\AppData\Local\Continuum\Anaconda3\lib\site-packages\scipy\optimize\lbfgsb.py in _minimize_lbfgsb(fun, x0, args, jac, bounds, disp, maxcor, ftol, gtol, eps, maxfun, maxiter, iprint, callback, **unknown_options)
    318                 # minimization routine wants f and g at the current x
    319                 # Overwrite f and g:
--> 320                 f, g = func_and_grad(x)
    321         elif task_str.startswith(b'NEW_X'):
    322             # new iteration

C:\Users\sbiggs\AppData\Local\Continuum\Anaconda3\lib\site-packages\scipy\optimize\lbfgsb.py in func_and_grad(x)
    265         def func_and_grad(x):
    266             f = fun(x, *args)
--> 267             g = _approx_fprime_helper(x, fun, epsilon, args=args, f0=f)
    268             return f, g
    269     else:

C:\Users\sbiggs\AppData\Local\Continuum\Anaconda3\lib\site-packages\scipy\optimize\optimize.py in _approx_fprime_helper(xk, f, epsilon, args, f0)
    556         ei[k] = 1.0
    557         d = epsilon * ei
--> 558         grad[k] = (f(*((xk + d,) + args)) - f0) / d[k]
    559         ei[k] = 0.0
    560     return grad

C:\Users\sbiggs\AppData\Local\Continuum\Anaconda3\lib\site-packages\scipy\optimize\optimize.py in function_wrapper(*wrapper_args)
    283     def function_wrapper(*wrapper_args):
    284         ncalls[0] += 1
--> 285         return function(*(wrapper_args + args))
    286 
    287     return ncalls, function_wrapper

<ipython-input-30-140b2338eb5f> in cost_function(dwell_times, debug)
      1 def cost_function(dwell_times, debug=False):
      2     dwell_times = np.reshape(dwell_times, (1, num_relevant_dwells))
----> 3     exposure = calculate_exposure(dwell_times)
      4 
      5     min_target = np.min(exposure[target_cube])

<ipython-input-23-d98bbfa96ff8> in calculate_exposure(dwell_times)
      1 def calculate_exposure(dwell_times):
----> 2     exposure = np.sum(dwell_times * exposure_per_unit_time, axis=1)
      3 
      4     return exposure

KeyboardInterrupt:

Presentation of results



In [ ]:

    
display_results(optimisation.result)

Display histogram of resulting dwell times



In [ ]:

    
plt.hist(optimisation.result)

Give overview of dwell times segmented by catheter



In [ ]:

    
for line_num in range(number_of_lines):
    
    print("Line number %d" % (line_num))
    
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    ax.scatter(
        x[target_cube], y[target_cube], z[target_cube], 
        alpha=0.2, s=2, color='blue')
    ax.scatter(
        x[avoid_cube], y[avoid_cube], z[avoid_cube], 
        alpha=0.2, s=2, color='red')
    
    ref = line_number_index == line_num
    ax.scatter(*np.transpose(relevant_dwell_positions[ref]), 
               c=line_colours[line_num], alpha=0.7)


    ax.set_xlim([np.min(x_), np.max(x_)])
    ax.set_ylim([np.min(y_), np.max(y_)])
    ax.set_zlim([np.min(z_), np.max(z_)])
    
    plt.show() 
    
    
    plt.scatter(np.arange(np.sum(ref)), optimisation.result[ref],
               c=line_colours[line_num])
    plt.xlabel("Dwell number")
    plt.ylabel("Units of dwell time")
    
    plt.ylim(bottom=0)
    plt.show()

Convert result into column vector for post analysis



In [ ]:

    
dwell_times = np.reshape(optimisation.result, (1, num_relevant_dwells))

Create a custom resolution exposure calculation function for post analysis



In [ ]:

    
def create_calculate_exposure(dwell_positions, 
                              x_new=None, y_new=None, z_new=None):
    x, y, z = np.meshgrid(
        x_new, y_new, z_new)
    x = np.ravel(x)
    y = np.ravel(y)
    z = np.ravel(z)

    distance_to_dwell_pos = np.array([
        np.sqrt(
            (x[i] - dwell_positions[:,0])**2 + 
            (y[i] - dwell_positions[:,1])**2 + 
            (z[i] - dwell_positions[:,2])**2
        )
        for i in range(len(x))
    ])

    exposure_per_unit_time = 1 / distance_to_dwell_pos**2
    
    def calculate_exposure(dwell_times, reshape=False):
        exposure = np.sum(dwell_times * exposure_per_unit_time, axis=1)
        if reshape:
            exposure = np.reshape(exposure, (
                len(x_new), len(y_new), len(z_new)))
        return exposure
        
    return calculate_exposure



In [ ]:

    
dx = 0.01; dy = 0.01; dz = 0.1

contour_x_ = np.arange(-1, 1 + dx, dx)
contour_y_ = np.arange(-1, 1 + dy, dy)
contour_z_ = z_ #np.arange(-0.7, 0.7 + dz, dz)

contour_calculate_exposure = create_calculate_exposure(
    relevant_dwell_positions, x_new=contour_x_, y_new=contour_y_, z_new=contour_z_)

contour_exposure_init = contour_calculate_exposure(dwell_times, reshape=True)



In [ ]:

    
vmin = 0
vmax = 200
num_contours = 41

fig = plt.figure(figsize=(10,0.3))
cbar_ax = fig.add_subplot(1,1,1)

bounds = np.linspace(vmin, vmax, num_contours)
norm = BoundaryNorm(bounds, green_cm.N)


cb = ColorbarBase(
    cbar_ax, norm=norm, cmap=green_cm,
    orientation='horizontal')

cb.ax.xaxis.set_ticks_position('top')
cb.ax.xaxis.set_label_position('top')

cb.set_clim([vmin, vmax])


contour_exposure = contour_exposure_init.copy()
contour_exposure[contour_exposure > vmax] = vmax


for i, z_value in enumerate(contour_z_):
    plt.figure(figsize=(6,6))
    
    c = plt.contourf(contour_x_, contour_y_,
        contour_exposure[:, :, i], num_contours,
        vmin=vmin, vmax=vmax, cmap=green_cm)
    
    reference = z[target_cube] == z_value
    x_target = x[target_cube][reference]
    y_target = y[target_cube][reference]
    plt.scatter(x_target, y_target, alpha=0.3, color='blue')
    
    reference = z[avoid_cube] == z_value
    x_avoid = x[avoid_cube][reference]
    y_avoid = y[avoid_cube][reference]
    plt.scatter(x_avoid, y_avoid, alpha=0.3, color='red')
    
    plt.title("Slice z = %0.1f" % (z_value))
    plt.xlim([np.min(contour_x_), np.max(contour_x_)])
    plt.ylim([np.min(contour_y_), np.max(contour_y_)])
    plt.show()

Analysis of results

Create analysis grid



In [ ]:

    
dx = 0.02
dy = 0.02
dz = 0.02

post_x_ = np.arange(-1, 1 + dx, dx)
post_y_ = np.arange(-1, 1 + dy, dy)
post_z_ = np.arange(-1, 1 + dz, dz)

post_x, post_y, post_z = np.meshgrid(
    post_x_, post_y_, post_z_)
post_x = np.ravel(post_x)
post_y = np.ravel(post_y)
post_z = np.ravel(post_z)

post_calculate_exposure = create_calculate_exposure(
    relevant_dwell_positions, x_new=post_x_, y_new=post_y_, z_new=post_z_)

post_exposure = post_calculate_exposure(dwell_times)

Define structures on new finer grid



In [ ]:

    
post_target_cube = target(post_x, post_y, post_z)
post_avoid_cube = avoid(post_x, post_y, post_z)

Make a function to plot DVHs given a reference



In [ ]:

    
def create_dvh(reference, exposure, ax=None):
    if ax is None:
        fig = plt.figure()
        ax = fig.add_subplot(1,1,1)
    
    results = exposure[reference]
    results[results>150] = 150
    hist = np.histogram(results, 100)
    
    freq = hist[0]
    bin_edge = hist[1]
    bin_mid = (bin_edge[1::] + bin_edge[:-1:])/2
    
    cumulative = np.cumsum(freq[::-1])
    cumulative = cumulative[::-1]
    bin_mid = np.append([0], bin_mid)
    
    cumulative = np.append(cumulative[0], cumulative)
    percent_cumulative = cumulative / cumulative[0] * 100

    ax.plot(bin_mid, percent_cumulative)

Plotting of relevant DVHs



In [ ]:

    
fig = plt.figure()
ax = fig.add_subplot(1,1,1)

create_dvh(post_target_cube, post_exposure, ax=ax)
create_dvh(post_avoid_cube, post_exposure, ax=ax)
create_dvh(~post_target_cube, post_exposure, ax=ax)

plt.legend(["Target", "Avoid", "Patient"])

plt.xlim([0, 150])

Add a random shift to each catheter, see the effect it has on target DVH.



In [ ]:

    
shift_uncertainty = 0.05



In [ ]:

    
line_displacement = np.random.normal(scale=shift_uncertainty, size=number_of_lines)



In [ ]:

    
shifted_dwell_positions = np.array([
    find_distance_coords(
        distance=relevant_dwell_distances[i] + line_displacement[line_number_index[i]], 
        line_num=line_number_index[i])
    for i in range(len(relevant_dwell_distances))
])



In [ ]:

    
shifted_calculate_exposure = create_calculate_exposure(
    shifted_dwell_positions, x_new=post_x_, y_new=post_y_, z_new=post_z_)

shifted_exposure = shifted_calculate_exposure(dwell_times)

With the shift the change in the DVH is the following. For the optimisation given here there is very little change in DVH.



In [ ]:

    
fig = plt.figure()
ax = fig.add_subplot(1,1,1)

create_dvh(post_target_cube, post_exposure, ax=ax)
create_dvh(post_target_cube, shifted_exposure, ax=ax)

plt.legend(["Original", "With position error"], loc='lower left')

plt.xlim([0, 110])
# plt.ylim([95, 100])



In [ ]:

    
np.min(shifted_exposure[post_target_cube])



In [ ]:

    
np.min(post_exposure[post_target_cube])



In [ ]:

    
fig = plt.figure()
ax = fig.add_subplot(1,1,1)

create_dvh(post_avoid_cube, post_exposure, ax=ax)
create_dvh(post_avoid_cube, shifted_exposure, ax=ax)

plt.legend(["Original", "With position error"], loc='upper right')

plt.xlim([0, 110])



In [ ]:



In [ ]:

Overview of shift results



In [ ]:

    
exposure = calculate_exposure(dwell_times)
np.min(exposure[target_cube])



In [ ]:

    
def return_cold_spot_in_random_displace():
    line_displacement = np.random.normal(scale=shift_uncertainty, size=number_of_lines)

    shifted_dwell_positions = np.array([
        find_distance_coords(
            distance=relevant_dwell_distances[i] + line_displacement[line_number_index[i]], 
            line_num=line_number_index[i])
        for i in range(len(relevant_dwell_distances))
    ])

    shifted_calculate_exposure = create_calculate_exposure(
        shifted_dwell_positions, x_new=x_, y_new=y_, z_new=z_)

    shifted_exposure = shifted_calculate_exposure(dwell_times)

    return np.min(shifted_exposure[target_cube])



In [ ]:

    
n = 1000
cold_spots = np.zeros(n)

for i in range(n):
    cold_spots[i] = return_cold_spot_in_random_displace()



In [ ]:

    
plt.hist(cold_spots,20);



In [ ]:



In [ ]:



In [ ]:



In [ ]:



In [ ]:



In [ ]: