DMLC Fluence 演算法教學 Notebook¶
本 notebook 示範兩種動態多葉準直儀（DMLC）的 fluence 生成演算法：
- Sliding Window（Convery-Rosenbloom / Spirou-Chui）— 葉片持續移動，
BeamType = DYNAMIC - Step-and-Shoot（Bortfeld / Engel 分解）— 葉片在 segment 間靜止，
BeamType = STATIC
輸入：target_fluence.png、demo.dcm
輸出：output_sliding.dcm、output_stepshoot.dcm
0. 環境準備¶
In [ ]:
# !pip install numpy pillow matplotlib scipy pydicom
import copy, datetime, logging
from dataclasses import dataclass, field
from typing import List, Optional, Tuple
import numpy as np
import pydicom
from PIL import Image
import matplotlib.pyplot as plt
from pydicom.dataset import Dataset
from pydicom.sequence import Sequence
from pydicom.uid import ExplicitVRLittleEndian, generate_uid
logging.basicConfig(level=logging.INFO, format='%(message)s')
np.set_printoptions(precision=2, suppress=True, linewidth=140)
plt.rcParams['figure.dpi'] = 100
print('pydicom', pydicom.__version__)
0.1 參數設定¶
In [ ]:
@dataclass
class Params:
rp_path: str = 'demo.dcm'
image_path: str = 'target_fluence.png'
out_sliding: str = 'output_sliding.dcm'
out_stepshoot: str = 'output_stepshoot.dcm'
field_x_mm: float = 100.0
field_y_mm: float = 100.0
field_center_x_mm: float = 0.0
field_center_y_mm: float = 0.0
pixel_mm: float = 2.0
invert: bool = True
gamma: float = 1.0
total_mu: float = 100.0
dose_rate_mu_min: Optional[float] = None
max_leaf_speed_mm_s: float = 25.0
leaf_speed_margin: float = 0.9
min_leaf_gap_mm: float = 1.0
cp_candidates: Tuple[int, ...] = (50, 100, 200, 400)
target_nrmse_pct: float = 2.0
max_segments: int = 60
min_segment_mu: float = 1.0
max_levels: int = 32
refine_weights: bool = True
beam_index: Optional[int] = None
new_study_uid: bool = False
p = Params()
print(p)
1. 讀取輸入：RT Plan (MLC 幾何)¶
In [ ]:
@dataclass
class MlcInfo:
device_type: str
boundaries: np.ndarray
jaw_types: List[str] = field(default_factory=list)
@property
def n_pairs(self): return len(self.boundaries) - 1
MLC_TYPES = ('MLCX', 'MLCY')
def _bld_type(bld): return str(bld.RTBeamLimitingDeviceType).upper()
def find_beam(rp, beam_index=None):
seq = rp.BeamSequence
if beam_index is not None: return beam_index, seq[beam_index]
for i, beam in enumerate(seq):
if any(_bld_type(b) in MLC_TYPES for b in beam.BeamLimitingDeviceSequence):
return i, beam
raise ValueError('RP 裡沒有任何帶 MLC 的 beam')
def read_mlc_info(beam) -> MlcInfo:
mlc_type, boundaries, jaws = None, None, []
for bld in beam.BeamLimitingDeviceSequence:
t = _bld_type(bld)
if t in MLC_TYPES:
b = getattr(bld, 'LeafPositionBoundaries', None)
if b is None: raise ValueError(f'{t} 缺少 LeafPositionBoundaries')
mlc_type = t
boundaries = np.array([float(v) for v in b])
else:
jaws.append(t)
if mlc_type is None: raise ValueError('此 beam 沒有 MLC 裝置')
return MlcInfo(mlc_type, boundaries, jaws)
def original_dose_rate(beam):
try:
v = getattr(beam.ControlPointSequence[0], 'DoseRateSet', None)
return float(v) if v is not None else None
except Exception:
return None
rp = pydicom.dcmread(p.rp_path)
beam_index, beam = find_beam(rp, p.beam_index)
mlc = read_mlc_info(beam)
dose_rate = p.dose_rate_mu_min or original_dose_rate(beam) or 600.0
print(f"beam #{beam.BeamNumber} '{getattr(beam,'BeamName','')}'")
print(f'MLC = {mlc.device_type}，{mlc.n_pairs} 對葉片，jaw = {mlc.jaw_types}')
print(f'dose rate = {dose_rate:.0f} MU/min')
2. 圖像 → 目標 Fluence Map¶
In [ ]:
@dataclass
class FluenceMap:
data: np.ndarray
x_edges: np.ndarray
y_edges: np.ndarray
@property
def x_centers(self): return 0.5*(self.x_edges[:-1] + self.x_edges[1:])
@property
def y_centers(self): return 0.5*(self.y_edges[:-1] + self.y_edges[1:])
def load_fluence(p: Params) -> FluenceMap:
nx = int(round(p.field_x_mm / p.pixel_mm))
ny = int(round(p.field_y_mm / p.pixel_mm))
img = Image.open(p.image_path)
if img.mode in ('RGBA', 'LA') or 'transparency' in img.info:
rgba = img.convert('RGBA')
bg = Image.new('RGBA', rgba.size, (255, 255, 255, 255))
img = Image.alpha_composite(bg, rgba)
img = img.convert('L').resize((nx, ny), Image.LANCZOS)
arr = np.asarray(img, dtype=np.float64) / 255.0
f = (1.0 - arr) if p.invert else arr
f = np.clip(f, 0.0, 1.0) ** p.gamma
f = f[::-1, :]
if f.max() <= 1e-9: raise ValueError('fluence 全為 0')
x0 = p.field_center_x_mm - p.field_x_mm/2
y0 = p.field_center_y_mm - p.field_y_mm/2
return FluenceMap(
data=np.ascontiguousarray(f),
x_edges=x0 + np.arange(nx+1)*(p.field_x_mm/nx),
y_edges=y0 + np.arange(ny+1)*(p.field_y_mm/ny))
fm = load_fluence(p)
print(f'fluence 形狀 = {fm.data.shape}，值域 = [{fm.data.min():.2f}, {fm.data.max():.2f}]')
plt.figure(figsize=(5, 4))
plt.imshow(fm.data, origin='lower', cmap='gray_r',
extent=[fm.x_edges[0], fm.x_edges[-1], fm.y_edges[0], fm.y_edges[-1]])
plt.colorbar(label='fluence (0-1)')
plt.xlabel('X (mm)'); plt.ylabel('Y (mm)')
plt.title('目標 Fluence Map')
plt.show()
2.1 抽取每片葉片的 1D 目標 fluence¶
In [ ]:
def leaf_row_profiles(data, row_edges, boundaries):
n_rows, n_cols = data.shape
dy = row_edges[1] - row_edges[0]
csum = np.vstack([np.zeros((1, n_cols)), np.cumsum(data, axis=0)])
lo = np.clip(boundaries[:-1], row_edges[0], row_edges[-1])
hi = np.clip(boundaries[1:], row_edges[0], row_edges[-1])
def cum_at(y):
pos = (y - row_edges[0]) / dy
i = np.clip(np.floor(pos).astype(int), 0, n_rows - 1)
frac = (pos - i)[:, None]
return csum[i] + frac*(csum[i+1] - csum[i])
area = cum_at(hi) - cum_at(lo)
width = (hi - lo) / dy
covered = width > 1e-6
profiles = np.zeros_like(area)
profiles[covered] = area[covered] / width[covered, None]
return profiles, covered
if mlc.device_type == 'MLCX':
data, row_edges = fm.data, fm.y_edges
row_centers, travel = fm.y_centers, fm.x_centers
else:
data, row_edges = fm.data.T, fm.x_edges
row_centers, travel = fm.x_centers, fm.y_centers
profiles, covered = leaf_row_profiles(data, row_edges, mlc.boundaries)
active = covered & (profiles.max(axis=1) > 1e-6)
print(f'需出束葉片對 = {active.sum()} / {len(active)}')
# 全景：所有 active 葉片的 1D profile
fig, axes = plt.subplots(1, 2, figsize=(14, 4), gridspec_kw={'width_ratios': [3, 1]})
active_rows = np.where(active)[0]
for r in active_rows:
axes[0].plot(travel, profiles[r], alpha=0.4)
axes[0].set_xlabel('travel (mm)'); axes[0].set_ylabel('fluence (0-1)')
axes[0].set_title(f'所有 {active.sum()} 片 active 葉片的 1D 目標 profile')
axes[0].grid(alpha=0.3)
# 側視圖：每片葉片的平均強度
row_mean = profiles.mean(axis=1)
axes[1].barh(np.arange(mlc.n_pairs), row_mean,
color=['C0' if a else 'lightgray' for a in active])
axes[1].set_xlabel('平均 fluence')
axes[1].set_ylabel('leaf pair #')
axes[1].set_title('每片葉片平均強度')
axes[1].invert_yaxis()
axes[1].grid(alpha=0.3, axis='x')
plt.tight_layout(); plt.show()
In [ ]:
def mu_per_mm(dose_rate_mu_min, max_speed_mm_s, margin=0.9):
return (dose_rate_mu_min / 60.0) / (max_speed_mm_s * margin)
g = mu_per_mm(dose_rate, p.max_leaf_speed_mm_s, p.leaf_speed_margin)
print(f'g = {g:.4f} MU/mm')
print(f'→ 掃過 {p.field_x_mm:.0f} mm 照野要 {g*p.field_x_mm:.1f} MU')
3.1 計算累積 MU 軌跡 T_A(x)、T_B(x)¶
In [ ]:
@dataclass
class Trajectories:
x: np.ndarray
t_a: np.ndarray
t_b: np.ndarray
active: np.ndarray
scale: float
total_mu: float
g: float
def compute_trajectories(profiles, active, x_centers, total_mu, g, min_gap_mm):
dx = float(x_centers[1] - x_centers[0])
pad_n = int(np.ceil(min_gap_mm / dx)) + 1
x = np.concatenate([
x_centers[0] - dx*np.arange(pad_n, 0, -1),
x_centers,
x_centers[-1] + dx*np.arange(1, pad_n+1)])
P = np.pad(profiles, ((0,0), (pad_n, pad_n)))
P[~active] = 0.0
d = np.diff(P, axis=1)
rise = np.maximum(d, 0.0).sum(axis=1)
travel_mu = g * (x[-1] - x[0])
if total_mu <= travel_mu * 1.001:
raise ValueError(f'total_mu={total_mu:.1f} 太小：光是掃過就要 {travel_mu:.1f} MU')
scale = (total_mu - travel_mu) / rise[active].max()
step = g * dx
zeros = np.zeros((P.shape[0], 1))
t_a = np.concatenate([zeros, np.cumsum(step + scale*np.maximum( d, 0.0), axis=1)], axis=1)
t_b = np.concatenate([zeros, np.cumsum(step + scale*np.maximum(-d, 0.0), axis=1)], axis=1)
return Trajectories(x=x, t_a=t_a, t_b=t_b, active=active.copy(),
scale=float(scale), total_mu=float(total_mu), g=float(g))
traj = compute_trajectories(profiles, active, travel, p.total_mu, g, p.min_leaf_gap_mm)
print(f'scale = {traj.scale:.3f}（fluence 1.0 → {traj.scale:.1f} MU）')
print(f'beam 總 MU = {traj.total_mu:.1f}')
3.2 視覺化軌跡（教學重點）¶
觀察 T_A − T_B 如何在每個位置 x 上剛好等於目標 fluence，以及所有葉片的軌跡疊起來的樣子。
In [ ]:
r = int(np.argmax(active))
target_rows = profiles * traj.scale
target_rows[~active] = 0.0
fig, axes = plt.subplots(1, 3, figsize=(18, 4))
# (a) 單一葉片：目標 fluence
axes[0].fill_between(travel, 0, target_rows[r], alpha=0.4, color='C0')
axes[0].plot(travel, target_rows[r], 'C0-', label=f'目標 × {traj.scale:.1f} MU')
axes[0].set_xlabel('travel (mm)'); axes[0].set_ylabel('MU')
axes[0].set_title(f'Row {r}：目標 fluence')
axes[0].legend(); axes[0].grid(alpha=0.3)
# (b) 單一葉片：T_A / T_B 軌跡
axes[1].plot(traj.x, traj.t_a[r], 'C3-', label='T_A(x) 左葉')
axes[1].plot(traj.x, traj.t_b[r], 'C2-', label='T_B(x) 右葉')
axes[1].fill_between(traj.x, traj.t_b[r], traj.t_a[r], color='C0', alpha=0.3,
label='T_A − T_B = fluence')
axes[1].set_xlabel('travel (mm)'); axes[1].set_ylabel('累積 MU')
axes[1].set_title('兩葉的累積 MU 軌跡')
axes[1].legend(); axes[1].grid(alpha=0.3)
# (c) 所有 active 葉片的 T_A / T_B 疊在一起
for rr in active_rows:
axes[2].plot(traj.x, traj.t_a[rr], 'C3-', alpha=0.25, lw=0.8)
axes[2].plot(traj.x, traj.t_b[rr], 'C2-', alpha=0.25, lw=0.8)
axes[2].plot([], [], 'C3-', label='T_A 全部')
axes[2].plot([], [], 'C2-', label='T_B 全部')
axes[2].set_xlabel('travel (mm)'); axes[2].set_ylabel('累積 MU')
axes[2].set_title('所有 active 葉片疊起來')
axes[2].legend(); axes[2].grid(alpha=0.3)
plt.tight_layout(); plt.show()
3.3 在每個 CP 取樣葉片位置¶
In [ ]:
def sample_leaf_positions(traj, t_mu, min_gap_mm, decimals=2):
n = traj.t_a.shape[0]
A = np.empty((len(t_mu), n)); B = np.empty((len(t_mu), n))
for r in range(n):
if traj.active[r]:
A[:, r] = np.interp(t_mu, traj.t_a[r], traj.x)
B[:, r] = np.interp(t_mu, traj.t_b[r], traj.x)
else:
A[:, r] = B[:, r] = traj.x[0]
q = 10 ** decimals
gap = np.ceil(min_gap_mm * q) / q
A, B = np.round(A, decimals), np.round(B, decimals)
A = np.round(np.minimum(A, B - gap), decimals)
return A, B
n_cp = 100
t_mu = np.linspace(0.0, traj.total_mu, n_cp)
A_slide, B_slide = sample_leaf_positions(traj, t_mu, p.min_leaf_gap_mm)
print(f'A shape = {A_slide.shape}')
print(f'CP0 A[:3]={A_slide[0,:3]}, B[:3]={B_slide[0,:3]}')
print(f'CP99 A[:3]={A_slide[-1,:3]}, B[:3]={B_slide[-1,:3]}')
3.3b Aperture 快照：5 個關鍵 CP 的葉片位置¶
類似 MLC 軟體畫面：每一列表示一片葉片，藍線是打開的區間 [A, B]。可以看到葉片如何從左往右掃過照野。
In [ ]:
fig, axes = plt.subplots(1, 5, figsize=(20, 5), sharey=True)
cp_shots = [0, n_cp//4, n_cp//2, 3*n_cp//4, n_cp-1]
for ax, k in zip(axes, cp_shots):
for rr in range(mlc.n_pairs):
if active[rr]:
ax.plot([A_slide[k, rr], B_slide[k, rr]], [rr, rr], 'C0-', lw=2.5)
else:
ax.plot([A_slide[k, rr], B_slide[k, rr]], [rr, rr], color='lightgray', lw=0.5)
ax.set_title(f'CP {k}\nCumW={k/(n_cp-1):.2f}')
ax.set_xlabel('travel (mm)')
ax.set_xlim(travel[0]-5, travel[-1]+5)
ax.grid(alpha=0.3)
axes[0].set_ylabel('leaf pair #')
plt.suptitle('Sliding window：葉片掃過的過程（藍線 = 打開區間）', y=1.02)
plt.tight_layout(); plt.show()
3.4 正向模擬與誤差¶
In [ ]:
def simulate_rows(A, B, t_mu, x_centers):
n = A.shape[1]
out = np.zeros((n, len(x_centers)))
for r in range(n):
if np.ptp(A[:, r]) == 0 and np.ptp(B[:, r]) == 0: continue
ta = np.interp(x_centers, A[:, r], t_mu)
tb = np.interp(x_centers, B[:, r], t_mu)
out[r] = np.maximum(ta - tb, 0.0)
return out
def compare(target_mu, sim_mu, active):
diff = (sim_mu - target_mu)[active]
peak = target_mu[active].max()
return {'nrmse_pct': float(np.sqrt(np.mean(diff**2))/peak*100),
'max_abs_err_pct': float(np.abs(diff).max()/peak*100)}
sim_slide = simulate_rows(A_slide, B_slide, t_mu, travel)
m_slide = compare(target_rows, sim_slide, active)
print(f"CP={n_cp}: NRMSE = {m_slide['nrmse_pct']:.2f}%，最大誤差 = {m_slide['max_abs_err_pct']:.2f}%")
3.4b 目標 vs 模擬：單葉片 + 全體¶
左：單一葉片的 1D 曲線對比；中：全體 row × travel 的模擬結果；右：殘差。
In [ ]:
fig, axes = plt.subplots(1, 3, figsize=(18, 4.5))
# (a) 單葉片對比
axes[0].fill_between(travel, 0, target_rows[r], alpha=0.3, color='C0', label='目標')
axes[0].plot(travel, sim_slide[r], 'C3-', lw=1.5, label='模擬')
axes[0].set_xlabel('travel (mm)'); axes[0].set_ylabel('MU')
axes[0].set_title(f'Row {r}：目標 vs 模擬')
axes[0].legend(); axes[0].grid(alpha=0.3)
# (b) 全體 heatmap
ext_t = [travel[0], travel[-1], -0.5, mlc.n_pairs-0.5]
vmax = float(target_rows.max())
im1 = axes[1].imshow(sim_slide, origin='lower', aspect='auto', cmap='gray_r',
vmin=0, vmax=vmax, extent=ext_t)
axes[1].set_xlabel('travel (mm)'); axes[1].set_ylabel('leaf pair #')
axes[1].set_title(f'Sliding 模擬結果（{n_cp} CP）')
fig.colorbar(im1, ax=axes[1], fraction=0.046, label='MU')
# (c) 殘差
res = sim_slide - target_rows
im2 = axes[2].imshow(res, origin='lower', aspect='auto', cmap='coolwarm',
vmin=-0.1*vmax, vmax=0.1*vmax, extent=ext_t)
axes[2].set_xlabel('travel (mm)'); axes[2].set_ylabel('leaf pair #')
axes[2].set_title('殘差（模擬 − 目標）')
fig.colorbar(im2, ax=axes[2], fraction=0.046, label='ΔMU')
plt.tight_layout(); plt.show()
3.5 CP 加密實驗：CP 越多誤差越小？¶
In [ ]:
print('CP 數 vs NRMSE：')
cand_results = []
best = None
cp_imgs = {}
for n in p.cp_candidates:
t = np.linspace(0.0, traj.total_mu, n)
A_, B_ = sample_leaf_positions(traj, t, p.min_leaf_gap_mm)
s = simulate_rows(A_, B_, t, travel)
mm = compare(target_rows, s, active)
cand_results.append((n, mm['nrmse_pct']))
cp_imgs[n] = s
print(f" CP={n:4d} NRMSE={mm['nrmse_pct']:5.2f}% 最大誤差={mm['max_abs_err_pct']:5.2f}%")
if best is None or mm['nrmse_pct'] < best[0]['nrmse_pct']:
best = (mm, t, A_, B_, s)
if mm['nrmse_pct'] <= p.target_nrmse_pct: break
m_slide, t_mu, A_slide, B_slide, sim_slide = best
print(f"\n最終採用：{len(t_mu)} CP，NRMSE={m_slide['nrmse_pct']:.2f}%")
fig, axes = plt.subplots(1, 2, figsize=(14, 4))
axes[0].plot([r[0] for r in cand_results], [r[1] for r in cand_results], 'o-')
axes[0].set_xscale('log')
axes[0].set_xlabel('control point 數'); axes[0].set_ylabel('NRMSE (%)')
axes[0].axhline(p.target_nrmse_pct, color='r', ls='--', label=f'目標 {p.target_nrmse_pct}%')
axes[0].legend(); axes[0].grid(alpha=0.3, which='both')
axes[0].set_title('CP 數對誤差的影響')
# 不同 CP 數的模擬結果對比
n_show = len(cp_imgs)
for i, (n, s) in enumerate(cp_imgs.items()):
mm = compare(target_rows, s, active)
axes[1].plot(travel, s[r], label=f'{n} CP (NRMSE {mm["nrmse_pct"]:.1f}%)', alpha=0.7)
axes[1].plot(travel, target_rows[r], 'k--', lw=2, label='目標')
axes[1].set_xlabel('travel (mm)'); axes[1].set_ylabel('MU')
axes[1].set_title(f'Row {r}：不同 CP 數的模擬結果')
axes[1].legend(fontsize=8); axes[1].grid(alpha=0.3)
plt.tight_layout(); plt.show()
In [ ]:
weights_slide = np.round(t_mu / traj.total_mu, 6)
weights_slide[0], weights_slide[-1] = 0.0, 1.0
print(f'weights[0]={weights_slide[0]}, weights[-1]={weights_slide[-1]}')
print(f'n_cp = {len(weights_slide)}')
4.1 一維分解示範（教學重點）¶
把 [0, 1, 2, 3, 2, 1, 0] 這種 1D profile 拆成 6 個 segment，每個 segment 對應「某一層」，全部疊起來剛好還原原 profile。
In [ ]:
def decompose_row(levels):
nx = len(levels)
d = np.diff(np.concatenate([[0], levels, [0]]))
edges = np.arange(nx + 1)
return (np.repeat(edges, np.maximum( d, 0)),
np.repeat(edges, np.maximum(-d, 0)))
demo = np.array([0, 1, 2, 3, 2, 1, 0])
L, R = decompose_row(demo)
print('levels =', demo)
print('左邊界 =', L)
print('右邊界 =', R)
print(f'→ {len(L)} 個 segment')
fig, axes = plt.subplots(3, 1, figsize=(11, 8.5))
xp = np.arange(len(demo))
# (a) 量化後的 profile
axes[0].bar(xp, demo, color='lightsteelblue', edgecolor='navy', width=0.8)
for i, v in enumerate(demo):
if v > 0: axes[0].text(i, v+0.08, str(v), ha='center', fontsize=9)
axes[0].set_ylabel('level')
axes[0].set_title('(a) 量化後的 1D profile I[j]')
axes[0].set_ylim(0, demo.max()+0.5)
axes[0].grid(alpha=0.3, axis='y')
# (b) 每個 segment 的區間
for k, (l, rr) in enumerate(zip(L, R)):
axes[1].barh(k, rr-l, left=l-0.5, height=0.8,
color=plt.cm.Set2(k), edgecolor='k')
axes[1].text(l-0.5+(rr-l)/2, k, f'[{l},{rr})', ha='center', va='center', fontsize=9)
axes[1].set_xlim(-0.5, len(demo)-0.5)
axes[1].set_ylim(-0.6, len(L)-0.4)
axes[1].set_ylabel('segment #')
axes[1].set_title(f'(b) 每個 segment 的開口 [L_k, R_k)（共 {len(L)} 個）')
axes[1].grid(alpha=0.3, axis='x')
# (c) Σ segments = levels
acc = np.zeros_like(demo)
for k, (l, rr) in enumerate(zip(L, R)):
seg = np.zeros_like(demo, dtype=int)
seg[l:rr] = 1
acc += seg
axes[2].bar(xp, demo, color='lightsteelblue', edgecolor='navy',
width=0.8, alpha=0.6, label='目標 levels')
axes[2].step(xp, acc, where='mid', color='C3', lw=2.5, label='Σ segments')
axes[2].set_ylabel('level'); axes[2].set_xlabel('pixel index')
axes[2].set_title('(c) 驗證：Σ 1[L_k ≤ x < R_k] = I(x) 精確成立')
axes[2].legend(); axes[2].grid(alpha=0.3, axis='y')
plt.tight_layout(); plt.show()
4.1b 量化誤差示範：層數越多越接近原 profile¶
拿一片真實葉片，展示不同層數的量化誤差。這也說明為什麼後面要「挑層數」。
In [ ]:
r_demo = int(np.argmax(active))
prof_r = profiles[r_demo]
fig, axes = plt.subplots(1, 4, figsize=(16, 4))
for ax, n in zip(axes, [4, 8, 16, 32]):
lv = np.rint(prof_r * n).astype(int)
ax.fill_between(travel, 0, prof_r * n, alpha=0.3, color='C0', label=f'連續 × {n}')
ax.step(travel, lv, where='mid', color='C3', lw=1.6, label='量化後')
err = prof_r * n - lv
nrmse = np.sqrt(np.mean(err**2)) / (prof_r * n).max() * 100
ax.set_title(f'n_levels = {n}\n量化 NRMSE = {nrmse:.1f}%')
ax.set_xlabel('travel (mm)'); ax.set_ylabel('level')
ax.grid(alpha=0.3); ax.legend(fontsize=8)
plt.suptitle(f'Row {r_demo}：不同量化層數的擬合效果')
plt.tight_layout(); plt.show()
4.2 對所有葉片做分解¶
In [ ]:
def quantize_rows(profiles, active, n_levels):
levels = np.rint(profiles * n_levels).astype(int)
levels[~active] = 0
return levels
def build_segments(levels, x_edges, min_gap_mm):
n_leaves = levels.shape[0]
rows = [decompose_row(levels[r]) for r in range(n_leaves)]
M = max(len(L) for L, _ in rows)
A = np.full((M, n_leaves), x_edges[0] - min_gap_mm)
B = np.full((M, n_leaves), x_edges[0])
is_open = np.zeros((M, n_leaves), dtype=bool)
for r, (L, R) in enumerate(rows):
k = len(L)
if k == 0: continue
a, b = x_edges[L], x_edges[R]
b = np.maximum(b, a + min_gap_mm)
A[:k, r], B[:k, r], is_open[:k, r] = a, b, True
return np.round(A, 2), np.round(B, 2), is_open
def segment_footprints(A, B, x_edges):
e0, e1 = x_edges[:-1], x_edges[1:]
dx = e1[0] - e0[0]
ov = np.minimum(B[:, :, None], e1) - np.maximum(A[:, :, None], e0)
return np.clip(ov, 0.0, dx) / dx
dx = float(travel[1] - travel[0])
x_edges = travel[0] - dx/2 + np.arange(len(travel)+1) * dx
n_test = [2, 4, 8, 12, 16, 24, 32, 48, 64]
seg_counts = []
for n in n_test:
lv = quantize_rows(profiles, active, n)
A_, B_, _ = build_segments(lv, x_edges, p.min_leaf_gap_mm)
seg_counts.append(A_.shape[0])
print(f' n_levels={n:3d} → M={A_.shape[0]:4d} segments')
fig, ax = plt.subplots(figsize=(9, 4))
ax.bar(range(len(n_test)), seg_counts, color='steelblue', edgecolor='navy')
ax.axhline(p.max_segments, color='r', ls='--', label=f'max_segments = {p.max_segments}')
ax.set_xticks(range(len(n_test))); ax.set_xticklabels(n_test)
ax.set_xlabel('n_levels（量化層數）'); ax.set_ylabel('segment 數 M')
ax.set_title('量化層數 vs segment 數')
ax.legend(); ax.grid(alpha=0.3, axis='y')
for i, c in enumerate(seg_counts):
ax.text(i, c+0.5, str(c), ha='center', fontsize=9)
plt.tight_layout(); plt.show()
4.2b 各葉片需要多少 segment？¶
不同葉片所需的 segment 數不同；整支 beam 的 M = 各列的最大值。
In [ ]:
# 用最終選定的層數算（見 4.3；這裡先用 16 示範）
n_demo = 16
lv_demo = quantize_rows(profiles, active, n_demo)
per_row = []
for r_ in range(mlc.n_pairs):
L_, R_ = decompose_row(lv_demo[r_])
per_row.append(len(L_))
per_row = np.array(per_row)
fig, ax = plt.subplots(figsize=(12, 4))
colors = ['C0' if a else 'lightgray' for a in active]
ax.bar(np.arange(mlc.n_pairs), per_row, color=colors, edgecolor='k')
ax.axhline(per_row.max(), color='r', ls='--',
label=f'最大值 M = {per_row.max()}')
ax.set_xlabel('leaf pair #'); ax.set_ylabel('該列需要的 segment 數')
ax.set_title(f'各葉片在 n_levels = {n_demo} 時所需的 segment 數')
ax.legend(); ax.grid(alpha=0.3, axis='y')
plt.tight_layout(); plt.show()
print(f'最多需要 {per_row.max()} 個 segment 的葉片 # = {int(np.argmax(per_row))}')
4.3 自動挑選層數¶
選最大 n_levels 使得 segment 數 ≤ max_segments 且每 segment MU ≥ min_segment_mu。
In [ ]:
chosen = None
for n in range(p.max_levels, 1, -1):
lv = quantize_rows(profiles, active, n)
A_, B_, is_open_ = build_segments(lv, x_edges, p.min_leaf_gap_mm)
M_ = A_.shape[0]
if M_ <= p.max_segments and p.total_mu / M_ >= p.min_segment_mu:
chosen = (n, lv, A_, B_, is_open_, M_)
break
if chosen is None:
raise ValueError('找不到符合條件的層數')
n_levels, levels, A_sas, B_sas, is_open_sas, M = chosen
print(f'選定：{n_levels} 層，{M} 個 segment，等權重約 {p.total_mu/M:.2f} MU/seg')
4.4 視覺化 segment¶
- 前 8 個 segment 的 aperture（哪些葉片打開）
- 每個 segment 打開幾對葉片
- 逐步累積：從 1 個 segment 累積到全部，觀察逼近目標的過程
In [ ]:
# 準備常數（供後續視覺化使用）
r_demo = int(np.argmax(active))
target0_local = profiles * (n_levels * p.total_mu / M)
target0_local[~active] = 0.0
w_eq = np.full(M, p.total_mu / M)
def plot_segments(A, B, is_open, n_show=8):
M_ = A.shape[0]; n = min(M_, n_show)
fig, axes = plt.subplots(n, 1, figsize=(10, 1.1*n), sharex=True)
if n == 1: axes = [axes]
for k in range(n):
ax = axes[k]
for rr in range(A.shape[1]):
if is_open[k, rr]:
ax.plot([A[k, rr], B[k, rr]], [rr, rr], 'C0-', lw=3)
ax.set_ylim(-1, A.shape[1])
ax.set_ylabel(f'seg {k}')
ax.grid(alpha=0.2)
axes[-1].set_xlabel('travel (mm)')
plt.tight_layout(); plt.show()
plot_segments(A_sas, B_sas, is_open_sas, n_show=8)
print(f'（只顯示前 8 個，共 {M} 個 segment）')
# 每 segment 打開幾個葉片對
open_count = is_open_sas.sum(axis=1)
fig, ax = plt.subplots(figsize=(14, 3.5))
ax.bar(np.arange(M), open_count, color='teal', edgecolor='k', alpha=0.8)
ax.set_xlabel('segment #'); ax.set_ylabel('打開的葉片對數')
ax.set_title(f'每個 segment 打開幾對葉片（共 {M} 個 segment）')
ax.grid(alpha=0.3, axis='y')
plt.tight_layout(); plt.show()
4.4b 逐步累積：從 1 個 segment 到全部¶
觀察 sum(前 k 個 segment 貢獻) 如何一步步逼近目標 profile。
In [ ]:
F_sas_full = segment_footprints(A_sas, B_sas, x_edges) # (M, n_leaves, nx)
F_demo = F_sas_full[:, r_demo, :] # (M, nx)
cums = sorted(set([1, max(1, M//8), max(1, M//4),
max(1, M//2), max(1, 3*M//4), M]))
cums = [c for c in cums if c > 0]
fig, axes = plt.subplots(2, 3, figsize=(15, 6.5))
for ax, c in zip(axes.ravel(), cums):
sim_c = np.tensordot(w_eq[:c], F_demo[:c], axes=1)
ax.fill_between(travel, 0, target0_local[r_demo], alpha=0.3, color='C0', label='目標')
ax.plot(travel, sim_c, 'C3-', lw=1.6, label=f'累積前 {c} seg')
ax.set_title(f'前 {c} / {M} 個 segment')
ax.set_xlabel('travel (mm)'); ax.set_ylabel('MU')
ax.legend(fontsize=8); ax.grid(alpha=0.3)
plt.suptitle(f'Row {r_demo}：逐步累積 segment 貢獻（越後面越接近目標）')
plt.tight_layout(); plt.show()
4.4c 全景圖：所有 segment 的葉片開放狀態¶
深色 = 該對葉片在這個 segment 打開。可以看到 MLC 在不同 segment 打開不同區域。
In [ ]:
n_show = min(M, 40)
fig, ax = plt.subplots(figsize=(14, max(3, 0.22 * mlc.n_pairs)))
im = ax.imshow(is_open_sas[:n_show].T, aspect='auto', cmap='Blues',
interpolation='nearest')
ax.set_xlabel('segment #（時間順序）'); ax.set_ylabel('leaf pair #')
ax.set_title(f'各 segment 的葉片開放狀態（前 {n_show}/{M} 個）')
ax.set_xticks(range(0, n_show, max(1, n_show//10)))
plt.colorbar(im, ax=ax, fraction=0.03, label='is_open')
plt.tight_layout(); plt.show()
4.5 權重擬合：等權重 vs NNLS¶
- 等權重：每個 segment
MU = total_mu / M - NNLS：最小化模擬與目標差異，剔除 MU 太小的 segment
In [ ]:
F_sas = segment_footprints(A_sas, B_sas, x_edges)
w_equal = np.full(M, p.total_mu / M)
target0 = profiles * (n_levels * p.total_mu / M)
target0[~active] = 0.0
sim_eq = np.tensordot(w_equal, F_sas, axes=1)
diff_eq = (sim_eq - target0)[active]
nrmse_eq = np.sqrt(np.mean(diff_eq**2)) / target0[active].max() * 100
print(f'等權重 NRMSE = {nrmse_eq:.2f}%')
def _nnls_weights(F, target, active, min_mu):
from scipy.optimize import nnls
keep = np.arange(F.shape[0])
b = target[active].ravel()
w = None; low = None
for _ in range(5):
Amat = F[keep][:, active, :].reshape(len(keep), -1).T
w, _ = nnls(Amat, b, maxiter=50*len(keep))
low = w < min_mu
if not low.any() or low.all(): break
keep = keep[~low]
if w is None or (low is not None and low.all()):
return None, None
return keep, w
nrmse_nn = None
if p.refine_weights:
keep, w_nn = _nnls_weights(F_sas, target0, active, p.min_segment_mu)
if keep is not None:
A_sas, B_sas, is_open_sas = A_sas[keep], B_sas[keep], is_open_sas[keep]
F_sas = F_sas[keep]
w_nn = w_nn * (p.total_mu / w_nn.sum())
sim_nn = np.tensordot(w_nn, F_sas, axes=1)
diff_nn = (sim_nn - target0)[active]
nrmse_nn = np.sqrt(np.mean(diff_nn**2)) / target0[active].max() * 100
print(f'NNLS：保留 {len(keep)}/{M} 個 segment，最小 MU={w_nn.min():.2f}，NRMSE={nrmse_nn:.2f}%')
w_sas, sim_sas, M = w_nn, sim_nn, len(keep)
else:
print('NNLS 全低於 min_mu，使用等權重')
w_sas, sim_sas = w_equal, sim_eq
else:
w_sas, sim_sas = w_equal, sim_eq
fig, axes = plt.subplots(1, 3, figsize=(16, 4))
axes[0].bar(np.arange(len(w_equal)), w_equal, color='steelblue', edgecolor='k', alpha=0.8)
axes[0].axhline(p.min_segment_mu, color='r', ls='--', label=f'min = {p.min_segment_mu}')
axes[0].set_xlabel('segment #'); axes[0].set_ylabel('MU')
axes[0].set_title(f'等權重（{len(w_equal)} 個，NRMSE={nrmse_eq:.1f}%）')
axes[0].legend(fontsize=8); axes[0].grid(alpha=0.3, axis='y')
axes[1].bar(np.arange(len(w_sas)), w_sas, color='coral', edgecolor='k', alpha=0.85)
axes[1].axhline(p.min_segment_mu, color='r', ls='--')
axes[1].set_xlabel('segment #'); axes[1].set_ylabel('MU')
title_nn = (f'NNLS（{len(w_sas)} 個，NRMSE={nrmse_nn:.1f}%）'
if nrmse_nn is not None else '等權重（未執行 NNLS）')
axes[1].set_title(title_nn); axes[1].grid(alpha=0.3, axis='y')
axes[2].fill_between(travel, 0, target0[r_demo], alpha=0.3, color='C0', label='目標')
axes[2].plot(travel, sim_eq[r_demo], 'C1-', lw=1.5, label=f'等權重 ({nrmse_eq:.1f}%)')
if nrmse_nn is not None:
axes[2].plot(travel, sim_sas[r_demo], 'C3-', lw=1.5, label=f'NNLS ({nrmse_nn:.1f}%)')
axes[2].set_xlabel('travel (mm)'); axes[2].set_ylabel('MU')
axes[2].set_title(f'Row {r_demo}：等權重 vs NNLS')
axes[2].legend(fontsize=8); axes[2].grid(alpha=0.3)
plt.tight_layout(); plt.show()
4.5b NNLS 修正了哪裡？¶
In [ ]:
if nrmse_nn is not None:
fig, axes = plt.subplots(1, 2, figsize=(14, 4))
ext_t = [travel[0], travel[-1], -0.5, mlc.n_pairs-0.5]
diff_eq_img = sim_eq - target0
diff_nn_img = sim_sas - target0
vmax_d = max(np.abs(diff_eq_img).max(), np.abs(diff_nn_img).max())
for ax, img, title in [
(axes[0], diff_eq_img, f'等權重 − 目標（NRMSE={nrmse_eq:.1f}%）'),
(axes[1], diff_nn_img, f'NNLS − 目標（NRMSE={nrmse_nn:.1f}%）')]:
im = ax.imshow(img, origin='lower', aspect='auto', cmap='coolwarm',
vmin=-vmax_d, vmax=vmax_d, extent=ext_t)
ax.set_xlabel('travel (mm)'); ax.set_ylabel('leaf pair #')
ax.set_title(title)
fig.colorbar(im, ax=ax, fraction=0.046, label='MU 誤差')
plt.suptitle('NNLS 明顯降低了等權重的誤差')
plt.tight_layout(); plt.show()
else:
print('未執行 NNLS，略過此圖')
In [ ]:
n_seg = len(w_sas)
cum = np.round(np.cumsum(w_sas) / w_sas.sum(), 6)
cum[-1] = 1.0
prev = np.concatenate([[0.0], cum[:-1]])
weights_sas = np.empty(2*n_seg)
weights_sas[0::2], weights_sas[1::2] = prev, cum
A_sas_cp = np.repeat(A_sas, 2, axis=0)
B_sas_cp = np.repeat(B_sas, 2, axis=0)
print(f'stepshoot: {n_seg} segments → {len(weights_sas)} CP')
print(f'weights[0:4] = {weights_sas[:4]}')
5. 兩種演算法並排比較¶
In [ ]:
def rows_to_grid(rows_mu, boundaries, row_centers):
idx = np.clip(np.searchsorted(boundaries, row_centers, side='right') - 1,
0, rows_mu.shape[0] - 1)
return rows_mu[idx]
delivered_slide = rows_to_grid(sim_slide, mlc.boundaries, row_centers)
delivered_sas = rows_to_grid(sim_sas, mlc.boundaries, row_centers)
if mlc.device_type == 'MLCY':
delivered_slide = delivered_slide.T
delivered_sas = delivered_sas.T
target_mu = fm.data * traj.scale
assert delivered_slide.shape == target_mu.shape, \
f'{delivered_slide.shape} vs {target_mu.shape}'
assert delivered_sas.shape == target_mu.shape
vmax = float(target_mu.max())
ext = [fm.x_edges[0], fm.x_edges[-1], fm.y_edges[0], fm.y_edges[-1]]
imgs = [
(target_mu, 'Target (MU)',
dict(cmap='gray_r', vmin=0, vmax=vmax)),
(delivered_slide, f'Sliding\n({len(weights_slide)} CP)',
dict(cmap='gray_r', vmin=0, vmax=vmax)),
(delivered_sas, f'Step-and-shoot\n({n_seg} seg)',
dict(cmap='gray_r', vmin=0, vmax=vmax)),
(delivered_slide - target_mu, 'Sliding − Target',
dict(cmap='coolwarm', vmin=-0.15*vmax, vmax=0.15*vmax)),
]
fig, ax = plt.subplots(1, 4, figsize=(18, 4))
for a, (img, title, kw) in zip(ax, imgs):
im = a.imshow(img, origin='lower', extent=ext, **kw)
a.set_title(title); a.set_xlabel('X (mm)'); a.set_ylabel('Y (mm)')
fig.colorbar(im, ax=a, fraction=0.046)
nr_slide = np.sqrt(np.mean((delivered_slide - target_mu)**2)) / vmax * 100
nr_sas = np.sqrt(np.mean((delivered_sas - target_mu)**2)) / vmax * 100
fig.suptitle(f'Sliding NRMSE={nr_slide:.2f}% | Step-and-shoot NRMSE={nr_sas:.2f}%')
fig.tight_layout(); plt.show()
6. 寫回 DICOM¶
In [ ]:
JAW_X = {'X', 'ASYMX'}
JAW_Y = {'Y', 'ASYMY'}
def _r(v, nd): return float(round(float(v), nd))
def _bld_item(dev_type, values):
d = Dataset()
d.RTBeamLimitingDeviceType = dev_type
d.LeafJawPositions = [_r(v, 2) for v in values]
return d
def _fill_cp0_devices(cp0, mlc, mlc_values, jaw_x, jaw_y):
seq = getattr(cp0, 'BeamLimitingDevicePositionSequence', None)
if seq is None:
seq = Sequence(); cp0.BeamLimitingDevicePositionSequence = seq
seen = set()
for item in seq:
t = str(item.RTBeamLimitingDeviceType).upper(); seen.add(t)
if t == mlc.device_type: item.LeafJawPositions = [_r(v, 2) for v in mlc_values]
elif t in JAW_X: item.LeafJawPositions = [_r(v, 2) for v in jaw_x]
elif t in JAW_Y: item.LeafJawPositions = [_r(v, 2) for v in jaw_y]
if mlc.device_type not in seen:
seq.append(_bld_item(mlc.device_type, mlc_values))
for t in mlc.jaw_types:
if t not in seen:
if t in JAW_X: seq.append(_bld_item(t, jaw_x))
elif t in JAW_Y: seq.append(_bld_item(t, jaw_y))
def build_plan(rp_template, beam_index, mlc, A, B, weights, total_mu, beam_type, plan_label):
new = copy.deepcopy(rp_template)
beam = new.BeamSequence[beam_index]
beam_number = int(beam.BeamNumber)
new.BeamSequence = Sequence([beam])
groups = []
for fg in new.FractionGroupSequence:
keep = [rb for rb in fg.ReferencedBeamSequence
if int(rb.ReferencedBeamNumber) == beam_number]
if not keep: continue
for rb in keep:
rb.BeamMeterset = _r(total_mu, 4)
for kw in ('BeamDose', 'BeamDoseSpecificationPoint'):
if hasattr(rb, kw): delattr(rb, kw)
fg.ReferencedBeamSequence = Sequence(keep)
fg.NumberOfBeams = len(keep)
groups.append(fg)
if not groups:
raise ValueError(f'FractionGroup 找不到 BeamNumber={beam_number}')
new.FractionGroupSequence = Sequence(groups)
now = datetime.datetime.now()
new.RTPlanDate = now.strftime('%Y%m%d')
new.RTPlanTime = now.strftime('%H%M%S')
new.RTPlanLabel = plan_label[:16]
new.SOPInstanceUID = generate_uid()
new.SeriesInstanceUID = generate_uid()
if p.new_study_uid: new.StudyInstanceUID = generate_uid()
new.ApprovalStatus = 'UNAPPROVED'
for kw in ('ReviewDate', 'ReviewTime', 'ReviewerName'):
if hasattr(new, kw): delattr(new, kw)
new.file_meta.MediaStorageSOPInstanceUID = new.SOPInstanceUID
new.file_meta.TransferSyntaxUID = ExplicitVRLittleEndian
orig_cps = beam.ControlPointSequence
cp0 = copy.deepcopy(orig_cps[0])
last_ref = getattr(orig_cps[-1], 'ReferencedDoseReferenceSequence', None)
last_ref = copy.deepcopy(last_ref) if last_ref is not None else None
n_cp = len(weights)
beam.BeamType = beam_type
beam.NumberOfControlPoints = n_cp
beam.FinalCumulativeMetersetWeight = 1.0
jaw_x = [p.field_center_x_mm - p.field_x_mm/2, p.field_center_x_mm + p.field_x_mm/2]
jaw_y = [p.field_center_y_mm - p.field_y_mm/2, p.field_center_y_mm + p.field_y_mm/2]
cp0.ControlPointIndex = 0
cp0.CumulativeMetersetWeight = 0.0
_fill_cp0_devices(cp0, mlc, list(A[0]) + list(B[0]), jaw_x, jaw_y)
cps = [cp0]
for k in range(1, n_cp):
cp = Dataset()
cp.ControlPointIndex = k
cp.CumulativeMetersetWeight = _r(weights[k], 6)
cp.BeamLimitingDevicePositionSequence = Sequence(
[_bld_item(mlc.device_type, list(A[k]) + list(B[k]))])
cps.append(cp)
if last_ref is not None and n_cp > 1:
cps[-1].ReferencedDoseReferenceSequence = last_ref
beam.ControlPointSequence = Sequence(cps)
return new
def save_plan(ds, path):
major = int(pydicom.__version__.split('.')[0])
if major >= 3:
ds.save_as(path, enforce_file_format=True)
else:
ds.is_little_endian = True
ds.is_implicit_VR = False
ds.save_as(path, write_like_original=False)
print('build_plan / save_plan 已定義')
6.1 寫出 Sliding Window 的 RT Plan¶
In [ ]:
rp_slide = build_plan(rp, beam_index, mlc,
A=A_slide, B=B_slide, weights=weights_slide,
total_mu=traj.total_mu, beam_type='DYNAMIC',
plan_label='MLC_SLIDING')
save_plan(rp_slide, p.out_sliding)
print(f'已寫出 {p.out_sliding}')
print(f" BeamType = {rp_slide.BeamSequence[0].BeamType}")
print(f" NumberOfControlPoints = {rp_slide.BeamSequence[0].NumberOfControlPoints}")
print(f" RTPlanLabel = {rp_slide.RTPlanLabel}")
6.2 寫出 Step-and-Shoot 的 RT Plan¶
In [ ]:
rp_sas = build_plan(rp, beam_index, mlc,
A=A_sas_cp, B=B_sas_cp, weights=weights_sas,
total_mu=p.total_mu, beam_type='STATIC',
plan_label='MLC_STEP')
save_plan(rp_sas, p.out_stepshoot)
print(f'已寫出 {p.out_stepshoot}')
print(f" BeamType = {rp_sas.BeamSequence[0].BeamType}")
print(f" NumberOfControlPoints = {rp_sas.BeamSequence[0].NumberOfControlPoints}")
print(f" RTPlanLabel = {rp_sas.RTPlanLabel}")
6.3 回讀驗證¶
In [ ]:
for path in (p.out_sliding, p.out_stepshoot):
ds = pydicom.dcmread(path)
b = ds.BeamSequence[0]
print(f'\n=== {path} ===')
print(f"BeamType={b.BeamType} CP數={b.NumberOfControlPoints}")
cps = b.ControlPointSequence
print(f" cp[0] Index={cps[0].ControlPointIndex} CumW={cps[0].CumulativeMetersetWeight}")
print(f" cp[1] Index={cps[1].ControlPointIndex} CumW={cps[1].CumulativeMetersetWeight}")
print(f" cp[-1] Index={cps[-1].ControlPointIndex} CumW={cps[-1].CumulativeMetersetWeight}")
mlc_item = [it for it in cps[0].BeamLimitingDevicePositionSequence
if str(it.RTBeamLimitingDeviceType).upper() == mlc.device_type][0]
print(f" cp[0] {mlc.device_type} 位置前 4 個 = {mlc_item.LeafJawPositions[:4]}")
7. 總結¶
| 項目 | Sliding Window | Step-and-Shoot |
|---|---|---|
BeamType |
DYNAMIC |
STATIC |
| 葉片運動 | 連續（CP 間線性內插） | 只在 segment 邊界跳動 |
| 每個 CP 的葉片位置 | 都不同 | 兩兩成對（同位置） |
| CP 數量 | 可任意加密 | 固定 2 × segment 數 |
| 誤差來源 | 取樣不足 | 量化層數 |
| MU 效率 | 較高 | 較低（segment 間 idle） |
| 演算法核心 | T_A − T_B = I |
1D 層級分解 Σ 1[L,R) = I |
延伸思考題¶
dose_rate_mu_min提高時，sliding 的誤差會如何變化？min_leaf_gap_mm增大時，兩種演算法各受什麼影響？- Step-and-shoot 為何需要 NNLS 修正？不做誤差差多少？
total_mu太小會發生什麼事？pixel_mm從 2.0 改 1.0 / 4.0 時，兩種演算法表現如何變？- 為什麼 sliding window 至少要一片葉片全速移動？
練習題¶
p.dose_rate_mu_min = 300再跑一次p.max_segments = 20，觀察 step-and-shoot 誤差變化p.invert = False，觀察 fluence mapp.min_leaf_gap_mm = 3.0，看哪種演算法受影響更大p.refine_weights = False，對比 4.5 的兩張圖