DICOM RTPLAN — MLC Fluence 互動式教學¶
這份 notebook 帶你一步步看懂放射治療計畫檔（RTPLAN）：
- 讀取 DICOM 檔案
- 觀察 Beam 與 Control Point（CP）的結構
- 把 MLC / Jaw 位置畫成二維的「照野」
- 依照每個 CP 的 MU 權重，累積出 fluence 分佈圖
- 用滑桿互動瀏覽每一個 CP 的射束形狀
需要的套件
pip install ipympl pydicom numpy matplotlib
Step 1｜匯入套件¶
%matplotlib widget 讓 matplotlib 畫的圖可以互動（滑桿、按鈕即時回應）。
In [ ]:
!pip install ipympl pydicom numpy matplotlib #安裝套件
In [ ]:
%matplotlib widget
import pydicom as dcm
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import transforms
from matplotlib.widgets import Slider, RadioButtons
Step 2｜讀取 DICOM 計畫檔¶
請把 .dcm 檔放在 notebook 同一個資料夾，或修改下面的路徑。
In [ ]:
DICOM_PATH = './neurp.dcm' # ← 改成你的檔案路徑
plan = dcm.dcmread(DICOM_PATH)
print(f"Modality : {plan.Modality}")
print(f"Patient ID : {getattr(plan, 'PatientID', 'N/A')}")
print(f"Beam 總數 : {len(plan.BeamSequence)}")
Step 3｜認識資料結構¶
一個治療計畫包含多支 Beam，每支 Beam 又由多個 Control Point（CP） 組成。
BeamMeterset：這支 Beam 總共要打多少 MU（機器單位）CumulativeMetersetWeight：0 到 1 的累積比例，兩個 CP 相減 = 該段的相對 MUBeamLimitingDevicePositionSequence：記錄當下 CP 的 MLC / Jaw 位置
💡 重要概念：DICOM 允許某個 CP 不重新定義 MLC 或 Jaw， 這時就代表「沿用上一個 CP」的設定。這是後面「位置繼承」的依據。
In [ ]:
for i, beam in enumerate(plan.BeamSequence):
meterset = plan.FractionGroupSequence[0].ReferencedBeamSequence[i].BeamMeterset
n_cp = len(beam.ControlPointSequence)
print(f"[{i}] {beam.BeamName:<24} CP={n_cp:<3} MU={float(meterset):.1f}")
In [ ]:
first_beam = plan.BeamSequence[0]
print(f"\nBeam 0「{first_beam.BeamName}」的前 5 個 Control Point：\n")
for i, cp in enumerate(first_beam.ControlPointSequence[:5]):
defined = [d.RTBeamLimitingDeviceType
for d in getattr(cp, 'BeamLimitingDevicePositionSequence', [])]
print(f" CP{i}: 累積MU比例={cp.CumulativeMetersetWeight:.3f} "
f"此CP定義的裝置={defined}")
4-1. 定義常數與轉換函式¶
In [ ]:
IMAGE_SIZE = 400 # 陣列邊長（像素）
BEV_MIN_MM = -200.0 # BEV 左邊界（mm）
BEV_MAX_MM = 200.0 # BEV 右邊界（mm）
MM_PER_PIXEL = (BEV_MAX_MM - BEV_MIN_MM) / IMAGE_SIZE
EXTENT = [BEV_MIN_MM, BEV_MAX_MM, BEV_MIN_MM, BEV_MAX_MM]
def mm_to_pixel(mm):
"""把 mm 座標轉成像素索引（四捨五入到最近的整數）。"""
return int(round((mm - BEV_MIN_MM) / MM_PER_PIXEL))
def clamp_pixel(index):
"""把像素索引限制在 [0, IMAGE_SIZE] 之間，避免陣列越界。"""
return max(0, min(IMAGE_SIZE, index))
print(f"解析度：1 pixel = {MM_PER_PIXEL} mm")
print(f"BEV 範圍：[{BEV_MIN_MM}, {BEV_MAX_MM}] mm")
4-2. 用一個小範例驗證 Y 軸方向¶
下面的 3×3 陣列，配上 origin='lower' 後，[0, 0] 會畫在左下角。
這樣陣列的 第一個 index 就是 Y、第二個 index 就是 X，
跟我們對 BEV 平面的直覺完全一致。
In [ ]:
demo = np.array([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
fig_demo, demo_ax = plt.subplots(figsize=(4, 4))
demo_ax.imshow(demo, cmap='YlOrRd', extent=[0, 3, 0, 3], origin='lower')
demo_ax.set_title("demo[0,0] 在左下角")
demo_ax.set_xlabel("X"); demo_ax.set_ylabel("Y")
for (i, j), val in np.ndenumerate(demo):
demo_ax.text(j + 0.5, i + 0.5, str(val),
ha='center', va='center', fontsize=14)
plt.show()
Step 5｜核心函式：建立每支 Beam 的 Fluence¶
這個函式對每個 Control Point 做四件事：
- 算這段的 MU — 用累積比例相減
- 畫出 MLC 遮罩 — 把每片葉子的位置填成布林陣列
- 畫出 Jaw 遮罩 — Jaw 是矩形，把四個邊界外的區域遮住
- 兩者相乘 × 段 MU — 得到這段打開的區域，再加到累積圖上
函式回傳 snapshot 列表，每個元素是「到這個 CP 為止」的累積 fluence 與當時的 MLC / Jaw 位置。
In [ ]:
def build_fluence_for_beam(beam_index):
"""為指定 Beam 建立每個 Control Point 的累積 fluence snapshot。"""
beam = plan.BeamSequence[beam_index]
beam_meterset = plan.FractionGroupSequence[0].ReferencedBeamSequence[beam_index].BeamMeterset
# --- 每一段的 MU 權重 ---
cumulative_weights = [cp.CumulativeMetersetWeight for cp in beam.ControlPointSequence]
n_cp = len(cumulative_weights)
segment_mu = [cumulative_weights[i] - cumulative_weights[i - 1] if i > 0
else cumulative_weights[i]
for i in range(n_cp)]
# --- 找出 MLC 裝置與葉片 Y 邊界 ---
leaf_boundaries = None
for device in beam.BeamLimitingDeviceSequence:
if device.RTBeamLimitingDeviceType in ('MLCX', 'MLC') \
and hasattr(device, 'LeafPositionBoundaries'):
leaf_boundaries = device.LeafPositionBoundaries
break
# --- 跨 CP 繼承用的狀態變數 ---
current_gantry_angle = 0.0
current_collimator_angle = 0.0
current_jaw_x = None
current_jaw_y = None
fluence_cumulative = np.zeros((IMAGE_SIZE, IMAGE_SIZE))
snapshots = []
for cp_index, cp in enumerate(beam.ControlPointSequence):
# ===== 1. 角度繼承 =====
if hasattr(cp, 'GantryAngle'):
current_gantry_angle = float(cp.GantryAngle)
if hasattr(cp, 'BeamLimitingDeviceAngle'):
current_collimator_angle = float(cp.BeamLimitingDeviceAngle)
elif hasattr(cp, 'CollimatorAngle'):
current_collimator_angle = float(cp.CollimatorAngle)
# ===== 2. 讀這個 CP 有定義的裝置位置 =====
mlc_raw = None
for device_position in getattr(cp, 'BeamLimitingDevicePositionSequence', []):
kind = device_position.RTBeamLimitingDeviceType
if kind in ('ASYMX', 'X'):
current_jaw_x = device_position.LeafJawPositions
elif kind in ('ASYMY', 'Y'):
current_jaw_y = device_position.LeafJawPositions
elif kind in ('MLCX', 'MLC'):
mlc_raw = device_position.LeafJawPositions
# ===== 3. 把 MLC 位置畫成遮罩 =====
blocked_by_mlc = np.zeros((IMAGE_SIZE, IMAGE_SIZE), dtype=bool)
mlc_array = np.array([])
expected_size = 0
if mlc_raw is not None and leaf_boundaries is not None:
mlc_array = np.array(mlc_raw)
n_leaves = len(leaf_boundaries) - 1
expected_size = n_leaves * 2
if mlc_array.size == expected_size:
for leaf in range(n_leaves):
row_start = clamp_pixel(mm_to_pixel(leaf_boundaries[leaf]))
row_end = clamp_pixel(mm_to_pixel(leaf_boundaries[leaf + 1]))
left_x = clamp_pixel(mm_to_pixel(mlc_array[leaf]))
right_x = clamp_pixel(mm_to_pixel(mlc_array[n_leaves + leaf]))
blocked_by_mlc[row_start:row_end, 0:left_x] = True
blocked_by_mlc[row_start:row_end, right_x:IMAGE_SIZE] = True
# ===== 4. 把 Jaw 位置畫成遮罩 =====
blocked_by_jaw = np.zeros((IMAGE_SIZE, IMAGE_SIZE), dtype=bool)
if current_jaw_x is not None:
x_left = clamp_pixel(mm_to_pixel(current_jaw_x[0]))
x_right = clamp_pixel(mm_to_pixel(current_jaw_x[1]))
blocked_by_jaw[:, 0:x_left] = True
blocked_by_jaw[:, x_right:IMAGE_SIZE] = True
if current_jaw_y is not None:
y_bottom = clamp_pixel(mm_to_pixel(current_jaw_y[0]))
y_top = clamp_pixel(mm_to_pixel(current_jaw_y[1]))
blocked_by_jaw[0:y_bottom, :] = True
blocked_by_jaw[y_top:IMAGE_SIZE, :] = True
# ===== 5. 計算這段 fluence 並累加 =====
open_field = ~blocked_by_mlc & ~blocked_by_jaw
segment_fluence = open_field * segment_mu[cp_index] * beam_meterset
fluence_cumulative += segment_fluence
# ===== 6. 記錄 snapshot =====
mlc_positions = None
if mlc_array.size == expected_size and expected_size > 0:
mlc_positions = mlc_array.reshape(2, expected_size // 2).tolist()
snapshots.append({
'fluence_cumulative': np.copy(fluence_cumulative),
'jaw_x': current_jaw_x,
'jaw_y': current_jaw_y,
'mlc_positions': mlc_positions,
'leaf_boundaries': leaf_boundaries,
'cumulative_mu_percent': cumulative_weights[cp_index] * 100,
'gantry_angle': current_gantry_angle,
'collimator_angle': current_collimator_angle,
})
return snapshots, beam_meterset
Step 6｜建立互動介面¶
這裡我們建立三個東西：
- 主圖：顯示 fluence 與 MLC / Jaw 輪廓
- Beam 選單（左側）：切換不同射束
- CP 滑桿（下方）：瀏覽單一射束內的所有 Control Point
6-1. 全域狀態變數¶
widget 的 callback 之間用全域變數共享狀態，這是 matplotlib 的標準做法。
In [ ]:
figure = None
axes = None
slider_axes = None
fluence_image = None
cp_slider = None
snapshots = []
beam_meterset = 0.0
beam_labels = []
current_beam = 0
outline_artists = []
6-2. 繪製單一 Control Point¶
每當滑桿改變時，重畫三樣東西：
- Fluence 彩色圖：用
set_data直接換內容，比重新建立快很多 - 紅框：Jaw 的邊界
- 黑框：MLC 每一片葉子的位置
旋轉角度來自 CollimatorAngle，用 rotate_deg_around(0, 0, -angle) 讓整個
圖形繞著 BEV 原點旋轉。
In [ ]:
def draw_control_point(_value):
"""當滑桿數值改變時，重畫當前 Control Point。"""
if not snapshots:
return
cp_index = int(cp_slider.val)
cp_index = max(0, min(cp_index, len(snapshots) - 1))
snap = snapshots[cp_index]
# 旋轉 + 資料座標 的合成轉換
rotation = transforms.Affine2D().rotate_deg_around(0, 0, -snap['collimator_angle'])
transform = rotation + axes.transData
# ===== Fluence 影像（原地更新） =====
fluence = snap['fluence_cumulative']
masked = np.ma.masked_where(fluence == 0.0, fluence)
fluence_image.set_data(masked)
fluence_image.set_clim(0, max(beam_meterset * 0.7, 1e-6))
fluence_image.set_transform(transform)
# ===== 清除舊輪廓 =====
for artist in outline_artists:
artist.remove()
outline_artists.clear()
# ===== 畫 Jaw（紅框） =====
if snap['jaw_x'] is not None and snap['jaw_y'] is not None:
x1, x2 = snap['jaw_x']
y1, y2 = snap['jaw_y']
xs = [x1, x2, x2, x1, x1]
ys = [y1, y1, y2, y2, y1]
line, = axes.plot(xs, ys, color='red', linewidth=3, transform=transform)
outline_artists.append(line)
# ===== 畫 MLC（黑框，一片一片畫） =====
mlc = snap['mlc_positions']
boundaries = snap['leaf_boundaries']
if mlc is not None and boundaries is not None:
left_bank, right_bank = mlc[0], mlc[1]
for leaf in range(len(left_bank)):
y0, y1 = boundaries[leaf], boundaries[leaf + 1]
xl, xr = left_bank[leaf], right_bank[leaf]
xs = [xl, xr, xr, xl, xl]
ys = [y0, y0, y1, y1, y0]
line, = axes.plot(xs, ys, color='black', linewidth=1, transform=transform)
outline_artists.append(line)
# ===== 更新標題 =====
n_cp = len(snapshots)
axes.set_title(
f"{beam_labels[current_beam]} | "
f"CP {cp_index + 1}/{n_cp} | "
f"Gantry {snap['gantry_angle']:.1f}° "
f"Collimator {snap['collimator_angle']:.1f}° | "
f"累積 MU {snap['cumulative_mu_percent']:.1f}%"
)
figure.canvas.draw_idle()
In [ ]:
def on_beam_selected(label):
"""RadioButtons 被點擊時，重新載入該 Beam 的資料。"""
global snapshots, beam_meterset, current_beam, cp_slider
# 清除舊輪廓
for artist in outline_artists:
artist.remove()
outline_artists.clear()
# 斷開舊滑桿事件（不這麼做的話換 beam 後滑桿會卡住）
if cp_slider is not None:
try:
cp_slider.disconnect_events()
except Exception:
pass
# 載入新 Beam 的資料
current_beam = int(label.split(':')[0])
snapshots, beam_meterset = build_fluence_for_beam(current_beam)
# 用新的 CP 數量重建滑桿
slider_axes.clear()
cp_slider = Slider(
ax=slider_axes, label='Control Point',
valmin=0, valmax=max(len(snapshots) - 1, 1),
valinit=0, valstep=1,
)
cp_slider.on_changed(draw_control_point)
# 第一次繪製
draw_control_point(0)
figure.canvas.draw_idle()
6-4. 啟動介面¶
In [ ]:
def launch_viewer():
global figure, axes, slider_axes, fluence_image, beam_labels
if plan.Modality != 'RTPLAN':
print("錯誤：這個 DICOM 檔不是 RTPLAN。")
return
beam_labels = [f"{i}: {b.BeamName}" for i, b in enumerate(plan.BeamSequence)]
# 建立主圖
figure, axes = plt.subplots(figsize=(6, 5))
plt.subplots_adjust(left=0.25, bottom=0.25)
axes.set_aspect('equal', adjustable='box')
axes.set_xlabel('X (mm)')
axes.set_ylabel('Y (mm)')
axes.set_xlim(BEV_MIN_MM, BEV_MAX_MM)
axes.set_ylim(BEV_MIN_MM, BEV_MAX_MM)
# 空白 fluence 圖（之後用 set_data 原地更新）
fluence_image = axes.imshow(
np.zeros((IMAGE_SIZE, IMAGE_SIZE)),
extent=EXTENT, origin='lower',
cmap='jet', alpha=0.7, vmax=1.0,
)
cbar = figure.colorbar(fluence_image, ax=axes, fraction=0.046, pad=0.04)
cbar.set_label('累積 Fluence（MU）')
# Beam 選擇按鈕（左側）
radio_axes = figure.add_axes([0.02, 0.25, 0.2, 0.2])
radio = RadioButtons(radio_axes, beam_labels, active=0)
radio.on_clicked(on_beam_selected)
# CP 滑桿（底部，先建立佔位）
slider_axes = figure.add_axes([0.25, 0.10, 0.70, 0.04])
on_beam_selected(beam_labels[0])
plt.show()
Step 7｜啟動！¶
In [ ]:
launch_viewer()
常見問題排查¶
| 症狀 | 原因 | 解法 |
|---|---|---|
| 打開後沒有任何視窗 | %matplotlib widget 沒有生效 |
確認有安裝 ipympl，重新執行 Step 1 |
| 滑桿拖不動 | 後端不是 widget | 執行 %matplotlib widget 後 Restart Kernel 再試 |
| 切換 Beam 後畫面卡住 | 舊滑桿事件尚未斷開 | 本版已自動處理；若仍異常請 Restart Kernel |
出現 AttributeError: 'Dataset' ... |
該機台 DICOM 缺少某些標籤 | 程式已用 hasattr / getattr 保護，若仍出現請換一個檔案測試 |
| 整張圖都是空白 | BeamMeterset 為 0，或沒有 FractionGroupSequence |
在 Step 3 先印出每支 Beam 的 MU 值確認 |
| 圖形位置與預期差 1 pixel | DICOM 中 BEV 原點與假設不同 | 調整 BEV_MIN_MM / BEV_MAX_MM 兩個常數 |
換機台時要改什麼？¶
只要改 Step 4-1 的兩個常數即可（其他程式碼完全不用動）：
BEV_MIN_MM = -100.0 # 例如 Neusoft 照野較小
BEV_MAX_MM = 100.0
如果解析度想更細，把 IMAGE_SIZE 加大（例如 800）即可。
延伸練習¶
- 換 colormap：把 Step 6-4 的
cmap='jet'改成'viridis'、'inferno'或'gray' - 存成圖片：
plt.imsave('cp20.png', snapshots[20]['fluence_cumulative'], cmap='jet', origin='lower') - 顯示單段 fluence：在
draw_control_point中把fluence_cumulative改成fluence_cumulative - snapshots[cp_index-1]['fluence_cumulative']（cp_index > 0時） - 加入 Gantry 旋轉：把
rotate_deg_around(0, 0, -snap['collimator_angle'])改成rotate_deg_around(0, 0, -snap['gantry_angle'] - snap['collimator_angle'])，可以看到從患者角度看的照野