Skip to content

Commit d78bb03

Browse files
gourav3017gourav3017
committed
Updated notebook using new package
1 parent f19c519 commit d78bb03

3 files changed

Lines changed: 252 additions & 126 deletions

File tree

‎Tutorial.ipynb‎

Lines changed: 142 additions & 68 deletions
Original file line numberDiff line numberDiff line change
@@ -1,15 +1,21 @@
11
{
22
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"id": "181b4f06",
6+
"metadata": {},
7+
"source": [
8+
"# <center> LowDimRT Tutorial </center>\n"
9+
]
10+
},
311
{
412
"cell_type": "code",
513
"execution_count": 1,
614
"id": "940f6eed",
715
"metadata": {},
816
"outputs": [],
917
"source": [
10-
"!pip install -r requirements.txt\n",
11-
"# import sys\n",
12-
"# sys.path.append('..')"
18+
"!pip install -r requirements.txt"
1319
]
1420
},
1521
{
@@ -23,8 +29,16 @@
2329
"from low_dim_rt import LowDimRT\n",
2430
"import numpy as np\n",
2531
"import cvxpy as cp\n",
26-
"import matplotlib.pyplot as plt\n",
27-
"# %run ex_wavelet.py"
32+
"import matplotlib.pyplot as plt"
33+
]
34+
},
35+
{
36+
"cell_type": "markdown",
37+
"id": "f77ebc50",
38+
"metadata": {},
39+
"source": [
40+
"### 1) Accessing the portpy data using DataExplorer class\n",
41+
" Note: you first need to download the patient database from the link provided in the GitHub page."
2842
]
2943
},
3044
{
@@ -44,45 +58,44 @@
4458
}
4559
],
4660
"source": [
47-
"# specify the patient data location\n",
48-
"# you first need to download the patient database from the link provided in the PortPy GitHub page\n",
49-
"data_dir = r'.\\data'\n",
50-
"# pick a patient from the existing patient list to get detailed info about the patient data (e.g., beams_dict, structures, ...)\n",
61+
"# specify the patient data location.\n",
62+
"data_dir = r'../data'\n",
63+
"\n",
64+
"# Use PortPy DataExplorer class to explore PortPy data and pick one of the patient\n",
65+
"data = pp.DataExplorer(data_dir=data_dir)\n",
5166
"patient_id = 'Lung_Patient_2'\n",
52-
"# create my_plan object for the planner beams_dict and select among the beams which are 30 degrees apart\n",
53-
"# for the customized beams_dict, you can pass the argument beam_ids\n",
54-
"my_plan = pp.Plan(patient_id, data_dir)\n"
67+
"data.patient_id = patient_id\n",
68+
"\n",
69+
"# Load ct and structure set for the above patient using CT and Structures class\n",
70+
"ct = pp.CT(data)\n",
71+
"structs = pp.Structures(data)\n",
72+
"\n",
73+
"# If the list of beams are not provided, it uses the beams selected manually\n",
74+
"# by a human expert planner for the patient (manually selected beams are stored in portpy data).\n",
75+
"# Create beams for the planner beams by default\n",
76+
"# for the customized beams, you can pass the argument beam_ids\n",
77+
"# e.g. beams = pp.Beams(data, beam_ids=[0,10,20,30,40,50,60])\n",
78+
"beams = pp.Beams(data)\n",
79+
"\n",
80+
"# create rinds based upon rind definition in optimization params\n",
81+
"protocol_name = 'Lung_2Gy_30Fx'\n",
82+
"opt_params = data.load_config_opt_params(protocol_name=protocol_name)\n",
83+
"structs.create_opt_structures(opt_params=opt_params)\n",
84+
"\n",
85+
"# load influence matrix based upon beams and structure set\n",
86+
"inf_matrix = pp.InfluenceMatrix(ct=ct, structs=structs, beams=beams)\n",
87+
"\n",
88+
"# load clinical criteria from the config files for which plan to be optimized\n",
89+
"clinical_criteria = pp.ClinicalCriteria(data, protocol_name=protocol_name)"
5590
]
5691
},
5792
{
58-
"cell_type": "code",
59-
"execution_count": 4,
60-
"id": "aaab565b",
93+
"cell_type": "markdown",
94+
"id": "27f0610c",
6195
"metadata": {},
62-
"outputs": [
63-
{
64-
"name": "stdout",
65-
"output_type": "stream",
66-
"text": [
67-
"creating rinds..\n",
68-
"rinds created!!\n"
69-
]
70-
}
71-
],
7296
"source": [
73-
"# Let us create rinds for creating reasonable dose fall off for the plan\n",
74-
"rind_max_dose = np.array([1.1, 1.05, 0.9, 0.85, 0.75]) * my_plan.get_prescription()\n",
75-
"rind_params = [{'rind_name': 'RIND_0', 'ref_structure': 'PTV', 'margin_start_mm': 0, 'margin_end_mm': 5,\n",
76-
" 'max_dose_gy': rind_max_dose[0]},\n",
77-
" {'rind_name': 'RIND_1', 'ref_structure': 'PTV', 'margin_start_mm': 5, 'margin_end_mm': 10,\n",
78-
" 'max_dose_gy': rind_max_dose[1]},\n",
79-
" {'rind_name': 'RIND_2', 'ref_structure': 'PTV', 'margin_start_mm': 10, 'margin_end_mm': 30,\n",
80-
" 'max_dose_gy': rind_max_dose[2]},\n",
81-
" {'rind_name': 'RIND_3', 'ref_structure': 'PTV', 'margin_start_mm': 30, 'margin_end_mm': 60,\n",
82-
" 'max_dose_gy': rind_max_dose[3]},\n",
83-
" {'rind_name': 'RIND_4', 'ref_structure': 'PTV', 'margin_start_mm': 60, 'margin_end_mm': 'inf',\n",
84-
" 'max_dose_gy': rind_max_dose[4]}]\n",
85-
"my_plan.add_rinds(rind_params=rind_params)"
97+
"### 2) Optimizing the plan without quadratic smoothness (With and without wavelet constraint)\n",
98+
"- Without wavelet constraint"
8699
]
87100
},
88101
{
@@ -103,10 +116,26 @@
103116
}
104117
],
105118
"source": [
106-
"# With no quadratic smoothness\n",
107-
"prob = pp.CvxPyProb(my_plan, smoothness_weight=0)\n",
108-
"prob.solve(solver='MOSEK', verbose=False)\n",
109-
"sol_no_quad_no_wav = prob.get_sol()"
119+
"# remove smoothness objective\n",
120+
"for i in range(len(opt_params['objective_functions'])):\n",
121+
" if opt_params['objective_functions'][i]['type'] == 'smoothness-quadratic':\n",
122+
" opt_params['objective_functions'][i]['weight'] = 0\n",
123+
"\n",
124+
"# create a plan using ct, structures, beams and influence matrix. Clinical criteria is optional\n",
125+
"my_plan = pp.Plan(ct, structs, beams, inf_matrix, clinical_criteria)\n",
126+
"\n",
127+
"# create cvxpy problem using the clinical criteria and optimization parameters\n",
128+
"opt = pp.Optimization(my_plan, opt_params=opt_params)\n",
129+
"opt.create_cvxpy_problem()\n",
130+
"sol_no_quad_no_wav = opt.solve(solver='MOSEK', verbose=False)"
131+
]
132+
},
133+
{
134+
"cell_type": "markdown",
135+
"id": "4ab09f28",
136+
"metadata": {},
137+
"source": [
138+
"- With wavelet constraint"
110139
]
111140
},
112141
{
@@ -120,9 +149,17 @@
120149
"wavelet_basis = LowDimRT.get_low_dim_basis(my_plan.inf_matrix, 'wavelet')\n",
121150
"# Smoothness Constraint\n",
122151
"y = cp.Variable(wavelet_basis.shape[1])\n",
123-
"prob.constraints += [wavelet_basis @ y == prob.x]\n",
124-
"prob.solve(solver='MOSEK', verbose=False)\n",
125-
"sol_no_quad_with_wav = prob.get_sol()"
152+
"opt.constraints += [wavelet_basis @ y == opt.vars['x']]\n",
153+
"sol_no_quad_with_wav = opt.solve(solver='MOSEK', verbose=False)"
154+
]
155+
},
156+
{
157+
"cell_type": "markdown",
158+
"id": "0342b80a",
159+
"metadata": {},
160+
"source": [
161+
"### 3) Optimizing with quadratic smoothness ( With and without wavelet constraint)\n",
162+
"- Without wavelet constraint"
126163
]
127164
},
128165
{
@@ -143,11 +180,25 @@
143180
}
144181
],
145182
"source": [
146-
"# create cvxpy problem using the clinical criteria\n",
147-
"prob = pp.CvxPyProb(my_plan, smoothness_weight=10)\n",
148-
"# run IMRT fluence map optimization using a low dimensional subspace for fluence map compression\n",
149-
"prob.solve(solver='MOSEK', verbose=False)\n",
150-
"sol_quad_no_wav = prob.get_sol()"
183+
"# set quadratic smoothness weight to 10\n",
184+
"for i in range(len(opt_params['objective_functions'])):\n",
185+
" if opt_params['objective_functions'][i]['type'] == 'smoothness-quadratic':\n",
186+
" opt_params['objective_functions'][i]['weight'] = 10\n",
187+
"\n",
188+
"# create cvxpy problem using the clinical criteria and optimization parameters\n",
189+
"opt = pp.Optimization(my_plan, opt_params=opt_params)\n",
190+
"opt.create_cvxpy_problem()\n",
191+
"\n",
192+
"# optimize the plan\n",
193+
"sol_quad_no_wav = opt.solve(solver='MOSEK', verbose=False)"
194+
]
195+
},
196+
{
197+
"cell_type": "markdown",
198+
"id": "1a584044",
199+
"metadata": {},
200+
"source": [
201+
"- With Wavelet constraint"
151202
]
152203
},
153204
{
@@ -157,11 +208,18 @@
157208
"metadata": {},
158209
"outputs": [],
159210
"source": [
160-
"# Smoothness Constraint\n",
211+
"# Wavelet Smoothness Constraint\n",
161212
"y = cp.Variable(wavelet_basis.shape[1])\n",
162-
"prob.constraints += [wavelet_basis @ y == prob.x]\n",
163-
"prob.solve(solver='MOSEK', verbose=False)\n",
164-
"sol_quad_with_wav = prob.get_sol()"
213+
"opt.constraints += [wavelet_basis @ y == opt.vars['x']]\n",
214+
"sol_quad_with_wav = opt.solve(solver='MOSEK', verbose=False)"
215+
]
216+
},
217+
{
218+
"cell_type": "markdown",
219+
"id": "c403b30e",
220+
"metadata": {},
221+
"source": [
222+
"### 4) Saving and loading plan"
165223
]
166224
},
167225
{
@@ -171,11 +229,11 @@
171229
"metadata": {},
172230
"outputs": [],
173231
"source": [
174-
"pp.save_plan(my_plan, plan_name='my_plan', path=r'C:\\temp')\n",
175-
"pp.save_optimal_sol(sol_no_quad_no_wav, sol_name='sol_no_quad_no_wav', path=r'C:\\temp')\n",
176-
"pp.save_optimal_sol(sol_no_quad_with_wav, sol_name='sol_no_quad_with_wav', path=r'C:\\temp')\n",
177-
"pp.save_optimal_sol(sol_quad_no_wav, sol_name='sol_quad_no_wav', path=r'C:\\temp')\n",
178-
"pp.save_optimal_sol(sol_quad_with_wav, sol_name='sol_quad_with_wav', path=r'C:\\temp')"
232+
"pp.save_plan(my_plan, plan_name='my_plan.pkl', path=r'C:\\temp')\n",
233+
"pp.save_optimal_sol(sol_no_quad_no_wav, sol_name='sol_no_quad_no_wav.pkl', path=r'C:\\temp')\n",
234+
"pp.save_optimal_sol(sol_no_quad_with_wav, sol_name='sol_no_quad_with_wav.pkl', path=r'C:\\temp')\n",
235+
"pp.save_optimal_sol(sol_quad_no_wav, sol_name='sol_quad_no_wav.pkl', path=r'C:\\temp')\n",
236+
"pp.save_optimal_sol(sol_quad_with_wav, sol_name='sol_quad_with_wav.pkl', path=r'C:\\temp')"
179237
]
180238
},
181239
{
@@ -192,6 +250,14 @@
192250
"# sol_quad_with_wav = pp.load_optimal_sol(sol_name='sol_quad_with_wav', path=r'C:\\temp')"
193251
]
194252
},
253+
{
254+
"cell_type": "markdown",
255+
"id": "4650e367",
256+
"metadata": {},
257+
"source": [
258+
"### 5) Visualize and compare the plans"
259+
]
260+
},
195261
{
196262
"cell_type": "code",
197263
"execution_count": 10,
@@ -223,13 +289,13 @@
223289
"# plot fluence 3D and 2D\n",
224290
"fig, ax = plt.subplots(1, 2, figsize=(18,6), subplot_kw={'projection': '3d'})\n",
225291
"fig.suptitle('Without Quadratic smoothness')\n",
226-
"pp.Visualize.plot_fluence_3d(sol=sol_no_quad_no_wav, beam_id=37, ax=ax[0], title='Without Wavelet')\n",
227-
"pp.Visualize.plot_fluence_3d(sol=sol_no_quad_with_wav, beam_id=37, ax=ax[1], title='With Wavelet')\n",
292+
"pp.Visualization.plot_fluence_3d(sol=sol_no_quad_no_wav, beam_id=37, ax=ax[0], title='Without Wavelet')\n",
293+
"pp.Visualization.plot_fluence_3d(sol=sol_no_quad_with_wav, beam_id=37, ax=ax[1], title='With Wavelet')\n",
228294
"\n",
229295
"fig, ax = plt.subplots(1, 2, figsize=(18,6), subplot_kw={'projection': '3d'})\n",
230296
"fig.suptitle('With Quadratic smoothness')\n",
231-
"pp.Visualize.plot_fluence_3d(sol=sol_quad_no_wav, beam_id=37, ax=ax[0], title='Without Wavelet')\n",
232-
"pp.Visualize.plot_fluence_3d(sol=sol_quad_with_wav, beam_id=37, ax=ax[1], title='With Wavelet')\n",
297+
"pp.Visualization.plot_fluence_3d(sol=sol_quad_no_wav, beam_id=37, ax=ax[0], title='Without Wavelet')\n",
298+
"pp.Visualization.plot_fluence_3d(sol=sol_quad_with_wav, beam_id=37, ax=ax[1], title='With Wavelet')\n",
233299
"\n",
234300
"plt.show()"
235301
]
@@ -255,20 +321,28 @@
255321
"# plot DVH for the structures in the given list. Default dose_1d is in Gy and volume is in relative scale(%).\n",
256322
"structs = ['PTV', 'ESOPHAGUS', 'HEART', 'CORD', 'LUNG_L', 'LUNG_R']\n",
257323
"fig, ax = plt.subplots(1, 2, figsize=(20, 8))\n",
258-
"ax0 = pp.Visualize.plot_dvh(my_plan, sol=sol_no_quad_no_wav, structs=structs, style='solid', ax=ax[0])\n",
259-
"ax0 = pp.Visualize.plot_dvh(my_plan, sol=sol_no_quad_with_wav, structs=structs, style='dashed', ax=ax0)\n",
324+
"ax0 = pp.Visualization.plot_dvh(my_plan, sol=sol_no_quad_no_wav, structs=structs, style='solid', ax=ax[0])\n",
325+
"ax0 = pp.Visualization.plot_dvh(my_plan, sol=sol_no_quad_with_wav, structs=structs, style='dashed', ax=ax0)\n",
260326
"fig.suptitle('DVH comparison')\n",
261327
"ax0.set_title('Without Quadratic smoothness \\n solid: Without Wavelet, Dash: With Wavelet')\n",
262328
"# plt.show()\n",
263329
"# print('\\n\\n')\n",
264330
"\n",
265331
"# fig, ax = plt.subplots(figsize=(12, 8))\n",
266-
"ax1 = pp.Visualize.plot_dvh(my_plan, sol=sol_quad_no_wav, structs=structs, style='solid', ax=ax[1])\n",
267-
"ax1 = pp.Visualize.plot_dvh(my_plan, sol=sol_quad_with_wav, structs=structs, style='dashed', ax=ax1)\n",
332+
"ax1 = pp.Visualization.plot_dvh(my_plan, sol=sol_quad_no_wav, structs=structs, style='solid', ax=ax[1])\n",
333+
"ax1 = pp.Visualization.plot_dvh(my_plan, sol=sol_quad_with_wav, structs=structs, style='dashed', ax=ax1)\n",
268334
"ax1.set_title('With Quadratic smoothness \\n solid: Without Wavelet, Dash: With Wavelet')\n",
269335
"plt.show()\n"
270336
]
271337
},
338+
{
339+
"cell_type": "markdown",
340+
"id": "f154abb6",
341+
"metadata": {},
342+
"source": [
343+
"### 6) Evaluate the plans"
344+
]
345+
},
272346
{
273347
"cell_type": "code",
274348
"execution_count": 13,
@@ -500,7 +574,7 @@
500574
],
501575
"source": [
502576
"# visualize plan metrics based upon clinical criteria\n",
503-
"pp.Visualize.plan_metrics(my_plan, sol=[sol_no_quad_no_wav, sol_no_quad_with_wav, sol_quad_no_wav, sol_quad_with_wav], sol_names=['no_quad_no_wav', 'no_quad_with_wav', 'quad_no_wav', 'quad_with_wav'])"
577+
"pp.Evaluation.plan_metrics(my_plan, sol=[sol_no_quad_no_wav, sol_no_quad_with_wav, sol_quad_no_wav, sol_quad_with_wav], sol_names=['no_quad_no_wav', 'no_quad_with_wav', 'quad_no_wav', 'quad_with_wav'])"
504578
]
505579
}
506580
],

0 commit comments

Comments
 (0)