Xiao Furen
2ae08ac2cd
Refactor the preprocessing and segmentation pipeline to handle AP orientation variations and improve anatomical boundary detection. Key changes include: - Implement automated AP orientation detection in `process_single_image` to handle prone scans by flipping CT and labels when necessary. - Enhance `segment_spinous_process` using a gap-based approach to identify the spinal canal, providing more stable thresholds for spinous process and vertebral body segmentation. - Improve optimization search space by using the vertebral body (VBODY) projection for x/z bounding box calculation instead of the whole bone. - Refactor `render_bone_figure` to unify 2D/3D visualization and support detailed anatomical coloring (VBODY, spinous process). - Update `cl_score_torch_xfr` with more robust penalty handling for out-of-bone and null-voxel regions. - Add `retry_robust` utility to handle transient NFS file system errors. - Update `xfr_preprocess.py` to include anatomical segmentation coloring in rotated level visualizations.
1137 lines
49 KiB
Python
1137 lines
49 KiB
Python
import numpy as np
|
||
import SimpleITK as sitk
|
||
from scipy.ndimage import center_of_mass, rotate, distance_transform_edt, gaussian_filter
|
||
import matplotlib.pyplot as plt
|
||
|
||
# 鏡稱面法線的左右（LR，=x 軸）分量大下。正常左右鏡稱面法線幾乎沿 x；
|
||
# 若 |a| < 0.7，代表搜尋掉進前後（coronal）面局部極大（椎體是厚實塊狀，
|
||
# 易自鏡射拿高分，例：0019 L3/L4），以該平面定義的「中線帶」變成斜切
|
||
# 板層，其後側中線缺如判定不可信 → diagnose 回 bad_mirror_plane、
|
||
# 不判 no_spinous（避免假 [NO-SP] 誤改椎體切分）。
|
||
MIRROR_MIN_LR = 0.7
|
||
|
||
def _best_vertical_split(mask2d):
|
||
"""
|
||
binary 2D 陣列（最後一軸 = x）：找讓 mask 與其鏡射重疊最大的垂直線 x = t。
|
||
c[s] = Σ_u A[u]·A[s-u] 是每列自卷積；批次 FFT 一次算出所有 s。
|
||
回傳 (score, s_best)；s = 2·t（t 可能為半整數）。
|
||
"""
|
||
n = mask2d.shape[-1]
|
||
flat = np.asarray(mask2d).reshape(-1, n)
|
||
rows = flat[flat.max(axis=1) > 0]
|
||
if rows.size == 0:
|
||
return 0.0, 0
|
||
L = 1 << (2 * n - 1).bit_length()
|
||
padded = np.zeros((rows.shape[0], L), dtype=np.float32)
|
||
padded[:, :n] = rows.astype(np.float32)
|
||
F = np.fft.rfft(padded, axis=1)
|
||
conv = np.fft.irfft(F * F, axis=1)[:, :2 * n - 1]
|
||
score = np.clip(conv.sum(axis=0), 0, None)
|
||
s_best = int(np.argmax(score))
|
||
return float(score[s_best]), s_best
|
||
|
||
|
||
def best_symmetry_axis_angle(proj, coarse_step=10.0, fine_step=1.0):
|
||
"""
|
||
2D binary 投影（axis0=row y, axis1=col x）的左右鏡稱軸搜尋：
|
||
將影像旋轉 θ 後，鏡稱軸變成垂直線（col = const），用 _best_vertical_split 評分；
|
||
粗搜 0..180°（coarse_step）再在 winner 附近細搜（fine_step）。
|
||
回傳 (theta_deg, score, s_best)。
|
||
"""
|
||
base = (np.asarray(proj) > 0).astype(np.float32)
|
||
best_th, best, best_s = 0.0, -1.0, 0
|
||
for th in np.arange(0.0, 180.0, coarse_step):
|
||
sc, s = _best_vertical_split(rotate(base, th, reshape=False, order=0) > 0.5)
|
||
if sc > best:
|
||
best_th, best, best_s = float(th), sc, s
|
||
for th in np.arange(best_th - coarse_step, best_th + coarse_step, fine_step):
|
||
th2 = float(th) % 180.0
|
||
sc, s = _best_vertical_split(rotate(base, th2, reshape=False, order=0) > 0.5)
|
||
if sc > best:
|
||
best_th, best, best_s = th2, sc, s
|
||
return best_th, best, best_s
|
||
|
||
|
||
def _mirror_hit_count(X, Y, Z, n, d, m, sz, sy, sx):
|
||
"""每個骨 voxel 對平面 n·p=d 鏡射後四捨五入到最近 voxel，回傳命中骨 voxel 的數量"""
|
||
nx, ny, nz = n
|
||
dist = X * nx + Y * ny + Z * nz - d
|
||
rx = (X - 2.0 * dist * nx).round().astype(np.int32)
|
||
ry = (Y - 2.0 * dist * ny).round().astype(np.int32)
|
||
rz = (Z - 2.0 * dist * nz).round().astype(np.int32)
|
||
ok = (rx >= 0) & (rx < sx) & (ry >= 0) & (ry < sy) & (rz >= 0) & (rz < sz)
|
||
if not ok.any():
|
||
return 0
|
||
return int(m[rz[ok], ry[ok], rx[ok]].sum())
|
||
|
||
|
||
def best_symmetry_plane(mask_zyx, phi_max=45.0, subsample=7,
|
||
coarse_step=7.5, fine_step=1.0):
|
||
"""
|
||
3D bone mask (z, y, x) 的最佳鏡稱面，一般平面方程 a·x + b·y + c·z = d
|
||
（voxel index 座標；(a,b,c) 為單位法線，方向任意，不限制平行 YZ 面）。
|
||
以「每個 bone voxel 鏡射後四捨五入到最近 voxel 的命中率」為分數，
|
||
參數化 n = R_y(phi)·R_z(theta)·(1,0,0)：
|
||
theta : 法線在 axial (x,y) 平面內的旋轉
|
||
phi : 法線出平面的傾斜，限制在 [-phi_max, +phi_max] 以確保仍切分左右
|
||
搜尋：theta 由 2D axial 投影粗定位，再 3D 三段式（粗→細→微細）
|
||
（粗/細階段用等距子樣 voxel 加速）。
|
||
回傳 dict：
|
||
plane : (a, b, c, d) -> a*x + b*y + c*z = d
|
||
normal : (a, b, c)
|
||
offset : d
|
||
theta_deg / phi_deg : 法線參數
|
||
ratio : 鏡射命中率（0..1）
|
||
u, v : 平面內兩個正交方向（供繪製用）
|
||
"""
|
||
m = np.asarray(mask_zyx) > 0
|
||
sz, sy, sx = m.shape
|
||
zz, yy, xx = np.nonzero(m)
|
||
if xx.size < 100:
|
||
n = np.array([1.0, 0.0, 0.0])
|
||
d = (sx - 1) / 2.0
|
||
return {'plane': (1.0, 0.0, 0.0, float(d)), 'normal': (1.0, 0.0, 0.0),
|
||
'offset': float(d), 'theta_deg': 0.0, 'phi_deg': 0.0,
|
||
'ratio': 1.0, 'u': (0.0, 1.0, 0.0), 'v': (0.0, 0.0, 1.0)}
|
||
xf = xx.astype(np.float32)
|
||
yf = yy.astype(np.float32)
|
||
zf = zz.astype(np.float32)
|
||
N = xf.size
|
||
c0 = np.array([(sx - 1) / 2.0, (sy - 1) / 2.0, (sz - 1) / 2.0])
|
||
Xs, Ys, Zs = xf[::subsample], yf[::subsample], zf[::subsample]
|
||
|
||
def n_of(theta_deg, phi_deg):
|
||
t = np.deg2rad(theta_deg)
|
||
p = np.deg2rad(phi_deg)
|
||
return np.array([np.cos(t) * np.cos(p), np.sin(t), -np.cos(t) * np.sin(p)])
|
||
|
||
# theta 初值：2D axial 投影搜尋（垂直面情形）
|
||
proj = m.max(axis=0)
|
||
theta0 = best_symmetry_axis_angle(proj.astype(np.float32))[0] if proj.sum() >= 10 else 0.0
|
||
|
||
best = None
|
||
|
||
def consider(theta, phi, d, full=False):
|
||
nonlocal best
|
||
n = n_of(theta, phi)
|
||
X, Y, Z = (xf, yf, zf) if full else (Xs, Ys, Zs)
|
||
s = _mirror_hit_count(X, Y, Z, n, d, m, sz, sy, sx)
|
||
if best is None or s > best['score']:
|
||
best = {'score': s, 'theta': float(theta), 'phi': float(phi), 'd': float(d)}
|
||
|
||
# stage 1：粗搜尋（子樣）
|
||
for th in np.arange(theta0 - coarse_step * 2, theta0 + coarse_step * 2 + 1e-9, coarse_step):
|
||
for ph in np.arange(-phi_max, phi_max + 1e-9, coarse_step):
|
||
n = n_of(th, ph)
|
||
d0 = float(n @ c0)
|
||
for dd in (-10.0, 0.0, 10.0):
|
||
consider(th, ph, d0 + dd)
|
||
# stage 2：細搜尋（子樣）
|
||
for th in np.arange(best['theta'] - 2 * fine_step * 2, best['theta'] + 2 * fine_step * 2 + 1e-9, fine_step):
|
||
for ph in np.arange(best['phi'] - 2 * fine_step * 2, best['phi'] + 2 * fine_step * 2 + 1e-9, fine_step):
|
||
for dd in (-3.0, -1.0, 0.0, 1.0, 3.0):
|
||
consider(th, ph, best['d'] + dd)
|
||
# stage 3：微細搜尋（全分辨率）
|
||
for th in np.arange(best['theta'] - fine_step, best['theta'] + fine_step + 1e-9, fine_step / 2):
|
||
for ph in np.arange(best['phi'] - fine_step, best['phi'] + fine_step + 1e-9, fine_step / 2):
|
||
for dd in (-1.0, -0.5, 0.0, 0.5, 1.0):
|
||
consider(th, ph, best['d'] + dd, full=True)
|
||
# stage 4：d 微調（全分辨率）
|
||
b = best
|
||
for dd in np.arange(-1.0, 1.0 + 1e-9, 0.25):
|
||
consider(b['theta'], b['phi'], b['d'] + dd, full=True)
|
||
|
||
n = n_of(best['theta'], best['phi'])
|
||
theta = best['theta']
|
||
d = best['d']
|
||
if theta > 90.0 or theta < -90.0: # 正規化到 (-90, 90]，法線翻轉時 d 取反
|
||
theta -= 180.0
|
||
n = -n
|
||
d = -d
|
||
ratio = _mirror_hit_count(xf, yf, zf, n, d, m, sz, sy, sx) / N
|
||
u = np.cross(n, [0.0, 0.0, 1.0])
|
||
if np.linalg.norm(u) < 0.1:
|
||
u = np.cross(n, [1.0, 0.0, 0.0])
|
||
u = u / np.linalg.norm(u)
|
||
v = np.cross(n, u)
|
||
return {
|
||
'plane': (float(n[0]), float(n[1]), float(n[2]), float(d)),
|
||
'normal': (float(n[0]), float(n[1]), float(n[2])),
|
||
'offset': float(d),
|
||
'theta_deg': float(theta),
|
||
'phi_deg': float(best['phi']),
|
||
'ratio': float(ratio),
|
||
'u': (float(u[0]), float(u[1]), float(u[2])),
|
||
'v': (float(v[0]), float(v[1]), float(v[2])),
|
||
}
|
||
|
||
def diagnose_spinous_process(mask_zyx, sym, band_frac=0.06, min_band=6.0,
|
||
deficit_min=4.0, rear_margin=3.0, rear3_max=20,
|
||
narrow_frac=0.4):
|
||
"""
|
||
診斷該椎體的棘突（中線後側構造）是否缺如
|
||
（先前手術如 laminectomy / 棘突切除所造成）。
|
||
|
||
棘突完整時，其尖端是全椎體最後側的骨構造且位於中線，因此：
|
||
deficit = 全骨最後側 AP - 中線帶最後側 AP ~= 0，
|
||
且中線窄帶 |s| <= narrow_frac*w 在最後側骨 rear_margin 內有充足骨 voxel。
|
||
切除後，中線帶的後側止於殘留後側要素，最後側骨偏到側方的
|
||
殘餘結構（椎板/關節突）→ deficit 大 + 後側中線空洞。
|
||
判定：no_spinous = deficit >= deficit_min 且 rear3 <= rear3_max。
|
||
|
||
回傳 dict：
|
||
n_bone : 骨 voxel 總數
|
||
band_w : 中線帶半寬 w（voxel）
|
||
deficit : 後側中線缺如量（voxel，0.5 mm/vox）
|
||
rear3 : 最後側 3 voxel 內、中線窄帶的骨 voxel 數
|
||
rear6 : 最後側 6 voxel 內、中線窄帶的骨 voxel 數
|
||
top_off : axial 投影（y,x）最後側一行相對鏡稱中線的偏移（voxel）
|
||
no_spinous : True = 後側中線缺如（棘突已切除）
|
||
reason : 'ok' / 'no_spinous' / 'too_few_bone' / 'too_few_midline'
|
||
/ 'bad_mirror_plane'
|
||
"""
|
||
m = np.asarray(mask_zyx) > 0
|
||
zz, yy, xx = np.nonzero(m)
|
||
out = {'n_bone': int(zz.size), 'band_w': None, 'deficit': 0.0,
|
||
'rear3': None, 'rear6': None, 'top_off': None,
|
||
'no_spinous': False, 'reason': 'ok'}
|
||
if zz.size < 100:
|
||
out['reason'] = 'too_few_bone'
|
||
return out
|
||
a, b, c, d = sym['plane']
|
||
if abs(a) < MIRROR_MIN_LR:
|
||
# 鏡稱面非左右為主（前後 coronal 局部極大）→ 缺如判定不可信，
|
||
# 按 ok 回傳（寧可漏報，不可假切除）
|
||
out['reason'] = 'bad_mirror_plane'
|
||
return out
|
||
X = xx.astype(np.float64)
|
||
Y = yy.astype(np.float64)
|
||
Z = zz.astype(np.float64)
|
||
s = X * a + Y * b + Z * c - d
|
||
u_ap = _ap_axis(sym) # +AP = 後側（y 小側）
|
||
ap = X * u_ap[0] + Y * u_ap[1] + Z * u_ap[2]
|
||
w = max(float(min_band), float(band_frac) * float(s.max() - s.min()))
|
||
mid = np.abs(s) <= w
|
||
aps = ap[mid]
|
||
out['band_w'] = float(w)
|
||
if aps.size < 50:
|
||
out['reason'] = 'too_few_midline'
|
||
return out
|
||
narrow = np.abs(s) <= w * float(narrow_frac)
|
||
ap_max_all = float(ap.max())
|
||
deficit = ap_max_all - float(aps.max())
|
||
rear3 = int(((ap >= ap_max_all - rear_margin) & narrow).sum())
|
||
rear6 = int(((ap >= ap_max_all - 2.0 * rear_margin) & narrow).sum())
|
||
out['deficit'] = float(deficit)
|
||
out['rear3'] = rear3
|
||
out['rear6'] = rear6
|
||
proj = m.max(axis=0) # (y, x)
|
||
ys_p, xs_p = np.where(proj > 0)
|
||
# 只在鏡稱面 x 分量主導時才算 top_off（最後側一行對中線的 x 偏移）；
|
||
# 脊椎在面內大角度旋轉時 |a| 小，(d-b*y-c*z)/a 會放大成無意義的巨值。
|
||
if ys_p.size and abs(a) >= 0.5 * max(abs(b), abs(c)):
|
||
y_min = int(ys_p.min())
|
||
zc = (m.shape[0] - 1) / 2.0
|
||
x_mid = (d - b * y_min - c * zc) / a
|
||
out['top_off'] = float(abs(xs_p[ys_p == y_min].mean() - x_mid))
|
||
if deficit >= float(deficit_min) and rear3 <= int(rear3_max):
|
||
out['no_spinous'] = True
|
||
out['reason'] = 'no_spinous'
|
||
return out
|
||
|
||
|
||
def segment_spinous_process(mask_zyx, sym, band_frac=0.06, min_band=6.0,
|
||
min_mass_frac=0.05, expand_cap=5):
|
||
"""
|
||
以 best_symmetry_plane 的結果 sym 從 3D bone mask (z, y, x) 切出棘突。
|
||
棘突是中線後側構造，利用鏡稱面 a·x+b·y+c·z=d 定義：
|
||
1) 中線帶：骨 voxel 到平面的有號距離 |s| <= w，
|
||
w = max(min_band, band_frac * s 全寬)
|
||
2) 前後方向：平面內兩軸 (u, v) 中 |y| 分量大者，
|
||
正規化成 +AP = 後側（本資料系 y 往前遞增，後側 = y 小側）
|
||
3) 中線帶的 AP 分佈呈兩大叢（椎體在前、椎弓/棘突在後），中間椎管
|
||
空隙（連續 <5% 峰值的安靜 bin）以「最長空隙」定位（= 椎管）：
|
||
棘突閾值取空隙後側端（AP >= 空隙後側端 = 後側叢，含中線椎弓）；
|
||
回傳的椎體閾值取空隙體側邊谷（供 segment_vertebral_body 使用）；
|
||
兩閾值被椎管隔開，椎體後側緣 fringe 不會被誤判成棘突。
|
||
無明顯空隙（如骨橋）fallback 取中線帶後側 15%。
|
||
4) 中線帶側緣補回（_expand_spinous_runs，expand_cap）：帶由鏡稱面定義，
|
||
棘突楔若略偏中線，側緣薄條會留在帶外成為 other bone；每條 (y, z)
|
||
線把棘突 run 向左右各補至多 expand_cap 個 bone voxel。
|
||
先做後側中線缺如判定（diagnose_spinous_process，laminectomy/棘突切除）；
|
||
缺如、或谷底後側叢只剩 <5% 殘片時 mode='no_spinous'、sp_mask=None，
|
||
由呼叫端（入口面 / 椎體）改用對應策略。
|
||
回傳 (sp_mask (z,y,x) bool, ap_thresh, info dict)；
|
||
資料過少或棘突缺如時 sp_mask = None（info['mode'] 說明原因）。
|
||
"""
|
||
m = np.asarray(mask_zyx) > 0
|
||
zz, yy, xx = np.nonzero(m)
|
||
info = {'n_bone': int(zz.size), 'n_sp': 0, 'ap_thresh': None,
|
||
'band_w': None, 'mode': 'empty'}
|
||
if zz.size < 50:
|
||
return None, None, info
|
||
a, b, c, d = sym['plane']
|
||
X = xx.astype(np.float64)
|
||
Y = yy.astype(np.float64)
|
||
Z = zz.astype(np.float64)
|
||
s = X * a + Y * b + Z * c - d
|
||
u = np.array(sym['u'])
|
||
v = np.array(sym['v'])
|
||
u_ap = u if abs(u[1]) >= abs(v[1]) else v
|
||
if u_ap[1] > 0:
|
||
u_ap = -u_ap # +AP = 後側（y 小側）
|
||
w = max(float(min_band), float(band_frac) * float(s.max() - s.min()))
|
||
mid = np.abs(s) <= w
|
||
ap = X * u_ap[0] + Y * u_ap[1] + Z * u_ap[2]
|
||
aps = ap[mid]
|
||
info['band_w'] = float(w)
|
||
if aps.size < 50:
|
||
info['mode'] = 'too_few_midline'
|
||
return None, None, info
|
||
|
||
# 後側中線缺如（laminectomy / 棘突切除）：最後側骨偏到側方殘餘、
|
||
# 中線帶後側無質量，谷底搜尋沒有意義 → no_spinous，
|
||
# 呼叫端（y_indices 入口面 / segment_vertebral_body）改用對應策略。
|
||
diag = diagnose_spinous_process(mask_zyx, sym, band_frac=band_frac, min_band=min_band)
|
||
info['deficit'] = diag['deficit']
|
||
info['rear3'] = diag['rear3']
|
||
info['rear6'] = diag['rear6']
|
||
if diag['no_spinous']:
|
||
info['mode'] = 'no_spinous'
|
||
info['top_off'] = diag['top_off']
|
||
return None, None, info
|
||
|
||
lo = int(np.floor(aps.min()))
|
||
hi = int(np.ceil(aps.max()))
|
||
th = None
|
||
th_sp = None
|
||
mode = 'fallback'
|
||
if hi - lo >= 10:
|
||
hist, edges = np.histogram(aps, bins=range(lo, hi + 1))
|
||
csum = np.concatenate([[0], np.cumsum(hist)])
|
||
total = csum[-1]
|
||
peak = hist.max()
|
||
# 收集連續安靜 bin（<5% 峰值）的「空隙」，要求兩側質量都夠
|
||
# （>= min_mass_frac * total），取最長者 = 椎管。
|
||
# （舊式單 bin score=min(前,後) 最大化：後側叢質量較小時恆落在
|
||
# 空隙的椎體側第一安靜 bin，椎體後側緣被傾斜鏡稱帶斜切出的
|
||
# 1~2 體素 fringe 會 >= 該閾值而誤判成棘突 → 圖上椎體內出現
|
||
# 棘突點。）
|
||
quiet = hist < 0.05 * peak
|
||
runs = []
|
||
i = 0
|
||
while i < len(hist):
|
||
if not quiet[i]:
|
||
i += 1
|
||
continue
|
||
j = i
|
||
while j + 1 < len(hist) and quiet[j + 1]:
|
||
j += 1
|
||
left = int(csum[i])
|
||
right = int(total - csum[j + 1])
|
||
if left >= min_mass_frac * total and right >= min_mass_frac * total:
|
||
runs.append((j - i + 1, min(left, right), i, j))
|
||
i = j + 1
|
||
if runs:
|
||
# 最長空隙勝（平手取兩側質量大者）；正常椎管是最長安靜區間
|
||
runs.sort(key=lambda r: (r[0], r[1]), reverse=True)
|
||
_, _, i0, i1 = runs[0]
|
||
post_frac = (total - csum[i1 + 1]) / total
|
||
if post_frac < 0.05 and abs(a) >= MIRROR_MIN_LR:
|
||
# 空隙後側叢只剩小殘片（<5% 中線帶質量）＝棘突幾乎全除
|
||
# （部分切除殘餘）：當作缺如，椎體切分不採用此閾值。
|
||
# 僅在鏡稱面左右為主時才算數（斜板層時 post_frac 不可信）
|
||
info['mode'] = 'no_spinous'
|
||
info['post_frac'] = float(post_frac)
|
||
return None, None, info
|
||
th = float(0.5 * (edges[i0] + edges[i0 + 1])) # 回傳值：體側邊谷，椎體 AP 切點
|
||
th_sp = float(edges[i1 + 1]) # 空隙後側端：棘突由此開始
|
||
mode = 'gap'
|
||
if th is None:
|
||
th = float(np.quantile(aps, 0.85))
|
||
sp_mask = np.zeros(m.shape, dtype=bool)
|
||
# 棘突用空隙後側端 th_sp（無空隙 fallback 時退回 quantile th）；
|
||
# 椎體切點（回傳 th）維持身體側，兩者在椎管兩端、mask 不相觸。
|
||
sel = mid & (ap >= (th_sp if th_sp is not None else th))
|
||
sp_mask[zz[sel], yy[sel], xx[sel]] = True
|
||
sp_mask = _expand_spinous_runs(sp_mask, m, cap=expand_cap)
|
||
info.update(n_sp=int(sp_mask.sum()), ap_thresh=th, ap_thresh_sp=th_sp, mode=mode)
|
||
return sp_mask, th, info
|
||
|
||
|
||
def _expand_spinous_runs(sp_mask, bone, cap=5):
|
||
"""對每條 (y, z) 水平線（x = 左右方向），把既有棘突 run 向左右兩側各延伸
|
||
至多 cap 個 bone voxel：補回落在中線帶 |s|<=w 之外的棘突楔側緣薄條
|
||
（帶由鏡稱面定義，楔略偏中線時側緣會被漏成 other bone）。
|
||
遇到非 bone、既有棘突 voxel 或達 cap 即停；只延伸既有 run，不新建
|
||
run（該線無棘突時絕不標記）。椎弓厚處最多少數 cap 個「棘突基部」
|
||
相鄰 bone 一併補入，屬棘突連續結構。"""
|
||
out = sp_mask.copy()
|
||
if not out.any():
|
||
return out
|
||
nz, ny, nx = out.shape
|
||
for y in range(ny):
|
||
sp_y = out[:, y, :]
|
||
if not sp_y.any():
|
||
continue
|
||
bone_y = bone[:, y, :]
|
||
for z in np.flatnonzero(sp_y.any(axis=1)):
|
||
idx = np.flatnonzero(sp_y[z])
|
||
runs = np.split(idx, np.flatnonzero(np.diff(idx) > 1) + 1)
|
||
for r in runs:
|
||
lo, hi = int(r[0]), int(r[-1])
|
||
x = lo - 1
|
||
added = 0
|
||
while x >= 0 and added < cap and bone_y[z, x] and not out[z, y, x]:
|
||
out[z, y, x] = True
|
||
x -= 1
|
||
added += 1
|
||
x = hi + 1
|
||
added = 0
|
||
while x < nx and added < cap and bone_y[z, x] and not out[z, y, x]:
|
||
out[z, y, x] = True
|
||
x += 1
|
||
added += 1
|
||
return out
|
||
|
||
def smooth_mask_sdf(mask_zyx, sigma=1.5):
|
||
"""對 3D boolean mask 做有號距離場（SDF）高斯平滑：
|
||
把表面沿法線做近似平均曲率流（mean curvature flow），圓化凸出的
|
||
鋸齒/尖刺、填進凹陷的缺口，但不會均勻收縮整體形状（與直接對 mask
|
||
做高斯Blur再 thresh 不同——後者會把薄結構一起磨細、體積明顯縮小）。
|
||
這是消除「三线性旋轉 > 0.5 之後仍剩餘的 voxel 級鋸齒」最有效的方法。
|
||
|
||
參數：
|
||
mask_zyx : (z, y, x) boolean / 0-1 陣列
|
||
sigma : 平滑強度（voxel）。1.5 ≈ 表面積 -4%、體積 -0.6%；
|
||
越大越圓，但薄結構（棘突尖）會被磨損越多。
|
||
回傳同 shape 的 boolean mask。
|
||
"""
|
||
m = np.asarray(mask_zyx) > 0
|
||
if not m.any() or m.all():
|
||
return m.copy()
|
||
d_in = distance_transform_edt(m)
|
||
d_out = distance_transform_edt(~m)
|
||
sdf = d_in - d_out
|
||
return gaussian_filter(sdf, sigma=float(sigma)) > 0.0
|
||
|
||
def shape_based_smooth(vol01, sigma=1.5):
|
||
"""Shape-based interpolation（SBI）平滑：以 0..1 體積的 0.5 等值面為
|
||
物件形狀（shape），計算有號距離場（SDF，內正外負），在距離域做高斯
|
||
插值（interpolation），再於 0 重新閾值，得到圓滑的 smooth mask。
|
||
與 smooth_mask_sdf 的差異：邊界層（0 < v < 1）內三线性場值滿足
|
||
v ≈ 0.5 + d（d = 以 voxel 為單位的有號距離），故直接以 v - 0.5 替換
|
||
EDT 量值，得到次體素精確的 SDF（層外仍由 distance_transform_edt 給出
|
||
精確距離）；表面梯度連續，不再是被 EDT 量化成 ±1 的硬階梯，平滑後
|
||
的邊界更貼合 soft boundary 的真實位置。
|
||
參數：
|
||
vol01 : (z, y, x) 0..1 浮點體積（例如旋轉後的 _binary_linear）
|
||
sigma : 高斯平滑強度（voxel），意義與 smooth_mask_sdf 相同。
|
||
回傳：同 shape 的 boolean mask。
|
||
"""
|
||
v = np.clip(np.asarray(vol01, dtype=np.float32), 0.0, 1.0)
|
||
b = v > 0.5
|
||
if not b.any() or b.all():
|
||
return b.copy()
|
||
d_in = distance_transform_edt(b)
|
||
d_out = distance_transform_edt(~b)
|
||
sdf = d_in - d_out
|
||
band = (v > 0) & (v < 1)
|
||
sdf[band] = v[band] - 0.5
|
||
return gaussian_filter(sdf, sigma=float(sigma)) > 0.0
|
||
|
||
def _ap_axis(sym):
|
||
"""鏡稱面內的前後（AP）軸：(u,v) 中 |y| 分量大者，正規化成 +AP = 後側（y 小側）
|
||
（本資料系 y 往前遞增，與 segment_spinous_process 同慣例）。"""
|
||
u = np.array(sym['u'])
|
||
v = np.array(sym['v'])
|
||
u_ap = u if abs(u[1]) >= abs(v[1]) else v
|
||
if u_ap[1] > 0:
|
||
u_ap = -u_ap
|
||
return u_ap
|
||
|
||
|
||
def anterior_y_side(mask_zyx, band_frac=0.06, min_band=6.0, min_ratio=1.3,
|
||
min_side_frac=0.05, min_voxels=100):
|
||
"""判定 index 系 (z,y,x) 中椎體前側（椎體塊所在側）在哪個 y 端：
|
||
'y_max' : 前側 = y 大側（與本套件各函式預設慣例一致：前 = y 大、後 = y 小）
|
||
'y_min' : 前側 = y 小側（前後翻轉，例：prone 伏位掃描個案；
|
||
呼叫端應把 CT 與 label 在 y 方向翻轉，使輸出方向與其他個案一致）
|
||
None : 無法判定（骨量太少、鏡稱面非左右為主（前後方向局部極大，
|
||
同 0019 L3/L4 情形）、無明顯椎管安靜區間、或前後質量差不夠）。
|
||
無法判定時呼叫端維持預設方向（寧可不翻轉，不誤翻轉）。
|
||
|
||
原理：先取最佳鏡稱面（best_symmetry_plane——左右對稱構造，其結果不受
|
||
前後翻轉影響），取該面的中線帶（與 segment_spinous_process 同 band 定義），
|
||
在帶內計算前後坐標（面內 |y| 分量大者為 AP 軸；此處用無向版本、
|
||
指向 y 大側），做 AP 直方圖：單椎體在帶內的 AP 分佈有兩大叢
|
||
（前側椎體塊、後側棘突/椎板），以椎管（最長安靜區間，與
|
||
segment_spinous_process 同一套閾值）分隔；椎體塊質量恆明顯大於
|
||
棘突/椎板（量測：正常 case ratio 1.4~2.5），質量大側 = 前側。
|
||
"""
|
||
m = np.asarray(mask_zyx) > 0
|
||
if int(m.sum()) < int(min_voxels):
|
||
return None
|
||
# 裁到骨頭 bbox：輸入若是整顆 volume（如 0.5mm 重取樣 label 的單層
|
||
# mask），best_symmetry_plane 的初始 search 中心 c0 = 體積盒中心會偏離
|
||
# 椎體；裁切後 c0 落在椎體上（與各輸出 bbox 裁切遮罩同條件）。
|
||
# 純平移不影響 y 端方向判定。
|
||
zz0, yy0, xx0 = np.nonzero(m)
|
||
z0, z1 = int(zz0.min()), int(zz0.max())
|
||
y0, y1 = int(yy0.min()), int(yy0.max())
|
||
x0, x1 = int(xx0.min()), int(xx0.max())
|
||
m = m[z0:z1 + 1, y0:y1 + 1, x0:x1 + 1]
|
||
try:
|
||
sym = best_symmetry_plane(m)
|
||
except Exception:
|
||
return None
|
||
if abs(sym['normal'][0]) < MIRROR_MIN_LR:
|
||
# 鏡稱面非左右為主（前後 coronal 局部極大）→ 中線帶失效，無法判定前後
|
||
return None
|
||
a, b, c, d = sym['plane']
|
||
zz, yy, xx = np.nonzero(m)
|
||
X = xx.astype(np.float64)
|
||
Y = yy.astype(np.float64)
|
||
Z = zz.astype(np.float64)
|
||
s = X * a + Y * b + Z * c - d
|
||
w = max(float(min_band), float(band_frac) * float(s.max() - s.min()))
|
||
band = np.abs(s) <= w
|
||
Xb, Yb, Zb = X[band], Y[band], Z[band]
|
||
if Xb.size < int(min_voxels):
|
||
return None
|
||
u = np.array(sym['u'])
|
||
v = np.array(sym['v'])
|
||
ap = u if abs(u[1]) >= abs(v[1]) else v
|
||
if ap[1] < 0:
|
||
ap = -ap # 無向：指向 y 大側
|
||
aproj = Xb * ap[0] + Yb * ap[1] + Zb * ap[2]
|
||
lo = int(np.floor(aproj.min()))
|
||
hi = int(np.ceil(aproj.max()))
|
||
if hi - lo < 10:
|
||
return None
|
||
hist, edges = np.histogram(aproj, bins=range(lo, hi + 1))
|
||
csum = np.concatenate([[0], np.cumsum(hist)])
|
||
total = csum[-1]
|
||
peak = hist.max()
|
||
# 最長安靜區間（<5% 峰值，兩側各 >= min_side_frac 質量）= 椎管
|
||
quiet = hist < 0.05 * peak
|
||
runs = []
|
||
i = 0
|
||
while i < len(hist):
|
||
if not quiet[i]:
|
||
i += 1
|
||
continue
|
||
j = i
|
||
while j + 1 < len(hist) and quiet[j + 1]:
|
||
j += 1
|
||
left = int(csum[i])
|
||
right = int(total - csum[j + 1])
|
||
if left >= min_side_frac * total and right >= min_side_frac * total:
|
||
runs.append((j - i + 1, min(left, right), i, j))
|
||
i = j + 1
|
||
if not runs:
|
||
return None
|
||
runs.sort(key=lambda r: (r[0], r[1]), reverse=True)
|
||
_, _, i0, i1 = runs[0]
|
||
m_low = int(csum[i0]) # 空隙 y 小側叢質量
|
||
m_high = int(total - csum[i1 + 1]) # 空隙 y 大側叢質量
|
||
if min(m_low, m_high) < min_side_frac * total:
|
||
return None
|
||
ratio = max(m_low, m_high) / float(max(1, min(m_low, m_high)))
|
||
if ratio < float(min_ratio):
|
||
return None
|
||
return 'y_max' if m_high > m_low else 'y_min'
|
||
|
||
|
||
def _full_ap_valley(aps, min_side_frac=0.15, max_ratio=0.85, smooth=3):
|
||
"""終板下骨體 AP 分佈的平滑谷底：最深相對谷底（sm[i] 對鄰近峰的最小比值），
|
||
要求兩側各有 >= min_side_frac 的質量（拒絕對小尾巴的偽谷底）。
|
||
回傳 (th 或 None, ratio 或 None)。"""
|
||
lo = int(np.floor(aps.min()))
|
||
hi = int(np.ceil(aps.max()))
|
||
if hi - lo < 10:
|
||
return None, None
|
||
hist, edges = np.histogram(aps, bins=range(lo, hi + 1))
|
||
sm = np.convolve(hist, np.ones(smooth) / float(smooth), mode='same')
|
||
csum = np.concatenate([[0], np.cumsum(hist)])
|
||
total = csum[-1]
|
||
best_i, best_ratio = None, 1.0
|
||
for i in range(1, len(hist) - 1):
|
||
if not (sm[i] <= sm[i - 1] and sm[i] <= sm[i + 1]):
|
||
continue
|
||
if csum[i] < min_side_frac * total or (total - csum[i + 1]) < min_side_frac * total:
|
||
continue
|
||
flank = min(float(sm[:i].max()), float(sm[i + 1:].max()))
|
||
if flank <= 0:
|
||
continue
|
||
ratio = float(sm[i]) / flank
|
||
if ratio <= max_ratio and ratio < best_ratio:
|
||
best_ratio, best_i = ratio, i
|
||
if best_i is None:
|
||
return None, None
|
||
return float(0.5 * (edges[best_i] + edges[best_i + 1])), best_ratio
|
||
|
||
|
||
def _posterior_min_threshold(aps, rear_frac=0.40, min_side_frac=0.08, smooth=3):
|
||
"""AP 分佈後側 `rear_frac` 區間內、平滑直方圖的最小值位置（該處後側質量
|
||
比例需 >= min_side_frac）：棘突缺如椎體中 _full_ap_valley 失敗時的
|
||
「椎體後側末端 = 體/弓最薄處」備援。
|
||
回傳 th 或 None。"""
|
||
lo = int(np.floor(aps.min()))
|
||
hi = int(np.ceil(aps.max()))
|
||
if hi - lo < 10:
|
||
return None
|
||
hist, edges = np.histogram(aps, bins=range(lo, hi + 1))
|
||
sm = np.convolve(hist, np.ones(smooth) / float(smooth), mode='same')
|
||
csum = np.concatenate([[0], np.cumsum(hist)])
|
||
total = csum[-1]
|
||
i0 = int(len(hist) * (1.0 - rear_frac))
|
||
best_i, best_val = None, np.inf
|
||
for i in range(i0, len(hist)):
|
||
if (total - csum[i + 1]) < min_side_frac * total:
|
||
continue
|
||
if sm[i] < best_val:
|
||
best_val, best_i = float(sm[i]), i
|
||
if best_i is None:
|
||
return None
|
||
return float(0.5 * (edges[best_i] + edges[best_i + 1]))
|
||
|
||
|
||
def segment_vertebral_body(mask_zyx, sym, endplate, ap_thresh, sp_mode, margin=2.0,
|
||
sliver_frac=0.20, lat_margin=3.0):
|
||
"""
|
||
以兩平面從 3D bone mask (z, y, x) 切出椎體（前側中央主體塊）：
|
||
1) best_upper_endplate_plane 的上終板面（法線朝上 a·x+b·y+c·z=d）：
|
||
只保留終板下側（身體側，e·p-d <= margin）的 bone voxel，
|
||
排除跨在終板上方的後側構造（椎板/棘突）。
|
||
2) best_symmetry_plane 的鏡稱面：提供平面內 AP 軸（_ap_axis，
|
||
+AP = 後側，與 segment_spinous_process 同慣例）；AP 切點取三層優先：
|
||
a) sp_mode == 'gap'：中線帶（椎管）的體/弓谷底 ap_thresh
|
||
（segment_spinous_process 回傳值，最可靠）；
|
||
b) sp_mode == 'no_spinous'（棘突已切除）：中線帶沒有空的體/弓谷底
|
||
可用；椎體後側末端 = 終板下整體 AP 分佈的局部最小（體/弓最薄處）。
|
||
切除後殘留後側要素會部分填滿椎管、谷底比正常椎體淺，
|
||
因此放寬 _full_ap_valley 條件（min_side_frac 0.15->0.10、
|
||
max_ratio 0.85->0.95）；再失敗則取後側 40% 區間的平滑谷底
|
||
（_posterior_min_threshold）。
|
||
c) 否則：終板下整體 AP 分佈的平滑谷底（_full_ap_valley，
|
||
處理中線骨橋等中線搜尋 fallback 的情形）；
|
||
d) 最後 fallback：AP 分佈 55 百分位（可能切進椎體內，會打 WARNING）。
|
||
椎體 = AP < 切點（切點之前側）且終板下側的 bone。
|
||
3) 側向包絡（lateral clip）：椎體是終板下的中央塊，橫突（及弓根
|
||
側緣）向鏡稱面法線方向延伸，其前緣恰好跨過 AP 切點（椎體後側
|
||
兩角處），只靠兩平面會把橫突前段算進椎體。在「離 AP 切點較遠」
|
||
的 AP 窗（排除靠切點後側 sliver_frac 的帶，帶內正是橫突前緣）
|
||
量出椎體自身的側向（鏡稱面）寬度，再把候選裁到該寬度
|
||
+ margin；橫突（遠超出椎體寬多達數厘米）被去除，椎體本體
|
||
（含最寬處，因其落在窗內或 margin 範圍）保留。
|
||
回傳 (vb_mask (z,y,x) bool, ap_thresh, info dict)；
|
||
資料不足或上終板面缺位時 vb_mask = None（info['mode'] 說明原因）。
|
||
info['lat_clip'] 記錄實際施加的側向切點（未施加時為 None）。
|
||
"""
|
||
m = np.asarray(mask_zyx) > 0
|
||
zz, yy, xx = np.nonzero(m)
|
||
info = {'n_bone': int(zz.size), 'n_vb': 0, 'ap_thresh': None, 'mode': 'empty'}
|
||
if zz.size < 50:
|
||
return None, None, info
|
||
if endplate is None:
|
||
info['mode'] = 'no_endplate'
|
||
return None, None, info
|
||
a, b, c, d = endplate['plane']
|
||
X = xx.astype(np.float64)
|
||
Y = yy.astype(np.float64)
|
||
Z = zz.astype(np.float64)
|
||
below = (X * a + Y * b + Z * c - d) <= float(margin)
|
||
u_ap = _ap_axis(sym)
|
||
ap = X * u_ap[0] + Y * u_ap[1] + Z * u_ap[2]
|
||
aps = ap[below]
|
||
if aps.size < 50:
|
||
info['mode'] = 'too_few_below'
|
||
return None, None, info
|
||
th, mode = None, 'quantile'
|
||
if sp_mode == 'gap' and ap_thresh is not None:
|
||
th, mode = float(ap_thresh), 'midline_gap'
|
||
elif sp_mode == 'no_spinous':
|
||
# 棘突已切除：用放寬條件的整體谷底（見 docstring b）
|
||
th, _ = _full_ap_valley(aps, min_side_frac=0.10, max_ratio=0.95)
|
||
if th is not None:
|
||
mode = 'nosp_gap'
|
||
else:
|
||
th = _posterior_min_threshold(aps)
|
||
if th is not None:
|
||
mode = 'nosp_post_min'
|
||
else:
|
||
th, _ = _full_ap_valley(aps)
|
||
if th is not None:
|
||
mode = 'full_gap'
|
||
if th is None:
|
||
th = float(np.quantile(aps, 0.55))
|
||
sel = below & (ap < th)
|
||
if not sel.any():
|
||
info['mode'] = 'no_body_voxels'
|
||
return None, None, info
|
||
# 側向包絡（見 docstring 3）：逐 z（endplate 法線層）在離 AP 切點較遠
|
||
# 的 AP 窗量該層椎體自身側向（鏡稱面 signed distance）寬度，裁掉橫突
|
||
# 前緣（其遠超出該層椎體寬）；層寬隨 z 變化（椎體不同高度寬度不同），
|
||
# 單一 3D 包絡會太寬（被最寬層撐大、中層橫突殘留）。量測不足的 z 用
|
||
# 相鄰 z 的封包插值（np.interp 端點延伸）。
|
||
# sliver_frac：靠切點後側、排除出量測窗的 AP 帶比例（橫突前緣所在）。
|
||
lat_clip = None
|
||
a_s, b_s, c_s, d_s = sym['plane']
|
||
lat = X * a_s + Y * b_s + Z * c_s - d_s
|
||
ap_anter = float(ap[sel].min())
|
||
ext = float(th - ap_anter)
|
||
if ext > 6.0:
|
||
body_sel = sel & (ap <= th - float(sliver_frac) * ext)
|
||
nz = m.shape[0]
|
||
zid = zz.astype(np.int64)
|
||
cnt = np.zeros(nz, dtype=np.int64)
|
||
np.add.at(cnt, zid[body_sel], 1)
|
||
zidx = np.flatnonzero(cnt >= 30)
|
||
if zidx.size > 0:
|
||
zv = zid[body_sel]
|
||
lv = lat[body_sel]
|
||
lo_z = np.full(nz, np.inf)
|
||
hi_z = np.full(nz, -np.inf)
|
||
np.minimum.at(lo_z, zv, lv)
|
||
np.maximum.at(hi_z, zv, lv)
|
||
lo_f = np.interp(np.arange(nz), zidx, lo_z[zidx])
|
||
hi_f = np.interp(np.arange(nz), zidx, hi_z[zidx])
|
||
m_lat = max(float(lat_margin), 0.04 * ext)
|
||
sel = sel & (lat >= lo_f[zid] - m_lat) & (lat <= hi_f[zid] + m_lat)
|
||
lat_clip = (float(lo_f.min()) - m_lat, float(hi_f.max()) + m_lat)
|
||
if not sel.any():
|
||
info['mode'] = 'no_body_voxels'
|
||
return None, None, info
|
||
vb_mask = np.zeros(m.shape, dtype=bool)
|
||
vb_mask[zz[sel], yy[sel], xx[sel]] = True
|
||
info.update(n_vb=int(sel.sum()), ap_thresh=th, mode=mode, lat_clip=lat_clip)
|
||
return vb_mask, th, info
|
||
|
||
def best_upper_endplate_plane(mask_zyx, angle_max=45.0, thresh=4.0,
|
||
n_iter=500, seed=42):
|
||
"""
|
||
3D bone mask (z, y, x) 的最佳「上終板」近似平面 a·x + b·y + c·z = d
|
||
（voxel index 座標；(a,b,c) 為朝上的單位法線）。
|
||
|
||
區域：只用該 level 的前側半邊（y 為 AP 方向，anterior = y 較大側）：
|
||
切點 = bone 的 y 範圍中點（y >= (y_min+y_max)//2），避開後方
|
||
棘突 / 弓根 / 椎管結構，與 2D superior endplate 定義一致。
|
||
|
||
兩階段穩健擬合（舊式單次 3D RANSAC 會被側旁結構 / 上層椎體侵入
|
||
的頂面污染擬出過陡平面，例：L1 擬到 21°）：
|
||
1) Stage A（矢狀線）：每個 (y, x) 欄位取最上方 bone voxel 作為頂面點，
|
||
再收斂成 per-y 最上 z，RANSAC 擬合矢狀 (y-z) 線 z = s·y + i，
|
||
取得穩健的矢狀斜率 s；
|
||
2) Stage B（側向線）：先把頂面點筛到矢狀線 ±max(2·thresh, 5) 帶內
|
||
（移除離帶的高/低污染點），再 RANSAC 擬合 (z − s·y − i) = c·x + k，
|
||
取得側向斜率 c；
|
||
3) 平面 z = s·y + c·x + (i+k) → 單位法線 (−c,−s,1)/√(1+s²+c²)、
|
||
d = (i+k)/√(1+s²+c²)、tilt_deg = 法線與 +z 軸夾角；
|
||
inlier_ratio = 全部前側半頂面點落在平面 ±thresh 的比例。
|
||
回傳 dict（資料不足時回傳 None）：
|
||
plane : (a, b, c, d)
|
||
normal / offset
|
||
tilt_deg : 法線與 +z 軸的夾角（上終板傾斜）
|
||
inlier_ratio : 頂面點落在平面 ±thresh 的比例
|
||
n_points / n_inliers
|
||
u, v : 平面內正交方向（供繪製用）
|
||
"""
|
||
from sklearn.linear_model import RANSACRegressor
|
||
|
||
m = np.asarray(mask_zyx) > 0
|
||
nz, ny, nx = m.shape
|
||
if not m.any():
|
||
return None
|
||
idx = np.where(m, np.arange(nz)[:, None, None], -1)
|
||
ztop = idx.max(axis=0) # (y, x) 每欄最上 z
|
||
# 前側半邊：bone y 範圍中點為切點（前側 = y >= 切點）
|
||
y_present = np.where(np.any(m, axis=(0, 2)))[0]
|
||
y_split = int((y_present[0] + y_present[-1]) // 2)
|
||
# ztop 是 (y, x)：條件作用在 y 軸（axis 0）
|
||
sel = (ztop >= 0) & (np.arange(ny)[:, None] >= y_split)
|
||
yy, xx = np.where(sel)
|
||
if xx.size < 20:
|
||
return None
|
||
xa = xx.astype(np.float64)
|
||
ya = yy.astype(np.float64)
|
||
za = ztop[yy, xx].astype(np.float64)
|
||
n = za.size
|
||
|
||
def _line1d(x1, y1, residual):
|
||
"""RANSAC 擬線 y = k·x + b；失敗時退回普通最小二乘。"""
|
||
try:
|
||
fit = RANSACRegressor(residual_threshold=float(residual),
|
||
max_trials=int(n_iter),
|
||
random_state=int(seed)).fit(x1.reshape(-1, 1), y1)
|
||
return float(fit.estimator_.coef_[0]), float(fit.estimator_.intercept_)
|
||
except Exception:
|
||
k, b = np.polyfit(x1, y1, 1)
|
||
return float(k), float(b)
|
||
|
||
# ---- Stage A：矢狀 (y-z) 線 → 穩健矢狀斜率 s ----
|
||
zt = np.full(ny, -1.0)
|
||
np.maximum.at(zt, yy, za)
|
||
yv = np.where(zt >= 0)[0]
|
||
if yv.size < 10:
|
||
return None
|
||
s, i = _line1d(yv.astype(np.float64), zt[yv], 5.0)
|
||
if abs(float(np.degrees(np.arctan(s)))) > angle_max:
|
||
return None
|
||
|
||
# ---- Stage B：固定 s，擬側向 (x) 線 (z − s·y − i) = c·x + k ----
|
||
# 合起來 z = s·y + c·x + (i + k) → n0·p = d0，n0 = (−c, −s, 1), d0 = i + k
|
||
res_sag = za - s * ya - i
|
||
band = np.abs(res_sag) <= max(2.0 * float(thresh), 5.0)
|
||
if int(band.sum()) < 20:
|
||
band = np.ones(n, dtype=bool)
|
||
c, k = _line1d(xa[band], res_sag[band], float(thresh))
|
||
d0 = i + k
|
||
|
||
n0 = np.array([-c, -s, 1.0])
|
||
norm0 = float(np.linalg.norm(n0))
|
||
nv = n0 / norm0
|
||
d = float(d0 / norm0)
|
||
tilt_deg = float(np.degrees(np.arccos(np.clip(nv[2], -1.0, 1.0))))
|
||
if tilt_deg > angle_max:
|
||
return None
|
||
|
||
dist = np.abs(-c * xa - s * ya + za - d0) / norm0
|
||
n_in = int(np.count_nonzero(dist <= float(thresh)))
|
||
u = np.cross(nv, [1.0, 0.0, 0.0])
|
||
u = u / np.linalg.norm(u)
|
||
v = np.cross(nv, u)
|
||
return {
|
||
'plane': (float(nv[0]), float(nv[1]), float(nv[2]), d),
|
||
'normal': (float(nv[0]), float(nv[1]), float(nv[2])),
|
||
'offset': d,
|
||
'tilt_deg': tilt_deg,
|
||
'inlier_ratio': float(n_in / n),
|
||
'n_points': int(n),
|
||
'n_inliers': n_in,
|
||
'u': (float(u[0]), float(u[1]), float(u[2])),
|
||
'v': (float(v[0]), float(v[1]), float(v[2])),
|
||
}
|
||
|
||
def azimuth_rotation(image, show_plt=False, save_plt=False, output_path=None):
|
||
|
||
img = sitk.ReadImage(image, sitk.sitkUInt8)
|
||
arr_zyx = sitk.GetArrayFromImage(img) # (z, y, x)
|
||
|
||
max_proj = np.max(arr_zyx, axis=0) # -> (y, x)
|
||
|
||
binary_proj = (max_proj > 0).astype(np.uint8)
|
||
|
||
ys, xs = np.where(binary_proj > 0)
|
||
if len(xs) < 10:
|
||
raise ValueError("Not enough foreground points")
|
||
|
||
# 2) centroid ←←← 這裡一定會定義 cx, cy
|
||
cy = ys.mean()
|
||
cx = xs.mean()
|
||
centroid = np.array([cy, cx])
|
||
|
||
cy = ys.mean()
|
||
cx = xs.mean()
|
||
centroid = np.array([cy, cx])
|
||
|
||
y_min = ys.min()
|
||
top_row_mask = (ys == y_min)
|
||
xs_top_row = xs[top_row_mask]
|
||
|
||
# 取這一排的中位數或平均值
|
||
x_center = int(np.median(xs_top_row)) # 或用 np.mean()
|
||
top_point = (y_min, x_center)
|
||
# print(f"最上排中心點: {top_point}")
|
||
|
||
# 計算從 top_point 到 centroid 的向量
|
||
dy = cy - y_min # y 方向的變化
|
||
dx = cx - x_center # x 方向的變化
|
||
|
||
# 計算與 y 軸的夾角
|
||
# 注意：影像座標系中 y 軸向下，所以要特別處理
|
||
angle_rad = np.arctan2(dx, dy) # 弧度
|
||
angle_deg = np.degrees(angle_rad) # 轉成角度
|
||
|
||
if show_plt:
|
||
|
||
fig = plt.figure(figsize=(12, 12))
|
||
|
||
# 視覺化時加上角度資訊
|
||
plt.imshow(binary_proj, cmap='gray')
|
||
|
||
plt.scatter(x_center, y_min, c='red', s=60, label='Top')
|
||
plt.scatter(cx, cy, c='yellow', s=60, label='Centroid')
|
||
plt.plot([x_center, cx], [y_min, cy], 'c-', lw=2,
|
||
label=f'Angle with y-axis: {angle_deg:.1f}°')
|
||
|
||
# 中矢狀面：最佳鏡稱面 a·x + b·y + c·z = d（法線方向任意，不平行 YZ 面）
|
||
# 畫該平面在體積中央 z 切片上的截線
|
||
sym3 = best_symmetry_plane(arr_zyx)
|
||
a3, b3, c3, d3 = sym3['plane']
|
||
zc_ = (arr_zyx.shape[0] - 1) / 2.0
|
||
rhs_ = d3 - c3 * zc_ # a*x + b*y = rhs
|
||
n2d = np.hypot(a3, b3)
|
||
x0_ = a3 * rhs_ / (n2d * n2d)
|
||
y0_ = b3 * rhs_ / (n2d * n2d)
|
||
u2_ = np.array([b3, -a3]) / n2d
|
||
tfg_ = (xs - x0_) * u2_[0] + (ys - y0_) * u2_[1]
|
||
L_ = float(np.abs(tfg_).max())
|
||
plt.plot([x0_ - L_ * u2_[0], x0_ + L_ * u2_[0]],
|
||
[y0_ - L_ * u2_[1], y0_ + L_ * u2_[1]],
|
||
color='magenta', lw=2, linestyle='--',
|
||
label=f'Mirror plane: {a3:+.2f}x {b3:+.2f}y {c3:+.2f}z = {d3:.1f}')
|
||
|
||
plt.legend()
|
||
plt.title("Top point + centroid-directed line")
|
||
plt.axis("off")
|
||
|
||
if save_plt:
|
||
if output_path is None:
|
||
output_path = "azimuth_rotation.png"
|
||
fig.savefig(output_path, dpi=200, bbox_inches="tight")
|
||
|
||
plt.show()
|
||
plt.close(fig)
|
||
|
||
return angle_deg
|
||
|
||
def split_spine_anterior_posterior(image_path, center_mode='com'):
|
||
"""
|
||
從 sagittal view 看，沿著 y 軸（前後方向）將脊椎切成前半部和後半部
|
||
|
||
Parameters:
|
||
image_path (str): 影像路徑
|
||
center_mode (str): 'com' 使用質心的 y 座標，'image' 使用圖片中心的 y 座標
|
||
|
||
Returns:
|
||
anterior, posterior: 前半部和後半部的 binary mask
|
||
"""
|
||
|
||
# Sagittal projection: 沿著 x 軸投影 -> (z, y)
|
||
img = sitk.ReadImage(image_path, sitk.sitkUInt8)
|
||
arr_zyx = sitk.GetArrayFromImage(img)
|
||
|
||
# Sagittal projection
|
||
max_proj_sagittal = np.max(arr_zyx, axis=2)
|
||
binary_proj = (max_proj_sagittal > 0).astype(np.uint8)
|
||
|
||
# 決定切割的 y 座標
|
||
if center_mode == 'com':
|
||
cz, cy = center_of_mass(binary_proj)
|
||
split_y = int(round(cy))
|
||
label = f'Center of Mass (y={split_y})'
|
||
|
||
elif center_mode == 'image':
|
||
split_y = binary_proj.shape[1] // 2
|
||
label = f'Image Center (y={split_y})'
|
||
|
||
# 檢查是否為數字，且範圍在 0 到 1 之間 (不含邊界)
|
||
elif isinstance(center_mode, (int, float)) and 0 < center_mode < 1:
|
||
ys = np.where(binary_proj > 0)[1]
|
||
|
||
if ys.size == 0: # 額外保險：如果投影是空的
|
||
split_y = binary_proj.shape[1] // 2
|
||
else:
|
||
y_min, y_max = ys.min(), ys.max()
|
||
split_y = int(round(y_min + center_mode * (y_max - y_min)))
|
||
|
||
label = f'Custom Ratio {center_mode} (y={split_y})'
|
||
|
||
else:
|
||
raise ValueError("center_mode 必須是 'com'、'image' 或介於 0 到 1 之間的浮點數 (例如 0.8)")
|
||
|
||
# 切割：anterior (y < split_y) 和 posterior (y >= split_y)
|
||
anterior = binary_proj.copy()
|
||
posterior = binary_proj.copy()
|
||
|
||
anterior[:, :split_y] = 0 # 保留spine前半部(image後半部)（y >= split_y）
|
||
posterior[:, split_y:] = 0 # 保留spine後半部(image前半部)（y < split_y）
|
||
|
||
return anterior, posterior, binary_proj
|
||
|
||
def analyze_vertebral_tilt_contour(image_path, edge_type='superior', show_plot=False, debug=False, save_plt=False, output_path=None):
|
||
"""
|
||
通過椎體前緣輪廓分析傾斜（可選上或下終板）
|
||
|
||
Parameters:
|
||
edge_type: 'superior' 上終板, 'inferior' 下終板, 'both' 兩者都分析
|
||
"""
|
||
|
||
# 切割出 anterior 部分
|
||
anterior, posterior, binary_proj = split_spine_anterior_posterior(image_path, center_mode='com')
|
||
|
||
zs, ys = np.where(anterior > 0)
|
||
|
||
if len(zs) == 0:
|
||
return None
|
||
|
||
from sklearn.linear_model import RANSACRegressor
|
||
|
||
results = {}
|
||
|
||
# === 根據 edge_type 決定要分析哪些邊 ===
|
||
edges_to_analyze = []
|
||
if edge_type == 'superior' or edge_type == 'both':
|
||
edges_to_analyze.append('superior')
|
||
if edge_type == 'inferior' or edge_type == 'both':
|
||
edges_to_analyze.append('inferior')
|
||
|
||
all_edge_points = {}
|
||
all_inliers = {}
|
||
all_outliers = {}
|
||
all_slopes = {}
|
||
all_intercepts = {}
|
||
all_angles = {}
|
||
|
||
for edge in edges_to_analyze:
|
||
# 對每個 y，找對應的邊緣點
|
||
edge_points = []
|
||
unique_ys = np.unique(ys)
|
||
|
||
for y in unique_ys:
|
||
z_at_y = zs[ys == y]
|
||
|
||
if edge == 'superior':
|
||
z_edge = z_at_y.max() # 最上面的點（z 最小）
|
||
else: # inferior
|
||
z_edge = z_at_y.min() # 最下面的點（z 最大）
|
||
|
||
edge_points.append([y, z_edge])
|
||
|
||
edge_points = np.array(edge_points)
|
||
all_edge_points[edge] = edge_points
|
||
|
||
if len(edge_points) < 10:
|
||
continue
|
||
|
||
# RANSAC 擬合
|
||
X = edge_points[:, 0].reshape(-1, 1)
|
||
y_data = edge_points[:, 1]
|
||
|
||
ransac = RANSACRegressor(
|
||
residual_threshold=5.0,
|
||
random_state=42
|
||
)
|
||
ransac.fit(X, y_data)
|
||
|
||
inlier_mask = ransac.inlier_mask_
|
||
outlier_mask = ~inlier_mask
|
||
|
||
edge_points_inliers = edge_points[inlier_mask]
|
||
edge_points_outliers = edge_points[outlier_mask]
|
||
|
||
all_inliers[edge] = edge_points_inliers
|
||
all_outliers[edge] = edge_points_outliers
|
||
|
||
# 獲取擬合結果
|
||
slope = ransac.estimator_.coef_[0]
|
||
intercept = ransac.estimator_.intercept_
|
||
|
||
all_slopes[edge] = slope
|
||
all_intercepts[edge] = intercept
|
||
|
||
# 計算 R²
|
||
y_pred = ransac.predict(edge_points_inliers[:, 0].reshape(-1, 1))
|
||
ss_res = np.sum((edge_points_inliers[:, 1] - y_pred) ** 2)
|
||
ss_tot = np.sum((edge_points_inliers[:, 1] - np.mean(edge_points_inliers[:, 1])) ** 2)
|
||
r_squared = 1 - (ss_res / ss_tot) if ss_tot > 0 else 0
|
||
|
||
# 計算傾斜角度
|
||
tilt_angle = np.degrees(np.arctan(slope))
|
||
all_angles[edge] = tilt_angle
|
||
|
||
if debug:
|
||
print(f"\n=== {edge.upper()} ENDPLATE ===")
|
||
print(f"Total points: {len(edge_points)}")
|
||
print(f"Inliers: {len(edge_points_inliers)}")
|
||
print(f"Outliers: {len(edge_points_outliers)}")
|
||
|
||
print(f"{edge.capitalize()} 終板傾斜角度: {tilt_angle:.2f}°")
|
||
print(f"斜率: {slope:.4f}, R²: {r_squared:.4f}")
|
||
|
||
results[edge] = {
|
||
'tilt_angle_deg': tilt_angle,
|
||
'slope': slope,
|
||
'intercept': intercept,
|
||
'r_squared': r_squared,
|
||
'n_inliers': len(edge_points_inliers),
|
||
'n_outliers': len(edge_points_outliers)
|
||
}
|
||
|
||
# Visualization
|
||
if show_plot:
|
||
n_edges = len(edges_to_analyze)
|
||
fig, axes = plt.subplots(n_edges, 2, figsize=(16, 6*n_edges))
|
||
|
||
if n_edges == 1:
|
||
axes = axes.reshape(1, -1)
|
||
|
||
colors = {'superior': 'red', 'inferior': 'cyan'}
|
||
|
||
for idx, edge in enumerate(edges_to_analyze):
|
||
edge_points = all_edge_points[edge]
|
||
edge_points_inliers = all_inliers[edge]
|
||
edge_points_outliers = all_outliers[edge]
|
||
slope = all_slopes[edge]
|
||
intercept = all_intercepts[edge]
|
||
tilt_angle = all_angles[edge]
|
||
color = colors[edge]
|
||
|
||
# 左圖：scatter plot
|
||
ax_left = axes[idx, 0]
|
||
|
||
if len(edge_points_outliers) > 0:
|
||
ax_left.scatter(edge_points_outliers[:, 0], edge_points_outliers[:, 1],
|
||
c='lightcoral', s=30, alpha=0.6, marker='x',
|
||
label=f'Outliers ({len(edge_points_outliers)})', zorder=3)
|
||
|
||
ax_left.scatter(edge_points_inliers[:, 0], edge_points_inliers[:, 1],
|
||
c=color, s=20, alpha=0.7,
|
||
label=f'Inliers ({len(edge_points_inliers)})', zorder=4)
|
||
|
||
# 擬合線
|
||
y_line = np.array([edge_points[:, 0].min(), edge_points[:, 0].max()])
|
||
z_line = slope * y_line + intercept
|
||
ax_left.plot(y_line, z_line, 'lime', linewidth=3,
|
||
label=f'Angle: {tilt_angle:.1f}°\nR²: {results[edge]["r_squared"]:.3f}',
|
||
zorder=5)
|
||
|
||
ax_left.set_xlabel('y')
|
||
ax_left.set_ylabel('z')
|
||
ax_left.set_title(f'{edge.capitalize()} Endplate Analysis\nTilt: {tilt_angle:.1f}°')
|
||
ax_left.invert_yaxis()
|
||
ax_left.legend()
|
||
ax_left.grid(True, alpha=0.3)
|
||
|
||
# 右圖：原始影像
|
||
ax_right = axes[idx, 1]
|
||
ax_right.imshow(anterior, cmap='gray', aspect='equal')
|
||
|
||
if len(edge_points_outliers) > 0:
|
||
ax_right.scatter(edge_points_outliers[:, 0], edge_points_outliers[:, 1],
|
||
c='red', s=40, alpha=0.7, marker='x', label='Outliers', zorder=4)
|
||
|
||
ax_right.scatter(edge_points_inliers[:, 0], edge_points_inliers[:, 1],
|
||
c=color, s=25, alpha=0.8, label='Inliers', zorder=3)
|
||
|
||
ax_right.plot(y_line, z_line, 'lime', linewidth=3, linestyle='--',
|
||
label=f'{edge.capitalize()}: {tilt_angle:.1f}°', zorder=5)
|
||
|
||
ax_right.set_title(f'Anterior Half - {edge.capitalize()} Edge')
|
||
ax_right.set_xlabel('y')
|
||
ax_right.set_ylabel('z')
|
||
ax_right.invert_yaxis()
|
||
ax_right.legend()
|
||
|
||
plt.tight_layout()
|
||
|
||
if save_plt:
|
||
if output_path is None:
|
||
output_path = "analyze_vertebral_tilt_contour.png"
|
||
fig.savefig(output_path, dpi=200, bbox_inches="tight")
|
||
|
||
plt.show()
|
||
plt.close(fig)
|
||
|
||
return results
|