Skip to content

Commit 6313423

Browse files
msneubauermsneubauer
committed
updates from student author
1 parent 30f9625 commit 6313423

1 file changed

Lines changed: 93 additions & 64 deletions

File tree

‎_sources/projects/Project_02_Radiotherapy.ipynb‎

Lines changed: 93 additions & 64 deletions
Original file line numberDiff line numberDiff line change
@@ -164,7 +164,7 @@
164164
"Your tasks:\n",
165165
"1) Add a line to the for loop that stacks `ct`, `ptv`, and `oars` (Organs at Risk) into a single tensor. Note that oars is a dictionary you may need to unpack. What's the final shape of this tensor? Draw an analogy to an RGB image. \n",
166166
"\n",
167-
"2) Visualize a single 2D axial slice (z-axis slice) of a patient's CT scan with the tumor(ptv) overlaid...as a sanity check."
167+
"2) Visualize a single 2D slice (axial, sagittal, or coronal) of a patient's CT scan with the tumor(ptv) overlaid...as a sanity check."
168168
]
169169
},
170170
{
@@ -241,7 +241,7 @@
241241
"metadata": {},
242242
"outputs": [],
243243
"source": [
244-
"#Your slice visualization here\n",
244+
"#Your slice visualization here: \n",
245245
"import matplotlib.pyplot as plt"
246246
]
247247
},
@@ -258,13 +258,13 @@
258258
"source": [
259259
"To simulate the effects of a naive beam (simple straight lines) passing through the patient, we provide a 'beam mask' that creates a beam onto a 2D slice of the patient scan. We assume there are 36 possible angles from which a beam can be fired (one every 10 degrees). `dose` is a 2D heatmap that represents how much radiation is delivered to each pixel in your slice after the passage of many beams.\n",
260260
"\n",
261-
"Create a Reinforcement Learning Model which selects beam angles that provide radiation dosage to the tumor while avoiding the spinal cord. \n",
261+
"Create a Reinforcement Learning Model which selects beam angles that provide radiation dosage to the tumor while avoiding critical organs. \n",
262262
"\n",
263263
"* Populate the training loop by: \n",
264264
" * Creating a policy that determines the next `angle` \n",
265265
" * Creating a Q-learning update \n",
266266
"\n",
267-
"* Define a reward function that heavily penalizes a beam that passes through the spinal cord and rewards a beam that hits the tumor/ptv (the output should be a score) \n",
267+
"* Define a reward function that heavily penalizes a beam that passes through an organ and rewards a beam that hits the tumor/ptv (the output should be a score) \n",
268268
"* Tweak parameters if you so desire\n",
269269
"* Print the top 3 beam angles after training is complete"
270270
]
@@ -299,35 +299,37 @@
299299
" end = (center + length * np.array([dy, dx])).astype(int)\n",
300300
" start = np.clip(start, 0, [h - 1, w - 1])\n",
301301
" end = np.clip(end, 0, [h - 1, w - 1])\n",
302-
" rr, cc = line_nd(start[0], start[1], end[0], end[1])\n",
302+
" rr, cc = line_nd(start, end)\n",
303303
" beam_mask[rr, cc] = 1\n",
304304
" return beam_mask\n",
305305
"\n",
306306
"\n",
307-
"for episode in range(episodes):\n",
308-
" #Randomly select a patient from patient_tensors HERE\n",
307+
" for episode in range(episodes):\n",
308+
" #Randomly select a patient from patient_tensors HERE\n",
309309
"\n",
310-
" #extracts a single 2D slice from the middle\n",
311-
" mid = ct.shape[2] // 2 \n",
312-
" ct_slice = ct[:, :, mid]\n",
313-
" ptv_slice = ptv[:, :, mid]\n",
314-
" #extract organ slices HERE\n",
310+
" ct,ptv =\n",
315311
"\n",
316-
" dose = np.zeros(slice_shape, dtype=np.float32)\n",
317-
" selected_angles = []\n",
312+
" #extracts a single 2D slice from the middle\n",
313+
" mid = ct.shape[2] // 2 \n",
314+
" ct_slice = ct[:, :, mid]\n",
315+
" ptv_slice = ptv[:, :, mid]\n",
316+
" #extract organ slices HERE\n",
318317
"\n",
319-
" for i in range(max_beams):\n",
320-
" #Your policy here\n",
318+
" dose = np.zeros(slice_shape, dtype=np.float32)\n",
319+
" selected_angles = []\n",
321320
"\n",
322-
" beam = generate_beam(angle, slice_shape)\n",
323-
" dose += beam.astype(np.float32) #adds 'radiation' to the pixels \n",
324-
" selected_angles.append(angle)\n",
321+
" for i in range(max_beams):\n",
322+
" #Your policy here\n",
325323
"\n",
326-
" \n",
327-
" reward_score = reward(dose, ptv, organs)\n",
324+
" beam = generate_beam(angle, slice_shape)\n",
325+
" dose += beam.astype(np.float32) #adds 'radiation' to the pixels \n",
326+
" selected_angles.append(angle)\n",
327+
"\n",
328+
" \n",
329+
" reward_score = reward(dose, ptv, organs)\n",
328330
"\n",
329-
" for angle in selected_angles:\n",
330-
" #Q-learning update here\n",
331+
" for angle in selected_angles:\n",
332+
" #Q-learning update here\n",
331333
"\n",
332334
"\n",
333335
"\n",
@@ -351,12 +353,9 @@
351353
"metadata": {},
352354
"outputs": [],
353355
"source": [
354-
"import matplotlib as plt\n",
355-
"\n",
356-
"\n",
357356
"patient = random.choice(patient_tensors)\n",
358357
"mid = patient.shape[3] // 2\n",
359-
"ct, ptv, cord = patient[0, :, :, mid], patient[1, :, :, mid], patient[2, :, :, mid]\n",
358+
"ct, ptv = patient[0, :, :, mid], patient[1, :, :, mid]\n",
360359
"shape = ct.shape\n",
361360
"dose = np.zeros(shape, dtype=np.float32)\n",
362361
"top_indices = sorted(range(len(Q)), key=lambda i: -Q[i])[:max_beams]\n",
@@ -388,7 +387,7 @@
388387
"cell_type": "markdown",
389388
"metadata": {},
390389
"source": [
391-
"Ready to boogie? We're about to imbue physics into our RL model. Let's start by making our beam of ionizing radiaiton more realistic using what we discussed in Question 01. Adjust the `decay` in `generate_true_beam` below so that it actually attenuates our beam (energy loss as it passes through tissue)"
390+
"Ready to boogie? We're about to imbue physics into our RL model. Let's start by making our beam of ionizing radiation more realistic using what we discussed in Question 01. Adjust the `decay` in `generate_true_beam` below so that it actually attenuates our beam (energy loss as it passes through tissue)"
392391
]
393392
},
394393
{
@@ -397,10 +396,17 @@
397396
"metadata": {},
398397
"outputs": [],
399398
"source": [
400-
"def generate_physical_beam(angle_deg, shape = (128,128), spread_sigma=2.0, attenuation_coeff=0.01):\n",
399+
"from scipy.ndimage import center_of_mass\n",
400+
"\n",
401+
"def generate_physical_beam(angle_deg, ptv_mask, shape = (128,128), spread_sigma=2.0, attenuation_coeff=0.01):\n",
401402
" h, w = shape\n",
402403
" beam_mask = np.zeros((h, w), dtype=np.float32)\n",
403-
" center = np.array([h // 2, w // 2])\n",
404+
" if isinstance(ptv_mask, torch.Tensor):\n",
405+
" ptv_mask_np = ptv_mask.detach().cpu().numpy()\n",
406+
" else:\n",
407+
" ptv_mask_np = ptv_mask\n",
408+
" tumor_center = np.array(center_of_mass(ptv_mask_np))\n",
409+
" center = tumor_center.astype(int)\n",
404410
" length = max(h, w)\n",
405411
" angle_rad = np.deg2rad(angle_deg)\n",
406412
"\n",
@@ -416,7 +422,7 @@
416422
" cc = np.clip(cc, 0, w-1)\n",
417423
"\n",
418424
" for i, (r, c) in enumerate(zip(rr, cc)):\n",
419-
" decay = #attenuatation by Beer-Lambert Law\n",
425+
" decay = 1*np.exp(-attenuation_coeff*i)\n",
420426
" beam_mask[r, c] += decay\n",
421427
"\n",
422428
" for dr in range(-3, 4): #gaussian spread via scattering\n",
@@ -426,9 +432,9 @@
426432
" if 0 <= r_spread < h and 0 <= c_spread < w:\n",
427433
" distance = np.sqrt(dr**2 + dc**2)\n",
428434
" spread_value = np.exp(- (distance**2) / (2 * spread_sigma**2))\n",
429-
" beam_mask[r_spread, c_spread] += decay * spread_value \n",
435+
" beam_mask[r_spread, c_spread] += decay * spread_value\n",
430436
"\n",
431-
" beam_mask = np.clip(beam_mask, 0, 1.0) \n",
437+
" beam_mask = np.clip(beam_mask, 0, 1.0)\n",
432438
" return beam_mask"
433439
]
434440
},
@@ -445,7 +451,7 @@
445451
"$$\n",
446452
"L_{\\text{Physics}} =\n",
447453
"\\lambda_{\\text{PTV}} \\cdot \\text{Underdose}_{\\text{PTV}} +\n",
448-
"\\lambda_{\\text{Spine}} \\cdot \\text{Overdose}_{\\text{Spine}}\n",
454+
"\\lambda_{\\text{Spine}} \\cdot \\text{Overdose}_{\\text{Organs}}\n",
449455
"$$\n",
450456
"Each $\\lambda$ assigns 'weight' to a term, tinkering with these may prove beneficial.\n",
451457
"The Underdose term penalizes not giving the tumor the maximum dose (1)\n",
@@ -458,12 +464,12 @@
458464
"- $\\max(0, 1 - D_{\\text{PTV}})$: penalizes only underdosed regions\n",
459465
"- $\\mathbb{E}[\\cdot]$: average \n",
460466
"\n",
461-
"The Overdose term penalizes giving the spine too much dose\n",
467+
"The Overdose term penalizes giving the organs too much dose\n",
462468
"$$\n",
463-
"\\text{Overdose}_{\\text{Spine}} = \\mathbb{E} \\left[ \\max \\left( 0,\\ D_{\\text{Spine}} - D_{\\text{safe}} \\right) \\right]\n",
469+
"\\text{Overdose}_{\\text{Organs}} = \\mathbb{E} \\left[ \\max \\left( 0,\\ D_{\\text{Organs}} - D_{\\text{safe}} \\right) \\right]\n",
464470
"$$\n",
465471
"\n",
466-
"- $D_{\\text{safe}}$: Safety threshold for organ (max dose it can receive). **In clinical settings, this is often 0.6**\n",
472+
"- $D_{\\text{safe}}$: Safety threshold for organ (max dose it can receive). **You can find these online**\n",
467473
"\n",
468474
"**Now create that network, like a true ML-infused medical physicist of the 21st century would!**"
469475
]
@@ -490,41 +496,64 @@
490496
"metadata": {},
491497
"outputs": [],
492498
"source": [
493-
"#your conclusions\n",
494-
"selected_angles = \n",
499+
"from matplotlib.lines import Line2D\n",
500+
"def visualize_beams_on_ct(patient, selected_angles, slice_index= None, shape = (128,128)):\n",
501+
" \n",
502+
" channel, height, width, depth = patient.shape\n",
495503
"\n",
496-
"# Choose a patient!\n",
497-
"patient = patient_tensors[0]\n",
504+
" if slice_index is None:\n",
505+
" slice_index = depth // 2 #default to middle slice\n",
498506
"\n",
507+
" ct_slice = patient[0, :, :, slice_index]\n",
508+
" ptv_slice = patient[1, :, :, slice_index]\n",
509+
" organ_slices = patient[2:, :, :, slice_index]\n",
499510
"\n",
511+
" dose = np.zeros(shape, dtype=np.float32)\n",
512+
" for angle in selected_angles:\n",
513+
" beam = generate_physical_beam(angle, ptv_slice, shape)\n",
514+
" dose += beam.astype(np.float32)\n",
515+
" dose_norm = dose / np.max(dose + 1e-5)\n",
500516
"\n",
501-
"ct_volume = patient[\"ct\"]\n",
502-
"ptv_volume = patient[\"ptv\"]\n",
503-
"spine_volume = patient[\"spine\"]\n",
504-
"mid = ct_volume.shape[2] // 2\n",
505-
"ct_slice = ct_volume[:, :, mid]\n",
506-
"ptv_slice = ptv_volume[:, :, mid]\n",
507-
"spine_slice = spine_volume[:, :, mid]\n",
517+
" plt.figure(figsize=(6, 6))\n",
518+
" plt.imshow(ct_slice, cmap='gray', alpha=0.6)\n",
519+
" plt.imshow(dose_norm, cmap='hot', alpha=0.4) #beam\n",
508520
"\n",
521+
" plt.contour(ptv_slice, colors='red', linewidths=1, label='PTV')\n",
522+
" for organ_slice in organ_slices:\n",
523+
" plt.contour(organ_slice, colors='blue', linewidths=1)\n",
524+
" \n",
525+
" plt.title(\"Radiotherapy Treatment Plan\")\n",
526+
" plt.axis('off')\n",
509527
"\n",
510-
"def visualize_beams_on_ct(ct_slice, ptv_slice, spine_slice, selected_angles, shape=(128, 128)):\n",
511-
" dose = np.zeros(shape, dtype=np.float32)\n",
512-
" for angle in selected_angles:\n",
513-
" beam = generate_physical_beam(angle, shape)\n",
514-
" dose += beam.astype(np.float32)\n",
515-
" dose_norm = dose / np.max(dose + 1e-5)\n",
516-
" plt.figure(figsize=(6, 6))\n",
517-
" plt.imshow(ct_slice, cmap='gray', alpha=0.6)\n",
518-
" plt.imshow(dose_norm, cmap='hot', alpha=0.4) #beam\n",
519-
" plt.contour(ptv_slice, colors='red', linewidths=1, label='PTV')\n",
520-
" plt.contour(spine_slice, colors='blue', linewidths=1, label='Spine')\n",
521-
" plt.title(\"Treatment Plan\")\n",
522-
" plt.axis('off')\n",
523-
" plt.colorbar(label=\"Relative Beam Intensity\")\n",
524-
" plt.show()\n",
528+
" legend_elements = [\n",
529+
" Line2D([0], [0], color='red', lw=2, label='Tumor'),\n",
530+
" Line2D([0], [0], color='blue', lw=2, label='Organs at risk')\n",
531+
" ]\n",
525532
"\n",
533+
" plt.legend(handles=legend_elements, loc='upper right')\n",
534+
" plt.show()\n"
535+
]
536+
},
537+
{
538+
"cell_type": "markdown",
539+
"metadata": {},
540+
"source": [
541+
"**(Optional adventures)** \n",
526542
"\n",
527-
"visualize_beams_on_ct(ct_slice, ptv_slice, spine_slice, selected_angles)\n"
543+
"If you'd like to spice things up and maybe even contribute to current medical physics research, here's some fun ideas:\n",
544+
"- The model & beam currently work on 2D slices, give it 3D awareness! \n",
545+
"- Imbibe more physics into the beam (e.g. consider penalizing the beam for deviating from real scattering physics in your reward function)\n",
546+
"- Incorporate consideration of the varying physical properties of tissue for each organ/bone/tumor region\n",
547+
"- Try finding datasets of real beam angles used in clinical settings and use a supervised learning approach to compare results.\n",
548+
"- Add more patients from OpenKBP and train for a long, long time. Got any supercomputers handy? \n",
549+
"\n"
550+
]
551+
},
552+
{
553+
"cell_type": "markdown",
554+
"metadata": {},
555+
"source": [
556+
"Thanks for doing my project! Feel free to reach out with any questions/comments on the notebook or beyond!! :)"
528557
]
529558
},
530559
{

0 commit comments

Comments
 (0)