Skip to content

Commit 03fcb6c

Browse files
gourav3017gourav3017
committed
1- Update SCP based VMAT optimization.
2- Get voxel coordinates 2- Use toml instead of setup.py for build 3- Fix bug in structures.py for "-" in struct name. 4- Fix visualization.py.
1 parent 59a1d7f commit 03fcb6c

File tree

16 files changed

+1567
-301
lines changed

16 files changed

+1567
-301
lines changed

‎README.md‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -14,8 +14,7 @@
1414
[![Total Downloads](https://static.pepy.tech/personalized-badge/portpy?period=total&units=international_system&left_color=grey&right_color=blue&left_text=total%20downloads)](https://pepy.tech/project/portpy?&left_text=totalusers)
1515
[![Monthly Downloads](https://static.pepy.tech/badge/portpy/month)](https://pepy.tech/project/portpy)
1616

17-
**WOC Participants:** To explore project details, please visit the WOC subfolder using this link: ([WOC Projects](https://github.com/PortPy-Project/PortPy/tree/master/WOC))
18-
17+
**Note**: If you have any questions about PortPy, please create an [Issue](https://github.com/PortPy-Project/PortPy/issues) on Github. If you are unable to do so for any reason, you may contact Gourav Jhanwar ([email protected])
1918

2019
# What is PortPy? <a name="What"></a>
2120

‎examples/python_files/vmat_scp_tutorial.py‎

Lines changed: 28 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,7 @@ def vmat_scp_tutorial():
2727
data = pp.DataExplorer(data_dir=data_dir)
2828

2929
# pick a patient from the existing patient list to get detailed info (e.g., beam angles, structures).
30-
data.patient_id = 'Lung_Patient_3'
30+
data.patient_id = 'Lung_Phantom_Patient_1'
3131

3232
# display the data of the patient in console or browser.
3333
data.display_patient_metadata()
@@ -57,6 +57,10 @@ def vmat_scp_tutorial():
5757
# Loading influence matrix
5858
inf_matrix = pp.InfluenceMatrix(ct=ct, structs=structs, beams=beams)
5959

60+
# scale influence matrix to get correct MU/deg for TPS import
61+
inf_matrix_scale_factor = vmat_opt_params['opt_parameters'].get('inf_matrix_scale_factor', 1)
62+
inf_matrix.A = inf_matrix.A * np.float32(inf_matrix_scale_factor)
63+
6064
"""
6165
2) Creating Arc (Arcs class)
6266
@@ -75,12 +79,33 @@ def vmat_scp_tutorial():
7579
3) Optimize VMAT plan using sequential convex programming
7680
7781
"""
82+
# If column generation flag is turned in opt_params. Generate initial feasible leaf positions for SCP based optimization
83+
if vmat_opt_params['opt_parameters']['initial_leaf_pos'].lower() == 'cg':
84+
vmat_opt = pp.VmatOptimizationColGen(my_plan=my_plan,
85+
opt_params=vmat_opt_params)
86+
87+
sol_col_gen = vmat_opt.run_col_gen_algo(solver='MOSEK', verbose=True, accept_unknown=True)
88+
dose_1d = inf_matrix.A @ sol_col_gen['optimal_intensity'] * my_plan.get_num_of_fractions()
89+
# # plot dvh for the above structures
90+
fig, ax = plt.subplots(figsize=(12, 8))
91+
struct_names = ['PTV', 'ESOPHAGUS', 'HEART', 'CORD', 'RIND_0', 'RIND_1', 'LUNGS_NOT_GTV', 'RECT_WALL', 'BLAD_WALL',
92+
'URETHRA']
93+
ax = pp.Visualization.plot_dvh(my_plan, dose_1d=dose_1d,
94+
struct_names=struct_names, ax=ax)
95+
ax.set_title('Initial Col gen dvh')
96+
plt.show(block=False)
97+
plt.close('all')
98+
# pp.save_optimal_sol(sol=sol_col_gen, sol_name='sol_col_gen', path=r'C:\Temp')
99+
100+
78101
# Initialize Optimization
79102
vmat_opt = pp.VmatScpOptimization(my_plan=my_plan,
80103
opt_params=vmat_opt_params)
81104
# Run Sequential convex algorithm for optimising the plan.
82105
# The final result will be stored in sol and convergence will store the convergence history (i.e., results of each iteration)
83-
sol, convergence = vmat_opt.run_sequential_cvx_algo(solver='MOSEK', verbose=True)
106+
convergence = vmat_opt.run_sequential_cvx_algo(solver='MOSEK', verbose=True)
107+
sol = convergence[vmat_opt.best_iteration]
108+
sol['inf_matrix'] = inf_matrix
84109

85110
# Visualize convergence. The convergence dataframe contains the following columns:
86111
df = pd.DataFrame(convergence, columns=['outer_iteration', 'inner_iteration', 'step_size_f_b', 'forward_backward', 'intermediate_obj_value', 'actual_obj_value', 'accept'])
@@ -96,7 +121,7 @@ def vmat_scp_tutorial():
96121
# plot dvh for the above structures
97122
fig, ax = plt.subplots(figsize=(12, 8))
98123
struct_names = ['PTV', 'ESOPHAGUS', 'HEART', 'CORD', 'LUNGS_NOT_GTV']
99-
pp.Visualization.plot_dvh(my_plan, sol=sol, struct_names=struct_names, style='dashed', ax=ax)
124+
pp.Visualization.plot_dvh(my_plan, sol=sol, struct_names=struct_names, style='solid', ax=ax)
100125
plt.show()
101126
print('Done')
102127

‎examples/vmat_scp_tutorial.ipynb‎

Lines changed: 31 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -85,7 +85,8 @@
8585
"source": [
8686
"import portpy.photon as pp\n",
8787
"import numpy as np\n",
88-
"import os"
88+
"import os\n",
89+
"import matplotlib.pyplot as plt"
8990
]
9091
},
9192
{
@@ -296,7 +297,12 @@
296297
" clinical_criteria=clinical_criteria)\n",
297298
"\n",
298299
"# Loading influence matrix\n",
299-
"inf_matrix = pp.InfluenceMatrix(ct=ct, structs=structs, beams=beams)"
300+
"inf_matrix = pp.InfluenceMatrix(ct=ct, structs=structs, beams=beams)\n",
301+
"\n",
302+
"# scale influence matrix to get correct MU/deg for TPS import\n",
303+
"inf_matrix_scale_factor = vmat_opt_params['opt_parameters'].get('inf_matrix_scale_factor', 1)\n",
304+
"print(\"inf_matrix_scale_factor: \", inf_matrix_scale_factor)\n",
305+
"inf_matrix.A = inf_matrix.A * np.float32(inf_matrix_scale_factor)"
300306
]
301307
},
302308
{
@@ -334,6 +340,12 @@
334340
"\n"
335341
]
336342
},
343+
{
344+
"metadata": {},
345+
"cell_type": "markdown",
346+
"source": "If initial_leaf_pos is set to 'CG'(Column Generation). Generate initial leaf positions using column generation for SCP based optimization. This plan is not necessarily optimal, but can provide good initial start for SCP based VMAT optimization. Default is set to 'CG'. Users can skip this by setting vmat_opt_params['opt_parameters']['initial_leaf_pos'] = 'random' or 'BEV'",
347+
"id": "80e58576f2b10775"
348+
},
337349
{
338350
"cell_type": "code",
339351
"execution_count": 18,
@@ -343,12 +355,28 @@
343355
},
344356
"outputs": [],
345357
"source": [
358+
"if vmat_opt_params['opt_parameters']['initial_leaf_pos'].lower() == 'cg':\n",
359+
" vmat_opt = pp.VmatOptimizationColGen(my_plan=my_plan,\n",
360+
" opt_params=vmat_opt_params)\n",
361+
"\n",
362+
" sol_col_gen = vmat_opt.run_col_gen_algo(solver='MOSEK', verbose=True, accept_unknown=True)\n",
363+
" dose_1d = inf_matrix.A @ sol_col_gen['optimal_intensity'] * my_plan.get_num_of_fractions()\n",
364+
" # # plot dvh for the above structures\n",
365+
" fig, ax = plt.subplots(figsize=(12, 8))\n",
366+
" struct_names = ['PTV', 'ESOPHAGUS', 'HEART', 'CORD', 'RIND_0', 'RIND_1', 'LUNGS_NOT_GTV']\n",
367+
" ax = pp.Visualization.plot_dvh(my_plan, dose_1d=dose_1d,\n",
368+
" struct_names=struct_names, ax=ax)\n",
369+
" ax.set_title('Initial Col gen dvh')\n",
370+
" plt.show(block=False)\n",
371+
" plt.close('all')\n",
346372
"\n",
347373
"# Initialize Optimization\n",
348374
"vmat_opt = pp.VmatScpOptimization(my_plan=my_plan,\n",
349375
" opt_params=vmat_opt_params)\n",
350376
"# Run Sequential convex algorithm for optimising the plan. The final result will be stored in sol and convergence will store the convergence history (i.e., results of each iteration)\n",
351-
"sol, convergence = vmat_opt.run_sequential_cvx_algo(solver='MOSEK', verbose=True)"
377+
"convergence = vmat_opt.run_sequential_cvx_algo(solver='MOSEK', verbose=True)\n",
378+
"sol = convergence[vmat_opt.best_iteration] # get best solution\n",
379+
"sol['inf_matrix'] = inf_matrix # reference to influence matrix object for further evaluation and visualization"
352380
]
353381
},
354382
{

‎portpy/config_files/optimization_params/optimization_params_Lung_2Gy_30Fx_vmat.json‎

Lines changed: 15 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -7,12 +7,6 @@
77
"weight": 10000,
88
"dose_perc": 100
99
},
10-
{
11-
"type": "linear-overdose",
12-
"structure_name": "CORD",
13-
"weight": 100,
14-
"dose_gy": 50
15-
},
1610
{
1711
"type": "quadratic-underdose",
1812
"structure_name": "PTV",
@@ -81,9 +75,13 @@
8175
"weight": 3
8276
},
8377
{
84-
"type": "smoothness-quadratic",
85-
"weight": 100
86-
}
78+
"type": "aperture_regularity_quadratic",
79+
"weight": 1000
80+
},
81+
{
82+
"type": "aperture_similarity_quadratic",
83+
"weight": 1000
84+
}
8785
],
8886
"constraints":[
8987
{
@@ -139,7 +137,7 @@
139137
],
140138
"opt_parameters":{
141139
"minimum_dynamic_leaf_gap_mm": 0.5,
142-
"initial_leaf_pos": "BEV",
140+
"initial_leaf_pos": "CG",
143141
"initial_step_size": 0,
144142
"step_size_increment": 1,
145143
"ss_termination": 1,
@@ -149,8 +147,8 @@
149147
"termination_gap": 2,
150148
"min_iteration_threshold": 10,
151149
"dose_threshold": 0.0001,
152-
"mu_max": 0.167,
153-
"mu_min": 0.001,
150+
"mu_max": 55,
151+
"mu_min": 0.5,
154152
"epsilon_mu": 1e-06,
155153
"forward_backward": 1,
156154
"total_iter_step1": 0,
@@ -159,8 +157,11 @@
159157
"first_beam_adj": 0,
160158
"second_beam_adj": 0.5,
161159
"last_beam_adj": 1,
162-
"w_ai_obj": 10,
163-
"w_as_obj": 20
160+
"smooth_delta": 0.008,
161+
"update_balanced_arc_score": 1,
162+
"inf_matrix_scale_factor": 0.6,
163+
"flag_full_matrix": 1,
164+
"sparsification": "rmr"
164165
}
165166

166167
}

‎portpy/photon/beam.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,7 @@ def __init__(self, data: DataExplorer, beam_ids: List[Union[int, str]] = None,
3838
beams_dict = data.load_data(meta_data=metadata['beams'], load_inf_matrix_full=load_inf_matrix_full)
3939
self.beams_dict = beams_dict
4040
self.preprocess_beams()
41+
self.patient_id = data.patient_id
4142

4243
def get_beamlet_idx_2d_finest_grid(self, beam_id: Union[int, str]) -> np.ndarray:
4344
"""

‎portpy/photon/structures.py‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -44,6 +44,7 @@ def __init__(self, data: DataExplorer) -> None:
4444
self.opt_voxels_dict['name'] = structures_dict['name']
4545
self._ct_voxel_resolution_xyz_mm = deepcopy(self.opt_voxels_dict['ct_voxel_resolution_xyz_mm'])
4646
self.preprocess_structures()
47+
self.patient_id = data.patient_id
4748

4849
def get_structures(self) -> list:
4950
"""
@@ -445,8 +446,9 @@ def evaluate_expression(self, expression):
445446
continue
446447
elif expression[i].isalpha() or expression[i].isdigit():
447448
j = i
448-
while j < len(expression) and (expression[j].isalpha() or expression[j] in '_.$^' or expression[j].isdigit()):
449-
j += 1
449+
while j < len(expression) and (expression[j].isalpha() or expression[j] in '-_.$^' or expression[j].isdigit()):
450+
if not ((j > 0 and expression[j - 1].isspace()) and (j + 1 < len(expression) and expression[j + 1].isspace())):
451+
j += 1
450452
struct = expression[i:j]
451453
values_stack.append(struct)
452454
i = j - 1

‎portpy/photon/utils/get_apt_reg_metric.py‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
# get smoothness function value using leaf positions in arcs
2+
def get_apt_reg_metric(my_plan, beam_ids=None):
3+
if beam_ids is None:
4+
beam_ids = [beam['beam_id'] for arc in my_plan.arcs.arcs_dict['arcs'] for beam in arc['vmat_opt']]
5+
arcs = my_plan.arcs
6+
apt_reg_obj = 0
7+
for arc in arcs.arcs_dict['arcs']:
8+
for beam in arc['vmat_opt']:
9+
if beam['beam_id'] in beam_ids:
10+
num_rows = beam['num_rows']
11+
for r, row in enumerate(beam['reduced_2d_grid']):
12+
if r < num_rows - 1:
13+
curr_left = beam['leaf_pos_left'][r] + 1
14+
curr_right = beam['leaf_pos_right'][r]
15+
next_left = beam['leaf_pos_left'][r + 1] + 1
16+
next_right = beam['leaf_pos_right'][r + 1]
17+
left_diff = abs(curr_left - next_left)
18+
right_diff = abs(curr_right - next_right)
19+
apt_reg_obj += left_diff ** 2 + right_diff ** 2
20+
# print('####apr_reg_obj: {} ####'.format(apt_reg_obj))
21+
return apt_reg_obj

‎portpy/photon/utils/get_sparse_only.py‎

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
1+
import numpy as np
2+
import math
3+
import scipy
4+
5+
6+
def get_sparse_only(A: np.ndarray, threshold_perc: float = 1, compression: str = 'naive', threshold_abs = None):
7+
"""
8+
Get sparse matrix using threshold and different methods
9+
:param A: matrix to be sparsified
10+
:param threshold_perc: threshold for matrix sparsification
11+
:param compression: Method of Sparsification
12+
:param threshold_abs: absolute threshold for matrix sparsification
13+
14+
:return: Sparse influence matrix
15+
"""
16+
threshold = np.max(A) * threshold_perc * 0.01
17+
if threshold_abs is not None:
18+
threshold = threshold_abs
19+
if compression == 'rmr':
20+
copy_matrix = A.copy()
21+
print('Generating sparse matrix using RMR...')
22+
np.apply_along_axis(row_operation, 1, copy_matrix, threshold)
23+
S = scipy.sparse.csr_matrix(copy_matrix)
24+
else:
25+
S = np.where(A >= threshold, A, 0)
26+
S = scipy.sparse.csr_matrix(S)
27+
return S
28+
29+
30+
def row_operation(copy_row, threshold):
31+
argzero = np.argwhere((np.abs(copy_row) <= threshold) * (copy_row != 0))
32+
argzero = argzero.reshape(len(argzero), )
33+
argzero_copy = copy_row[argzero]
34+
copy_row[argzero] = 0
35+
sum = np.sum(argzero_copy)
36+
if sum != 0:
37+
k = math.ceil(sum / threshold)
38+
39+
indices = np.random.choice(argzero, k, p=argzero_copy / sum, replace=True)
40+
np.add.at(copy_row, indices, sum / k)
41+

‎portpy/photon/utils/write_rt_plan_vmat.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -203,7 +203,7 @@ def write_rt_plan_vmat(my_plan: Plan, out_rt_plan_file: str, in_rt_plan_file: st
203203
for a, arc in enumerate(arcs.arcs_dict['arcs']):
204204
arc_id = arc['arc_id']
205205
beam_dataset = create_beam(my_plan=my_plan, arc_id=arc_id, beam_number=a+1)
206-
mu = [arc['vmat_opt'][i]['best_beam_weight'] for i in range(len(arc['vmat_opt']))]
206+
mu = [arc['vmat_opt'][i]['best_beam_weight']*100 for i in range(len(arc['vmat_opt']))] # multiply it by 100 to match eclipse mu/deg
207207
mu[0] = 0 # make first beam 0 mu
208208
total_mu = sum(mu)
209209
meterset_weight = 0

‎portpy/photon/visualization.py‎

Lines changed: 9 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -123,15 +123,15 @@ def plot_dvh(my_plan: Plan, sol: dict = None, dose_1d: np.ndarray = None, struct
123123
elif dose_scale == 'Relative(%)':
124124
x = x / pres * 100
125125
max_dose = np.maximum(max_dose, x[-1])
126-
ax.set_xlabel(r'Dose ($\%$)', fontsize=fontsize)
126+
ax.set_xlabel(r'Dose (%)', fontsize=fontsize)
127127

128128
if volume_scale == 'Absolute(cc)':
129129
y = y * my_plan.structures.get_volume_cc(all_orgs[i]) / 100
130130
max_vol = np.maximum(max_vol, y[1] * 100)
131131
ax.set_ylabel('Volume (cc)', fontsize=fontsize)
132132
elif volume_scale == 'Relative(%)':
133133
max_vol = np.maximum(max_vol, y[0] * 100)
134-
ax.set_ylabel(r'Fractional Volume ($\%$)', fontsize=fontsize)
134+
ax.set_ylabel(r'Fractional Volume (%)', fontsize=fontsize)
135135
ax.plot(x, 100 * y, linestyle=style, linewidth=width, color=colors[count], label=struct_names[i])
136136
count = count + 1
137137
# legend.append(struct_names[i])
@@ -290,15 +290,15 @@ def plot_robust_dvh(my_plan: Plan, sol: dict = None, dose_1d_list: list = None,
290290
elif dose_scale == 'Relative(%)':
291291
max_dose = np.maximum(max_dose, d_max_mat[-1])
292292
max_dose = max_dose / pres * 100
293-
ax.set_xlabel(r'Dose ($\%$)', fontsize=fontsize)
293+
ax.set_xlabel(r'Dose (%)', fontsize=fontsize)
294294

295295
if volume_scale == 'Absolute(cc)':
296296
y = y * my_plan.structures.get_volume_cc(all_orgs[i]) / 100
297297
max_vol = np.maximum(max_vol, y[1] * 100)
298298
ax.set_ylabel('Volume (cc)', fontsize=fontsize)
299299
elif volume_scale == 'Relative(%)':
300300
max_vol = np.maximum(max_vol, y[0] * 100)
301-
ax.set_ylabel(r'Fractional Volume ($\%$)', fontsize=fontsize)
301+
ax.set_ylabel(r'Fractional Volume (%)', fontsize=fontsize)
302302
# ax.plot(x, 100 * y, linestyle=style, linewidth=width, color=colors[count])
303303

304304
# ax.plot(d_min_mat, 100 * y, linestyle='dotted', linewidth=width*0.5, color=colors[count])
@@ -501,7 +501,7 @@ def plot_2d_slice(my_plan: Plan = None, sol: dict = None, dose_1d: np.ndarray =
501501
# adjust the main plot to make room for the legends
502502
plt.subplots_adjust(left=0.2)
503503
dose_3d = []
504-
if sol is not None:
504+
if sol is not None or dose_1d is not None:
505505
show_dose = True
506506
if show_dose:
507507
if sol is not None:
@@ -656,12 +656,13 @@ def get_colors():
656656
return colors
657657

658658
@staticmethod
659-
def surface_plot(matrix: np.ndarray, **kwargs):
659+
def surface_plot(matrix, ax=None, figsize=(8, 8), **kwargs):
660660
# acquire the cartesian coordinate matrices from the matrix
661661
# x is cols, y is rows
662662
(x, y) = np.meshgrid(np.arange(matrix.shape[0]), np.arange(matrix.shape[1]))
663-
fig = plt.figure()
664-
ax = fig.add_subplot(111, projection='3d')
663+
664+
if ax is None:
665+
fig, ax = plt.subplots(figsize=figsize, subplot_kw=dict(projection='3d'))
665666
surf = ax.plot_surface(x, y, np.transpose(matrix), **kwargs)
666667
return ax, surf
667668

@@ -706,14 +707,3 @@ def is_notebook() -> bool:
706707
return False # Probably standard Python interpreter
707708
except:
708709
return False # Probably standard Python interpreter
709-
710-
@staticmethod
711-
def surface_plot(matrix, ax=None, figsize=(8, 8), **kwargs):
712-
# acquire the cartesian coordinate matrices from the matrix
713-
# x is cols, y is rows
714-
(x, y) = np.meshgrid(np.arange(matrix.shape[0]), np.arange(matrix.shape[1]))
715-
716-
if ax is None:
717-
fig, ax = plt.subplots(figsize=figsize, subplot_kw=dict(projection='3d'))
718-
surf = ax.plot_surface(x, y, np.transpose(matrix), **kwargs)
719-
return ax, surf

0 commit comments

Comments
 (0)