From 2ae08ac2cdcd962b5ad17adea723f3c246196fe0 Mon Sep 17 00:00:00 2001 From: Xiao Furen Date: Mon, 7 Sep 2026 18:46:06 +0800 Subject: [PATCH] refactor(imaging): improve orientation detection and segmentation robustness 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. --- config/constant.py | 2 +- core/optimizer.py | 167 ++++---- core/scoring.py | 18 +- imaging/orientation.py | 198 +++++++++- imaging/preprocessing.py | 75 ++++ utils/helpers.py | 18 +- visualization/res_bone_figure.py | 504 ++++++++++++++++++++--- visualization/res_plot_3d.py | 658 ++----------------------------- xfr_preprocess.py | 20 +- xfr_reprocess_ap.py | 128 ++++++ 10 files changed, 976 insertions(+), 812 deletions(-) create mode 100644 xfr_reprocess_ap.py diff --git a/config/constant.py b/config/constant.py index 9641f2c..8e70484 100644 --- a/config/constant.py +++ b/config/constant.py @@ -24,7 +24,7 @@ ALLOWED_DIAMETERS = [ # 7.5, ] ALLOWED_LENGTHS = [ - 25, + # 25, 30, 35, 40, diff --git a/core/optimizer.py b/core/optimizer.py index 9ffe872..5e30f0d 100644 --- a/core/optimizer.py +++ b/core/optimizer.py @@ -14,7 +14,7 @@ from core.cylinder import generate_cylinder_n_torch, generate_cylinder_o_torch, from core.intersection import center_line_intersections_torch from core.scoring import compute_overlap_ratio_from_cylinder_mask, is_solution_ok, cl_score_torch_xfr from config.constant import OVERLAP_THRESH -from visualization.res_plot_3d import res_plt_2_torch +from visualization.res_bone_figure import render_bone_figure LATERAL_REFINE_MIN_IN_BONE = 0.97 @@ -25,7 +25,7 @@ VBODY_ENTRY_EDGE = 4 def _makedirs_retry(path, retries=5, delay=0.5): - """NFS 上建目錄重試（同 res_plot_3d._retry_robust 處理的瞬時錯誤）""" + """NFS 上建目錄重試（同 utils.helpers.retry_robust 處理的瞬時錯誤）""" for i in range(retries): try: os.makedirs(path, exist_ok=True) @@ -259,7 +259,7 @@ def run_pso_torch_xfr( # 上終板面與椎體都在「完整 mask（SP 移除前）」上計算： # - 終板面由前側頂面擬合，SP 移除不改變平面； - # - 椎體與 res_plt_2_torch 的 gold 顯示完全同 input / 同參數， + # - 椎體與 render_bone_figure 的 gold 顯示完全同 input / 同參數， # 確保顯示出來的椎體就是 loss 裡 VBODY 獎勵的區域。 # 棘突缺如（laminectomy，mode='no_spinous'）時不該把殘留後側要素 # 當「棘突」移除（會鏟進椎體後側），入口面維持完整 mask。 @@ -321,31 +321,32 @@ def run_pso_torch_xfr( # flat_min_index = np.argmin(y_indices) # z_border, x_border = np.unravel_index(flat_min_index, y_indices.shape) - # 脊椎中線：整段 (全體積) 骨頭 x 範圍的中點。 - # 不取單一行的原因 (0005 L5)：椎體軸狀面旋轉時單行只罩到單側骨塊 - # (x 1..76 / W=224 → x_mid≈0.17W)，L/R 兩個 band 被壓到同一側。 - # 不用鏡稱對稱軸的原因 (0001 L4)：逐切面對稱軸會被肋、後側要素 - # 左右不對稱與椎體傾斜牽引 (69.0 vs 範圍中點 74.5)，把 R band 內緣 - # (x_mid+0.1W) 拉進中線棘突/椎板區，R 側入口落在棘突上 (太靠內後)。 - # Laminectomy 只移除中線後側要素，左右極端 x 位置不變， - # 所以範圍中點同樣不受其影響，作為 L/R band 分割線比對稱軸穩定。 - x_with_nonzero = np.where(np.any(image2_array != 0, axis=(0, 1)))[0] - x1 = x_with_nonzero[0] - x2 = x_with_nonzero[-1] + # x/z 搜索空間邊界：把 VBODY 椎體投影到 xz 平面，取該投影的 bounding box， + # 再在 x 向切 L/R 兩 band + 中央缺口、z 向留 10%~90%（見下方 CBT bounds）。 + # CBT 入口點應落在椎體上（左/右 band），而非跨整段骨頭：整段骨頭包含 + # 肋、後側要素、橫突等極端，x 範圍比椎體寬，會把 L/R band 往外推。 + # （不取單一行 / 不用鏡稱對稱軸的原因同前，舊註保留於下）。 + # VBODY mask 為 (z,y,x)；np.any(..., axis=1) 折疊 y 得 xz 投影 (z,x)。 + # x 範圍 = 有 VBODY 的欄（投影 axis=0 是 z，沿 z 做 any） + # z 範圍 = 有 VBODY 的行（投影 axis=1 是 x，沿 x 做 any） + # VBODY 分割失敗（None / 全 0）時退回整段骨頭 x/z 範圍。 + if vb_mask_np is not None and vb_mask_np.any(): + vb_xz = np.any(vb_mask_np, axis=1) # (z, x) 投影 + x1 = int(np.where(vb_xz.any(axis=0))[0][0]) + x2 = int(np.where(vb_xz.any(axis=0))[0][-1]) + z1 = int(np.where(vb_xz.any(axis=1))[0][0]) + z2 = int(np.where(vb_xz.any(axis=1))[0][-1]) + print(f"[BOUNDS] VBODY xz projection: z[{z1},{z2}] x[{x1},{x2}] " + f"(z_height={z2 - z1}, x_width={x2 - x1})") + else: + print("[BOUNDS] VBODY unavailable, falling back to whole-bone xz range") + x_with_nonzero = np.where(np.any(image2_array, axis=(0, 1)))[0] + x1 = int(x_with_nonzero[0]) + x2 = int(x_with_nonzero[-1]) + z1 = int(z_with_nonzero[0]) + z2 = int(z_with_nonzero[-1]) x_width = x2 - x1 - - # print(x1,x2) - # exit() - - # x_mid = (x1 + x2) / 2 - # x1 = x_mid-image_shape[2]*.1 - # x2 = x_mid+image_shape[2]*.1 - - z_sum = np.sum(image2_array, axis=(1, 2)) - z_with_nonzero = np.where(z_sum > 0)[0] - z1 = z_with_nonzero[0] - z2 = z_with_nonzero[-1] - z_height = z2-z1 + z_height = z2 - z1 # print(x1,x2) # exit() @@ -363,14 +364,14 @@ def run_pso_torch_xfr( # z_bounds = (0, image_shape[0]-1) # z_bounds = (z1, (z1+z2)/2) # z_bounds = (.1*image_shape[0], .8*image_shape[0]) - z_bounds = (z1+z_height*.1, z1+z_height*.9) + z_bounds = (0, z1+z_height*.8) # x_bounds_right = (image_shape[2]/2 + image_shape[2]/10, image_shape[2] - 1) # x_bounds_left = (0, image_shape[2]/2 - image_shape[2]/10 - 1) # x_bounds_right = (x2, image_shape[2]*.9) # x_bounds_left = (image_shape[2]*.1, x1) - x_bounds_right = (x1+x_width*.6, +x_width*.9) - x_bounds_left = (x1+x_width*.1, +x_width*.4) + x_bounds_left = (x1+x_width*.1, x1+x_width*.4) + x_bounds_right = (x1+x_width*.6, x1+x_width*.9) # 脊椎若被體積邊界切到（真正偏心、骨頭貼著左/右邊緣）， # 對應那側的 x band 下限會 >= 上限，PSO 會丟 "upper-bound must be greater"。 @@ -392,8 +393,8 @@ def run_pso_torch_xfr( # azimuth_bounds_l = ((98-azi), (120-azi)) # azimuth_bounds_r = ((60-azi), (82-azi)) # altitude_bounds = ((60-alt), (70-alt)) - azimuth_bounds_l = (98, 110) - azimuth_bounds_r = (70, 82) + azimuth_bounds_l = (98, 105) + azimuth_bounds_r = (75, 82) altitude_bounds_l = (60, 65) altitude_bounds_r = (60, 65) @@ -571,26 +572,20 @@ def run_pso_torch_xfr( final_length_r = length def _plot_combined(d_l, l_l, d_r, l_r, pos_l, pos_r, total_t): - res_plt_2_torch( + # volume / level 由 image2_path 反推（…/{vol}/rotated/{level}_*.nii.gz） + render_bone_figure( + None, None, spine_tensor, cortical_tensor, - image_shape, - image2_path, folder, - label_str, - d_l, - l_l, - d_r, - l_r, - pos_l, - pos_r, - swarm_size, - max_iter, - total_t, - spacing, - CBT, - device, - grid, + spacing=spacing, + way='CBT' if CBT else 'TPS', + best_position_l=pos_l, best_position_r=pos_r, + diameter_l=d_l, length_l=l_l, + diameter_r=d_r, length_r=l_r, + image2_path=image2_path, + device=device, grid=grid, + swarm_size=swarm_size, max_iter=max_iter, total_time=total_t, ) if side == 'both': @@ -606,7 +601,7 @@ def run_pso_torch_xfr( # 兩側可能在不同 GPU worker：各自把結果寫 _.json # （先寫 tmp 再 os.replace，對端讀到的一定是完整檔）。先完成者看不到 # 對端檔就跳過；後完成者看到兩側齊了、搶到 plot lock（O_EXCL， - # 確保合併輸出只跑一次）才載入對端結果跑 res_plt_2_torch + # 確保合併輸出只跑一次）才載入對端結果跑 render_bone_figure # （3D 圖 + CSV 兩行，與 'both' 模式相同）。json / lock 保留供事後 # 檢查；若該側流程死在 plotting 中段，該 (volume, level) 重跑即可 # （run_id 是新的一次，不會互相干擾）。 @@ -956,26 +951,20 @@ def run_pso_torch( final_diameter_r = diameter final_length_r = length - res_plt_2_torch( - spine_tensor, - cortical_tensor, - image_shape, - image2_path, - folder, - label_str, - final_diameter_l, - final_length_l, - final_diameter_r, - final_length_r, - best_position_l, - best_position_r, - swarm_size, - max_iter, - total_time, - spacing, - CBT, - device, - grid) + render_bone_figure( + None, None, + spine_tensor, + cortical_tensor, + folder, + spacing=spacing, + way='CBT' if CBT else 'TPS', + best_position_l=best_position_l, best_position_r=best_position_r, + diameter_l=final_diameter_l, length_l=final_length_l, + diameter_r=final_diameter_r, length_r=final_length_r, + image2_path=image2_path, + device=device, grid=grid, + swarm_size=swarm_size, max_iter=max_iter, total_time=total_time, + ) return best_position_l, best_loss_l, best_position_r, best_loss_r, total_time @@ -990,7 +979,7 @@ from config.constant import ALLOWED_DIAMETERS, ALLOWED_LENGTHS from core.cylinder import generate_cylinder_n_torch, snap_to_discrete_values, create_coordinate_grid from core.scoring import compute_overlap_ratio_from_cylinder_mask, is_solution_ok from config.constant import OVERLAP_THRESH -from visualization.res_plot_3d import res_plt_2_torch +from visualization.res_bone_figure import render_bone_figure def run_de_torch( label_str: str, @@ -1154,12 +1143,20 @@ def run_de_torch( final_diameter_r = best_position_r[5] if optimize_size else diameter final_length_r = best_position_r[6] if optimize_size else length - res_plt_2_torch( - spine_tensor, cortical_tensor, image_shape, image2_path, 'Output', label_str, - final_diameter_l, final_length_l, final_diameter_r, final_length_r, - best_position_l, best_position_r, swarm_size, max_iter, total_time, spacing, CBT, device, grid + render_bone_figure( + None, None, + spine_tensor, cortical_tensor, + 'Output', + spacing=spacing, + way='CBT' if CBT else 'TPS', + best_position_l=best_position_l, best_position_r=best_position_r, + diameter_l=final_diameter_l, length_l=final_length_l, + diameter_r=final_diameter_r, length_r=final_length_r, + image2_path=image2_path, + device=device, grid=grid, + swarm_size=swarm_size, max_iter=max_iter, total_time=total_time, ) - + return best_position_l, best_loss_l, best_position_r, best_loss_r, total_time def run_nm_torch( @@ -1324,10 +1321,18 @@ def run_nm_torch( final_diameter_r = best_position_r[5] if optimize_size else diameter final_length_r = best_position_r[6] if optimize_size else length - res_plt_2_torch( - spine_tensor, cortical_tensor, image_shape, image2_path, 'Output', label_str, - final_diameter_l, final_length_l, final_diameter_r, final_length_r, - best_position_l, best_position_r, swarm_size, max_iter, total_time, spacing, CBT, device, grid + render_bone_figure( + None, None, + spine_tensor, cortical_tensor, + 'Output', + spacing=spacing, + way='CBT' if CBT else 'TPS', + best_position_l=best_position_l, best_position_r=best_position_r, + diameter_l=final_diameter_l, length_l=final_length_l, + diameter_r=final_diameter_r, length_r=final_length_r, + image2_path=image2_path, + device=device, grid=grid, + swarm_size=swarm_size, max_iter=max_iter, total_time=total_time, ) - + return best_position_l, best_loss_l, best_position_r, best_loss_r, total_time \ No newline at end of file diff --git a/core/scoring.py b/core/scoring.py index 06bd593..cd47c62 100644 --- a/core/scoring.py +++ b/core/scoring.py @@ -17,7 +17,7 @@ def cl_score_torch_xfr( 漸進式評分：優先確保找到骨頭，再改善細節 """ cyl_total = cylinder_torch.sum().item() - overlap = ((cortical_tensor == 1) & (cylinder_torch == 1)).sum().item() + overlap = ((cortical_tensor == 1) & (cylinder_torch == 1)).sum().item() # in cortical # VBODY 獎勵：螺絲落在 (cortical + VBODY) 內的 voxel 每個 100 分 # （100 分項由純 cortical 擴展到 cortical∪VBODY；cortical voxel 分數不變）， # 其中落在 VBODY 的 voxel 每個再加 10 分 @@ -46,14 +46,14 @@ def cl_score_torch_xfr( # if in_bone == 0: # return float(not_in_bone*200) - score += 20 * in_bone # 10 實在太低 - score += 100 * overlap - score += 100 * in_corti_vb # (cortical + VBODY) 每 voxel 100 分 - score += 50 * in_vbody # VBODY 每 voxel 再加 10 分 - # score -= 2000 * max(0, not_in_bone-10) - # score -= 1000 * max(0, null_vox2-10) - score -= 1000 * not_in_bone - score -= 2000 * null_vox2 + score += 10 * in_bone # 10 實在太低 + score += 100 * overlap # in cortical + score += 100 * in_corti_vb # (cortical + VBODY) 每 voxel 再加分 + score += 50 * in_vbody # VBODY 每 voxel 再加分 + # score -= 1000 * not_in_bone + score -= 1000 * max(0, not_in_bone-10) + # score -= 2000 * null_vox2 + score -= 2000 * max(0, null_vox2-10) return float(-score) diff --git a/imaging/orientation.py b/imaging/orientation.py index 11ad440..b5bcf4c 100644 --- a/imaging/orientation.py +++ b/imaging/orientation.py @@ -249,9 +249,12 @@ def segment_spinous_process(mask_zyx, sym, band_frac=0.06, min_band=6.0, w = max(min_band, band_frac * s 全寬) 2) 前後方向：平面內兩軸 (u, v) 中 |y| 分量大者， 正規化成 +AP = 後側（本資料系 y 往前遞增，後側 = y 小側） - 3) 中線帶的 AP 分佈呈兩大叢（椎體在前、椎弓/棘突在後）， - 以兩叢間的 AP 谷底為界，AP >= 谷底 的中線帶 voxel = 棘突（含中線椎弓）； - 無明顯谷底（如骨橋）fallback 取中線帶後側 15%。 +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。 @@ -301,39 +304,58 @@ def segment_spinous_process(mask_zyx, sym, band_frac=0.06, min_band=6.0, 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() - best_i, best_score = None, -1.0 - for i in range(len(hist)): - if hist[i] >= 0.05 * peak: + # 收集連續安靜 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 - if csum[i] < min_mass_frac * total or (total - csum[i + 1]) < min_mass_frac * total: - continue - score = min(csum[i], total - csum[i + 1]) - if score > best_score: - best_score, best_i = score, i - if best_i is not None: - post_frac = (total - csum[best_i + 1]) / total + 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% 中線帶質量）＝棘突幾乎全除 + # 空隙後側叢只剩小殘片（<5% 中線帶質量）＝棘突幾乎全除 # （部分切除殘餘）：當作缺如，椎體切分不採用此閾值。 # 僅在鏡稱面左右為主時才算數（斜板層時 post_frac 不可信） info['mode'] = 'no_spinous' info['post_frac'] = float(post_frac) return None, None, info - th = float(0.5 * (edges[best_i] + edges[best_i + 1])) + 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) - sel = mid & (ap >= th) + # 棘突用空隙後側端 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, mode=mode) + info.update(n_sp=int(sp_mask.sum()), ap_thresh=th, ap_thresh_sp=th_sp, mode=mode) return sp_mask, th, info @@ -429,6 +451,98 @@ def _ap_axis(sym): 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 的質量（拒絕對小尾巴的偽谷底）。 @@ -483,7 +597,8 @@ def _posterior_min_threshold(aps, rear_frac=0.40, min_side_frac=0.08, smooth=3): 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): +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）： @@ -501,10 +616,18 @@ def segment_vertebral_body(mask_zyx, sym, endplate, ap_thresh, sp_mode, margin=2 （_posterior_min_threshold）。 c) 否則：終板下整體 AP 分佈的平滑谷底（_full_ap_valley， 處理中線骨橋等中線搜尋 fallback 的情形）； - d) 最後 fallback：AP 分佈 55 百分位（可能切進椎體內，會打 WARNING）。 - 椎體 = AP < 切點（切點之前側）且終板下側的 bone。 +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) @@ -544,12 +667,45 @@ def segment_vertebral_body(mask_zyx, sym, endplate, ap_thresh, sp_mode, margin=2 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) + 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, diff --git a/imaging/preprocessing.py b/imaging/preprocessing.py index 2bc297d..97a4499 100644 --- a/imaging/preprocessing.py +++ b/imaging/preprocessing.py @@ -1,4 +1,5 @@ import os +import numpy as np import SimpleITK as sitk from imaging.resample import resample_img from imaging.affine import standardize_affine @@ -7,6 +8,18 @@ import json import glob from config.constant import LABEL_MAP from imaging.nifti_io import sitk_to_nibabel, nibabel_to_sitk +from imaging.orientation import anterior_y_side + + +def flip_y_sitk(img): + """index 系 y 軸（array axis 1，前後方向）翻轉： + 只翻數據、spacing/origin/direction 等幾何不變，即整顆體積的前後 + 朝向在 index 系翻轉（y 大側 <-> y 小側）。供 supine / prone 個案 + 統一前後慣例（前側 = y 大側）用。""" + arr = np.flip(sitk.GetArrayFromImage(img), axis=1) + out = sitk.GetImageFromArray(arr) + out.CopyInformation(img) + return out @@ -111,6 +124,68 @@ def process_single_image(image_path, label_path, output_dir_base=None, max_z_spa resampled_sitk_img = resample_img(image, out_spacing=[0.5, 0.5, 0.5], is_label=False) resampled_sitk_lbl = resample_img(label, out_spacing=[0.5, 0.5, 0.5], is_label=True) + # 前後（AP）方向判定：本流程最終輸出慣例是前側 = y 大側（後 = y 小側）。 + # 但 prone（伏位）掃描個案經 standardize_affine 後前側落在 y 小側 + # （例：CTSpine1K colon 0003、0075、0460...，全 dataset 約 1%）， + # 不修正時上終板 / 棘突 / 椎體分割與 rotated/ 對齊全部反掉。 + # + # 判定（對每個 allowed level 的 0.5mm label 個別做、再票決——不能 + # 直接 union 多層：腰椎前凸（lordosis）下各層椎體在 AP 投影會散開， + # 椎管「空隙」被其它層的骨填掉）： +# 1) anterior_y_side 回傳工作體積（0.5mm 重取樣、standardize_affine + # 之前）grid 中椎體塊所在的 y 端（y_min / y_max / None）； + # 2) standardize_affine（nibabel 端）會在輸出 affine y 分量 < 0 時 + # 再翻一次 y。注意 nibabel affine 與 SimpleITK direction 的 y 分 + # 量符號相反（NIfTI RAS <-> SITK LPS），所以 + # standardize_affine 會翻 y <=> direction[4] > 0； + # 最終 y 端 = pre_side（若會翻 y 則 y_min<->y_max 互換）； + # 3) 最終前側會落 y 小側時，現在先對 CT 與 label 做 y 翻轉 + # （純 index 翻轉、幾何不變），與 standardize_affine 的翻轉 + # 組成淨效果，使所有輸出與其它個案同方向。 + # 無法判定 / 投票平手時維持原方向（寧可不翻、不誤翻）。 + # 判定結果存入 metadata db（ap_flip），重跑免重算。 + ap_flip = (meta or {}).get("ap_flip") + if ap_flip is None: + # standardize_affine 是否會翻轉 y（nibabel affine[1,1] < 0， + # 等价於 sitk direction[4] > 0，兩者符號相反） + std_flips_y = resampled_sitk_img.GetDirection()[4] > 0 + arr_lbl = sitk.GetArrayFromImage(resampled_sitk_lbl) + votes = [] + for n in allowed_label_list: + side = anterior_y_side(arr_lbl == n) + if side is not None: + votes.append(side) + n_min = votes.count("y_min") + n_max = votes.count("y_max") + if n_min > n_max: + pre_side = "y_min" + elif n_max > n_min: + pre_side = "y_max" + else: + pre_side = None + if pre_side is None: + ap_flip = False + print(f"AP orientation undetermined for {name} (votes={votes}); " + f"proceeding with default orientation (anterior = large y)") + else: + final_side = (pre_side if not std_flips_y + else ("y_min" if pre_side == "y_max" else "y_max")) + ap_flip = final_side == "y_min" + print(f"AP orientation for {name}: pre={pre_side} " + f"(std_flips_y={std_flips_y}) -> final={final_side} " + f"[votes={votes}], flip={ap_flip}") + if metadata_cache is not None: + metadata_cache.put(name, {"ap_flip": bool(ap_flip)}) + if ap_flip: + # label（原解析度 raw label）也要翻：seg_bone 的主遮罩鏈 + #（_binary / SMD / _binary_sdf）是用 original_label（= label） + # 算的，不是用 0.5mm resampled label；漏翻時翻轉不生效。 + label = flip_y_sitk(label) + resampled_sitk_img = flip_y_sitk(resampled_sitk_img) + resampled_sitk_lbl = flip_y_sitk(resampled_sitk_lbl) + print(f"AP orientation corrected for {name}: CT and label flipped " + f"along y; outputs unified to anterior = large y") + # 建立每個檔案的輸出資料夾 file_name = os.path.basename(image_path) name = file_name.replace(".nii.gz", "") diff --git a/utils/helpers.py b/utils/helpers.py index bd7383d..b89295d 100644 --- a/utils/helpers.py +++ b/utils/helpers.py @@ -1,7 +1,23 @@ -import os +import errno +import os +import time + import numpy as np import matplotlib.pyplot as plt + +def retry_robust(fn, *args, retries=20, delay=0.5, **kwargs): + """對 ENOENT/EEXIST 重試：NFS 上輸出樹被外部刪除（或多 worker 併發建同一 + output 目錄）會有短暫的 ENOENT 窗口，重試可恢復；其他錯誤直接丟出。""" + for i in range(retries): + try: + return fn(*args, **kwargs) + except OSError as e: + if e.errno not in (errno.ENOENT, errno.EEXIST) or i == retries - 1: + raise + time.sleep(delay) + + def get_unique_filepath(path: str) -> str: """ 如果檔案已存在，自動加 _1, _2... 避免覆蓋 diff --git a/visualization/res_bone_figure.py b/visualization/res_bone_figure.py index 9914f54..fe490a6 100644 --- a/visualization/res_bone_figure.py +++ b/visualization/res_bone_figure.py @@ -1,3 +1,4 @@ +import csv import os from datetime import datetime @@ -11,9 +12,10 @@ import numpy as np import SimpleITK as sitk from scipy.ndimage import map_coordinates -from imaging.orientation import (best_symmetry_plane, best_upper_endplate_plane, - segment_spinous_process, segment_vertebral_body) -from utils.helpers import get_unique_filepath +from imaging.orientation import (azimuth_rotation, analyze_vertebral_tilt_contour, + best_symmetry_plane, best_upper_endplate_plane, + segment_spinous_process, segment_vertebral_body) +from utils.helpers import get_unique_filepath, retry_robust, save_with_unique_name # 體積吸收渲染（Beer-Lambert），與 res_plot_3d 相同： @@ -235,38 +237,105 @@ def _rotate_plane_params(plane, R, c_xyz): return out +def _mask_to_array(m): + """骨頭 / 皮質遮罩：path (nifti) / (z,y,x) ndarray / torch tensor + -> (z,y,x) bool；None / 缺檔 / 空 -> None。""" + if m is None: + return None + if isinstance(m, str): + return _load_mask(m) + arr = m.cpu().numpy() if hasattr(m, 'cpu') else m + arr = np.asarray(arr) + if arr.size == 0: + return None + b = arr > 0 + return b if int(b.sum()) > 0 else None + + def render_bone_figure(volume_id, level, binary_path, cortical_path, - base_folder="/mnt/1248/open2/cyrou/Output", - spacing=(0.5, 0.5, 0.5), way="CBT", - planes_only=False, output_path=None, rotation=None): - """繪製單一 (volume, level) 骨頭 X-ray 圖（不畫螺絲）。 + base_folder="/mnt/1248/open2/cyrou/Output", + spacing=(0.5, 0.5, 0.5), way="CBT", + planes_only=False, output_path=None, rotation=None, + best_position_l=None, best_position_r=None, + diameter_l=None, length_l=None, + diameter_r=None, length_r=None, + image2_path=None, device=None, grid=None, + swarm_size=None, max_iter=None, total_time=None, + write_csv=True): + """統一骨頭 X-ray 四視角圖（合併原 render_bone_figure + res_plt_2_torch）。 - 四視角：預設/axial/coronal/sagittal。 - 內容：皮質 vs 鬆質吸收骨、椎體(gold)、棘突(purple)、 - 中矢狀鏡稱面(orange)、上終板面(green)。無圓柱/中心線。 + 四視角：預設 / axial（俯視 XY）/ 冠狀（後視）/ 矢狀。 + 內容：皮質 vs 鬆質吸收骨、椎體(gold)、棘突(purple，獨立層覆蓋椎體)、 + 中矢狀鏡稱面(orange)、上終板面(green)；螺絲模式另畫中心線(紅) + + 圓柱（L darkcyan / R blue，o 層粉）。 + 繪製採固定分層（不依深度排序）：基底骨 < VBODY < 棘突 < 終板 < 鏡稱面 < 螺絲。 - planes_only=True：只畫骨頭 + 中矢狀鏡稱面 + 上終板面， - 不做棘突 / 椎體（VBODY）分割。 - output_path：若給定，直接存到該路徑（含自動加 _1/_2 防覆蓋）； - 否則存到 base_folder/{date}/{volume_id}/{volume_id} {level}_{way}.png。 - rotation：(R, c_xyz)。給定時，把骨頭點雲與兩個平面都依 R 旋轉（繞 c_xyz）， - 用於畫「對齊後（rotated）」的平面圖；R 作用於 (x,y,z) 向量。 - cortical_path：皮質遮罩路徑，或 (z,y,x) 0/1 陣列（須與 binary_path 同 grid）； - None / 缺檔時全部視為鬆質骨。 + 骨骼輸入： + binary_path 骨頭遮罩：path (nifti) 或 (z,y,x) ndarray / torch tensor + cortical_path 皮質遮罩：同型（path / ndarray / tensor）；None / 缺檔 + 時全部視為鬆質骨 + volume_id/level 可為 None：由路徑反推（…/{vol}/rotated/{level}_*.nii.gz， + rotated 的上一層 = vol） + image2_path 未給定時 = binary_path（path 情形）；TPS 模式用其算 2D + 參考 az/alt（Azimuth/Altitude 用相對角） + + 平面 / 分割： + planes_only=True 只畫平面，不做棘突 / 椎體（VBODY）分割 + rotation=(R, c_xyz) 點雲與平面同依 R 旋轉（繞 c_xyz），用於「對齊後 + （rotated）」的平面圖；R 作用於 (x,y,z) 向量 + （螺絲模式一律做分割：R 側 loss 有 VBODY 獎勵、gold + 顯示需椎體，planes_only 被忽略） + + 螺絲（可選，給定 best_position_l/r 時啟用，optimizer 的輸出）： + best_position = (z, y, x, az, alt[, d, L])，須搭配該側 diameter/length； + device/grid 供圓柱生成；swarm_size/max_iter/total_time 供圖面註記； + write_csv 時在 base_folder/{date}/{vol}/output.csv 依欄位名稱 append + L/R 兩行（舊 schema 自動重映射）。 + + 輸出： + output_path 給定 -> 存該路徑（含自動加 _1/_2 防覆蓋）；None 時 + 螺絲模式: base_folder/{date}/{vol}/{level}_{way}_L{d}_{l}_R{d}_{l}_{swarm}_{iter}.png + 否則: base_folder/{date}/{vol}/{vol} {level}_{way}.png 回傳存檔路徑；無有效骨頭遮罩時回傳 None。 """ - spine = _load_mask(binary_path) + # ---- 路徑反推：volume_id / level / image2_path ---- + if image2_path is None and isinstance(binary_path, str): + image2_path = binary_path + if (volume_id is None or level is None) and image2_path: + parent = os.path.dirname(image2_path) + vol_p = (os.path.basename(os.path.dirname(parent)) + if os.path.basename(parent) == 'rotated' + else os.path.basename(parent)) + volume_id = volume_id or vol_p + level = level or os.path.basename(image2_path).split('_')[0] + if volume_id is None or level is None: + raise ValueError('volume_id / level unknown: 需提供 image2_path，' + '或顯式傳 volume_id 與 level') + + # ---- 螺絲：決定畫哪些側 ---- + side_cfg = {'L': (best_position_l, diameter_l, length_l, 'Left'), + 'R': (best_position_r, diameter_r, length_r, 'Right')} + screw_sides = [] + for s in ('L', 'R'): + pos, d, L, _cn = side_cfg[s] + if pos is None: + continue + if d is None or L is None: + raise ValueError(f'{s} 側：best_position 須搭配 diameter / length') + screw_sides.append(s) + screw_mode = bool(screw_sides) + CBT = str(way).upper() == 'CBT' + + # ---- 骨骼 / 皮質 ---- + spine = _mask_to_array(binary_path) if spine is None: print(f"[skip] {volume_id} {level}: 無/空骨頭遮罩 {binary_path}") return None - - if isinstance(cortical_path, np.ndarray): - cortical = cortical_path > 0 - else: - cortical = _load_mask(cortical_path) + cortical = _mask_to_array(cortical_path) if cortical is None: cortical = np.zeros_like(spine) + image_shape = spine.shape voxel_mm = float(spacing[0]) alpha_cortical = 1.0 - np.exp(-BONE_MU_CORTICAL * voxel_mm) alpha_trabecular = 1.0 - np.exp(-BONE_MU_TRABECULAR * voxel_mm) @@ -279,26 +348,57 @@ def render_bone_figure(volume_id, level, binary_path, cortical_path, sym = best_symmetry_plane(spine) symp = best_upper_endplate_plane(spine) - if planes_only: + if planes_only and not screw_mode: # 只做方向平面，不做棘突 / 椎體分割 + vb_mask = None sp_corti = sp_trab = None vb_corti = vb_trab = None sp_info = {} vb_info = {} else: + # 棘突：鏡稱面中線帶（|s|<=w）且在 AP 谷底之後側；棘突缺如 + #（先前 laminectomy / 棘突切除）時 sp_mask=None，不標示 sp_mask, sp_th, sp_info = segment_spinous_process(spine, sym) - if sp_mask is not None and sp_mask.any(): + if sp_info.get('mode') == 'no_spinous': + if screw_mode: + top_off = f"{sp_info['top_off']:.1f}" \ + if sp_info.get('top_off') is not None else 'n/a' + print(f"[NO-SP] 中線後側缺如（先前 laminectomy / 棘突切除）: " + f"deficit={sp_info.get('deficit', float('nan')):.1f} voxel " + f"({sp_info.get('deficit', 0.0) * 0.5:.1f} mm), " + f"rear3={sp_info.get('rear3')} voxel, top_off={top_off} " + f"-> 不標示棘突；椎體用放寬後側谷底切分") + sp_corti = sp_trab = None + elif sp_mask is not None and sp_mask.any(): sp_corti = sp_mask[z_corti, y_corti, x_corti] sp_trab = sp_mask[z_trab, y_trab, x_trab] + if screw_mode: + sp_n_bone = max(int(spine.sum()), 1) + print(f"[SPINOUS] n={sp_info['n_sp']} " + f"({100.0 * sp_info['n_sp'] / sp_n_bone:.1f}% of bone) " + f"band=+/-{sp_info['band_w']:.1f} voxel " + f"AP>={sp_info['ap_thresh']:.1f} mode={sp_info['mode']}") else: sp_corti = sp_trab = None - vb_mask, vb_th, vb_info = segment_vertebral_body(spine, sym, symp, sp_th, sp_info["mode"]) + # 椎體：上終板之下 + 中線 AP 谷底之前側（對齊基準系下與 loss 的 + # VBODY 獎勵完全同 input / 同參數，顯示的椎體＝計分的椎體） + vb_mask, vb_th, vb_info = segment_vertebral_body(spine, sym, symp, + sp_th, sp_info["mode"]) if vb_mask is not None and vb_mask.any(): vb_corti = vb_mask[z_corti, y_corti, x_corti] vb_trab = vb_mask[z_trab, y_trab, x_trab] + if screw_mode: + print(f"[VBODY] n={vb_info['n_vb']} " + f"({100.0 * vb_info['n_vb'] / max(int(spine.sum()), 1):.1f}% of bone) " + f"AP<{vb_info['ap_thresh']:.1f} mode={vb_info['mode']}") + if vb_info['mode'] == 'quantile': + print(f"[VBODY] WARNING: 未找到體/弓後側谷底，閾值退回 55 百分位 " + f"（可能切進椎體內），建議人工核對該 level 的椎體邊界") else: vb_corti = vb_trab = None + if screw_mode: + print(f"[VBODY] skipped: {vb_info['mode']}") # ---- 骨頭點雲（體積吸收）---- x_bone = np.concatenate([x_corti, x_trab]) @@ -317,16 +417,23 @@ def render_bone_figure(volume_id, level, binary_path, cortical_path, bone_rgba = bone_rgba[::BONE_SUBSAMPLE] bone_size = bone_size[::BONE_SUBSAMPLE] + # VBODY / 棘突的 voxel flags（corti+trab 接合陣列上）；兩 mask 重合時 + # 歸棘突（與 label map 2 覆蓋 1 一致） + vb_flag = None + sp_flag = None if vb_corti is not None: - vb_flag = np.concatenate([vb_corti, vb_trab]) + vb_flag = np.concatenate([vb_corti, vb_trab]).astype(bool) if BONE_SUBSAMPLE > 1: vb_flag = vb_flag[::BONE_SUBSAMPLE] - bone_rgba[vb_flag] = to_rgba("gold", 0.95) if sp_corti is not None: - sp_flag = np.concatenate([sp_corti, sp_trab]) + sp_flag = np.concatenate([sp_corti, sp_trab]).astype(bool) if BONE_SUBSAMPLE > 1: sp_flag = sp_flag[::BONE_SUBSAMPLE] - bone_rgba[sp_flag] = to_rgba("purple", 0.95) + if vb_flag is None: + vb_flag = np.zeros(len(x_bone), dtype=bool) + if sp_flag is None: + sp_flag = np.zeros(len(x_bone), dtype=bool) + vb_flag &= ~sp_flag # ---- 旋轉對齊（若給定）：骨頭點雲與平面同旋轉 ---- if rotation is not None: @@ -340,6 +447,146 @@ def render_bone_figure(volume_id, level, binary_path, cortical_path, if symp is not None: symp = _rotate_plane_params(symp, R, c_xyz) + # ---- 拆三層：基底骨 / VBODY (gold) / 棘突 (purple) ---- + # x_bone 保留完整點雲（含 VBODY / SP）供下方平面 patch 算範圍 + #（VBODY 前側是整顆骨最前緣，剔除後綠色終板 patch 會縮小）； + # mpl 3D scatter 在同一 collection 內依深度排序 markers， + # 「棘突覆蓋椎體、螺絲覆蓋骨頭」改以固定 zorder 分層達成（見 _fill_ax） + base_idx = ~(vb_flag | sp_flag) + x_base, y_base, z_base = x_bone[base_idx], y_bone[base_idx], z_bone[base_idx] + rgba_base, size_base = bone_rgba[base_idx], bone_size[base_idx] + vb_idx = vb_flag & ~sp_flag + x_vb, y_vb, z_vb = x_bone[vb_idx], y_bone[vb_idx], z_bone[vb_idx] + x_sp, y_sp, z_sp = x_bone[sp_flag], y_bone[sp_flag], z_bone[sp_flag] + + # ---- 螺絲：圓柱 + 中心線 + loss（lazy import torch / core.*； + # 無螺絲的 preprocess 路徑不會載入 torch）---- + side_info = {} + side_azlat, side_acep = {}, {} + theta_v = tau_y = tau_x = float('nan') + azi = alt = float('nan') + x_screw = y_screw = z_screw = None + screw_rgba = screw_size = None + if screw_mode: + if device is None: + raise ValueError('螺絲模式需要 device（torch device）') + spacing = list(spacing) # core.cylinder 以 list 比對 spacing + import torch + from core.cylinder import (generate_cylinder_n_torch, + generate_cylinder_o_torch) + from core.intersection import center_line_intersections_torch + from core.scoring import cl_score_torch, cl_score_torch_xfr + + # TPS 用 2D 參考角（CBT 無 2D 參考，az/alt 為 nan） + if not CBT: + if image2_path: + azi = float(azimuth_rotation(image2_path)) + alt = float(analyze_vertebral_tilt_contour( + image2_path, edge_type='superior', + show_plot=False, debug=False)['superior']['tilt_angle_deg']) + else: + print('[warn] TPS 模式缺 image2_path；' + 'Azimuth/Altitude 退回原始角') + + spine_tensor = torch.from_numpy(spine.astype(np.uint8)).to(device) + cortical_tensor = torch.from_numpy(cortical.astype(np.uint8)).to(device) + vbody_tensor = None + if vb_mask is not None and vb_mask.any(): + # R 側 loss 使用與 PSO 目標函數相同的 VBODY 獎勵 + #（回報分數與優化一致） + vbody_tensor = torch.from_numpy(vb_mask.astype(np.uint8)).to(device=device) + + # 角度註記（與 CSV 同參數）： + # Azimuth_Lateral = 螺絲在鏡稱面內相對 AP 軸的發散角（+ = L 側往外，− = R 側） + # Altitude_Endplate = 螺絲相對上終板面的仰角 + sym_n = np.asarray(sym['normal'], dtype=float) + theta_v = float(np.degrees(np.arctan2(sym_n[1], sym_n[0]))) + if symp is not None: + e_n = np.asarray(symp['normal'], dtype=float) + e_n = e_n / np.linalg.norm(e_n) + tau_y = float(np.degrees(np.arctan2(e_n[1], e_n[2]))) + tau_x = float(np.degrees(np.arctan2(e_n[0], e_n[2]))) + else: + e_n = None + + def _rel_angles(az_deg, alt_deg): + az_r = np.radians(az_deg) + alt_r = np.radians(alt_deg) + d_v = np.array([np.cos(az_r) * np.sin(alt_r), + np.sin(az_r) * np.sin(alt_r), + np.cos(alt_r)]) + az_lateral = az_deg - theta_v - 90.0 + alt_cep = (90.0 - float(np.degrees(np.arccos( + np.clip(d_v @ e_n, -1.0, 1.0)))) + if e_n is not None else float('nan')) + return az_lateral, alt_cep + + for s in screw_sides: + pos, d, L, _cn = side_cfg[s] + # 螺絲方向向量（與 generate_cylinder_n_torch 同慣例）： + # d = (cos(az)·sin(alt), sin(az)·sin(alt), cos(alt))，alt = 相對 +z 的極角 + cyl_n = generate_cylinder_n_torch(d, L, pos[0], pos[1], pos[2], + pos[3], pos[4], + image_shape, spacing, device, grid) + cyl_o = generate_cylinder_o_torch(d, L, pos[0], pos[1], pos[2], + pos[3], pos[4], + image_shape, spacing, device, grid) + inter, line_mask = center_line_intersections_torch( + pos[0], pos[1], pos[2], pos[3], pos[4], + int(L), spine_tensor, spacing, device) + if s == 'L': + loss = cl_score_torch(cortical_tensor, spine_tensor, + cyl_n, cyl_o, inter) + else: + loss = cl_score_torch_xfr(cortical_tensor, spine_tensor, + cyl_n, cyl_o, inter, + vbody_tensor=vbody_tensor) + cyl_points = int(torch.sum(cyl_n).item()) + ovc = (100.0 * int(((cortical_tensor == 1) & (cyl_n == 1)).sum().item()) + / cyl_points) if cyl_points else 0.0 + ovb = (100.0 * int(((spine_tensor == 1) & (cyl_n == 1)).sum().item()) + / cyl_points) if cyl_points else 0.0 + side_info[s] = { + 'line': np.where(line_mask.cpu().numpy() == 1), + 'cyl_n': np.where(cyl_n.cpu().numpy() == 1), + 'cyl_o': np.where(cyl_o.cpu().numpy() == 1), + 'loss': float(loss), + 'inter': inter, + 'cyl_points': cyl_points, + 'ovc': ovc, + 'ovb': ovb, + } + side_azlat[s], side_acep[s] = _rel_angles(float(pos[3]), float(pos[4])) + + # 螺絲點雲（存在的側接合）：中心線紅、圓柱 n = L darkcyan / R blue、o = 粉。 + # 點序固定為舊 res_plt_2_torch 的接合序 [L線, R線, L_n, L_o, R_n, R_o]： + # mpl 3D 的 per-point 深度排序對近同深 markers 以輸入序決 tie，換序會讓 + # 同一點集產生亞像素級抗鋸齒邊界差（A/B 實測 ~0.07% 邊緣像素）， + # 固定點序保持與舊圖位元級一致。 + _CYL_COLOR = {'L': 'darkcyan', 'R': 'blue'} + parts_x, parts_y, parts_z, parts_c, parts_s = [], [], [], [], [] + for s in ('L', 'R'): # 中心線 + if s not in side_info: + continue + z_p, y_p, x_p = side_info[s]['line'] + parts_x.append(x_p); parts_y.append(y_p); parts_z.append(z_p) + parts_c.append(_rgba_block(len(x_p), 'r', 1.0)) + parts_s.append(np.full(len(x_p), 3)) + for s in ('L', 'R'): # 圓柱：各側 n 後 o（L_n, L_o, R_n, R_o） + if s not in side_info: + continue + for key, color, size in (('cyl_n', None, 36), ('cyl_o', 'pink', 36)): + z_p, y_p, x_p = side_info[s][key] + parts_x.append(x_p); parts_y.append(y_p); parts_z.append(z_p) + parts_c.append(_rgba_block(len(x_p), _CYL_COLOR[s] if color is None else color, 1.0)) + parts_s.append(np.full(len(x_p), size)) + if parts_x: + x_screw = np.concatenate(parts_x) + y_screw = np.concatenate(parts_y) + z_screw = np.concatenate(parts_z) + screw_rgba = np.concatenate(parts_c) + screw_size = np.concatenate(parts_s) + # ---- 中矢狀（鏡稱）平面 ---- _a, _b, _c, _d = sym["plane"] _n = np.array([_a, _b, _c]) @@ -376,22 +623,41 @@ def render_bone_figure(volume_id, level, binary_path, cortical_path, fig = plt.figure(figsize=(12, 12)) legend_handles = [] + for s in ("L", "R"): + if s in side_info: + legend_handles.append(Line2D([], [], marker="o", ls="", ms=6, + color="darkcyan" if s == "L" else "blue", + label=f"Cylinder({s})")) if vb_corti is not None: legend_handles.append(Line2D([], [], marker="o", ls="", ms=6, color="gold", label="VertebralBody")) if sp_corti is not None: - legend_handles.append(Line2D([], [], marker="o", ls="", ms=6, color="purple", label="Spinous")) + legend_handles.append(Line2D([], [], marker="o", ls="", ms=6, color="purple", label="SpinousProcess")) def _fill_ax(ax): + # 固定分層（關 depth zorder，否則半透明骨頭會被重繪到螺絲上方）： + # 基底骨(5) < VBODY gold(6) < 棘突 purple(6.5) < 終板(7) < 鏡稱面(8) < 螺絲(10) ax.computed_zorder = False - sc_bone = ax.scatter(x_bone, y_bone, z_bone, c=bone_rgba, s=bone_size, marker="o") + sc_bone = ax.scatter(x_base, y_base, z_base, c=rgba_base, s=size_base, marker="o") sc_bone.set_zorder(5) - plane = ax.plot_surface(_Xp, _Yp, _Zp, color="orange", alpha=0.30, - linewidth=1.0, edgecolor="orange", rstride=1, cstride=1) - plane.set_zorder(8) + if x_vb.size: + sc_vb = ax.scatter(x_vb, y_vb, z_vb, c=to_rgba("gold", 0.95), + s=BONE_MARKER_SIZE, marker="o") + sc_vb.set_zorder(6) + if x_sp.size: + sc_sp = ax.scatter(x_sp, y_sp, z_sp, c=to_rgba("purple", 0.95), + s=BONE_MARKER_SIZE, marker="o") + sc_sp.set_zorder(6.5) if _EX is not None: ep = ax.plot_surface(_EX, _EY, _EZ, color="green", alpha=0.35, linewidth=1.0, edgecolor="green", rstride=1, cstride=1) ep.set_zorder(7) + plane = ax.plot_surface(_Xp, _Yp, _Zp, color="orange", alpha=0.30, + linewidth=1.0, edgecolor="orange", rstride=1, cstride=1) + plane.set_zorder(8) + if x_screw is not None: + sc_screw = ax.scatter(x_screw, y_screw, z_screw, + c=screw_rgba, s=screw_size, marker="o") + sc_screw.set_zorder(10) ax1 = fig.add_subplot(221, projection="3d") _fill_ax(ax1) @@ -407,7 +673,8 @@ def render_bone_figure(volume_id, level, binary_path, cortical_path, ax2.legend(handles=legend_handles) ax3 = fig.add_subplot(223, projection="3d") - ax3.view_init(elev=0, azim=90, roll=0) + # 後視圖：相機在 −y 後側，x 軸畫面左小右大 + ax3.view_init(elev=0, azim=-90, roll=0) _fill_ax(ax3) ax3.set_xlabel("X-axis"); ax3.set_ylabel("Y-axis"); ax3.set_zlabel("Z-axis") set_axes_equal_3d(ax3) @@ -419,33 +686,156 @@ def render_bone_figure(volume_id, level, binary_path, cortical_path, set_axes_equal_3d(ax4) label_str = f"{volume_id} {level}" - base_tag = "planes only" if planes_only else "no screws" - if rotation is not None: - base_tag += ", rotated" - fig.text(0.5, 0.98, f"{label_str} ({base_tag})", ha="center", fontsize=15) + if screw_mode: + fig.text(0.5, 0.98, f"{label_str} Best Position", ha="center", fontsize=15) - ratio = float(sym.get("ratio", float("nan"))) - if planes_only: - ep_ratio = float(symp.get("inlier_ratio", float("nan"))) if symp is not None else float("nan") - info = (f"sym_ratio={ratio:.3f} " - f"endplane_ratio={ep_ratio:.3f} " - f"endplate={'yes' if symp is not None else 'no'}") + d_l = f"{diameter_l} mm, {length_l} mm" if diameter_l is not None else '' + d_r = f"{diameter_r} mm, {length_r} mm" if diameter_r is not None else '' + t_s = f"Total time = {total_time:.2f} s" if total_time is not None else '' + fig.text( + 0.5, 0.44, + f"L: Diameter = {d_l}, " + f"R: Diameter = {d_r}, " + f"Swarm size = {swarm_size}, Iteration = {max_iter}, {t_s}", + ha="center", fontsize=12 + ) + + def _ang_segs(az_v, alt_v, azlat_v, acep_v): + # CBT 沒有 2D 參考面，Azimuth/Altitude 直接用最佳化出的原始角； + # TPS 沿用 2D 參考之相對角。終板面擬合失敗（nan）時該段自動略過。 + segs = [f"Azimuth = {az_v:.2f}", f"Altitude = {alt_v:.2f}"] + if np.isfinite(azlat_v): + segs.append(f"Azimuth_Lateral = {azlat_v:.2f}") + if np.isfinite(acep_v): + segs.append(f"Altitude_Endplate = {acep_v:.2f}") + return ', '.join(segs) + + def _side_footer(s, y): + pos, _d, _L, cn_name = side_cfg[s] + if pos is None: + return + if CBT: + segs = _ang_segs(float(pos[3]), float(pos[4]), + side_azlat[s], side_acep[s]) + else: + segs = _ang_segs(90.0 - float(pos[3]) - azi, + 90.0 - float(pos[4]) - alt, + side_azlat[s], side_acep[s]) + di = side_info[s] + cb_ratio = di['ovc'] / di['ovb'] if di['ovb'] else 0.0 + fig.text( + 0.5, y, + f"{cn_name} : Position = ({pos[2]:.2f}, {pos[1]:.2f}, {pos[0]:.2f}), " + f"{segs}, " + f"Intersection = {di['inter']}, " + f"Score = {di['ovc']:.2f} / {di['ovb']:.2f} / {cb_ratio:.2f}", + ha="center", fontsize=8 + ) + + _side_footer('L', 0.03) + _side_footer('R', 0.01) else: - info = (f"sym_ratio={ratio:.3f} " - f"spinous={sp_info.get('n_sp', 0)} ({sp_info.get('mode', '?')}) " - f"vertebral_body={vb_info.get('n_vb', 0)} ({vb_info.get('mode', '?')})") - fig.text(0.5, 0.03, info, ha="center", fontsize=8) + base_tag = "planes only" if planes_only else "no screws" + if rotation is not None: + base_tag += ", rotated" + fig.text(0.5, 0.98, f"{label_str} ({base_tag})", ha="center", fontsize=15) + + ratio = float(sym.get("ratio", float("nan"))) + if planes_only: + ep_ratio = float(symp.get("inlier_ratio", float("nan"))) if symp is not None else float("nan") + info = (f"sym_ratio={ratio:.3f} " + f"endplane_ratio={ep_ratio:.3f} " + f"endplate={'yes' if symp is not None else 'no'}") + else: + info = (f"sym_ratio={ratio:.3f} " + f"spinous={sp_info.get('n_sp', 0)} ({sp_info.get('mode', '?')}) " + f"vertebral_body={vb_info.get('n_vb', 0)} ({vb_info.get('mode', '?')})") + fig.text(0.5, 0.03, info, ha="center", fontsize=8) fig.tight_layout() + date_str = datetime.now().strftime("%Y%m%d") + output_folder = os.path.join(base_folder, date_str, volume_id) + retry_robust(os.makedirs, output_folder, exist_ok=True) + + # ---- CSV（螺絲模式）：L/R 兩行；舊 schema 檔案依欄位名稱重映射後改寫， + # 避免 append 欄位錯位 ---- + if screw_mode and write_csv: + csv_path = os.path.join(output_folder, 'output.csv') + # CBT 模式下恆為 nan 的 2D 參考欄不寫入 CSV + headers = [ + 'Label', 'Side', 'Diameter', 'Length', 'Swarm_Size', 'Max_Iter', + 'Position_XYZ', 'Raw_Azimuth', 'Raw_Altitude', + 'Intersections', 'Best_Loss', 'cyl_points', 'Overlap_Cortical', 'Overlap_Bone', + 'Cortical_Bone_Ratio', + 'Sym_Theta_v_deg', 'Endplate_Tau_y_deg', 'Endplate_Tau_x_deg', + 'Azimuth_Lateral_deg', + 'Altitude_Cephalad_Endplate_deg', + 'Total_Time' + ] + + def _fmt(v): + return '' if not np.isfinite(v) else f"{float(v):.2f}" + + file_exists = os.path.isfile(csv_path) + if file_exists: + with retry_robust(open, csv_path, newline='') as f: + old_rows = [row for row in csv.reader(f) if any(c.strip() for c in row)] + if not old_rows or old_rows[0] != headers: + old_h = old_rows[0] if old_rows else None + with retry_robust(open, csv_path, 'w', newline='') as f: + w = csv.writer(f) + w.writerow(headers) + for r in (old_rows[1:] if old_rows else []): + if old_h: + d = dict(zip(old_h, r)) + w.writerow([d.get(h, '') for h in headers]) + else: + w.writerow(r + [''] * max(0, len(headers) - len(r))) + + try: + with retry_robust(open, csv_path, 'a', newline='') as csvfile: + writer = csv.writer(csvfile) + if not file_exists: + writer.writerow(headers) + for s in ('L', 'R'): + pos, d, L, _cn = side_cfg[s] + if pos is None: + continue + di = side_info[s] + writer.writerow([ + label_str, s, d, L, swarm_size, max_iter, + f"({pos[2]:.2f}, {pos[1]:.2f}, {pos[0]:.2f})", + f"{pos[3]:.2f}", f"{pos[4]:.2f}", + di['inter'], f"{di['loss']:.2f}", di['cyl_points'], + f"{di['ovc']:.2f}", f"{di['ovb']:.2f}", + f"{(di['ovc'] / di['ovb'] if di['ovb'] else 0):.2f}", + _fmt(theta_v), _fmt(tau_y), _fmt(tau_x), + _fmt(side_azlat[s]), _fmt(side_acep[s]), + f"{total_time:.2f}" if total_time is not None else '' + ]) + print(f"[CSV Saved] {csv_path}") + except Exception as e: + print(f"[Error] Failed to write CSV: {e}") + if output_path is not None: path = get_unique_filepath(output_path) + elif screw_mode: + path = save_with_unique_name( + output_folder, level, way, + diameter_l if diameter_l is not None else '', + length_l if length_l is not None else '', + diameter_r if diameter_r is not None else '', + length_r if length_r is not None else '', + swarm_size if swarm_size is not None else '', + max_iter if max_iter is not None else '', + ) else: - date_str = datetime.now().strftime("%Y%m%d") - output_folder = os.path.join(base_folder, date_str, volume_id) - output_file = os.path.join(output_folder, f"{volume_id} {level}_{way}.png") - path = get_unique_filepath(output_file) - os.makedirs(os.path.dirname(path) or ".", exist_ok=True) - fig.savefig(path, dpi=200, bbox_inches="tight") + path = get_unique_filepath( + os.path.join(output_folder, f"{volume_id} {level}_{way}.png")) + + retry_robust(os.makedirs, os.path.dirname(path) or ".", exist_ok=True) + retry_robust(fig.savefig, path, dpi=200, bbox_inches="tight") + print("[Saved figure]", path) plt.close(fig) return path \ No newline at end of file diff --git a/visualization/res_plot_3d.py b/visualization/res_plot_3d.py index 8c9cb94..1158b5b 100644 --- a/visualization/res_plot_3d.py +++ b/visualization/res_plot_3d.py @@ -1,47 +1,14 @@ import torch import numpy as np -import matplotlib.pyplot as plt -from matplotlib.colors import to_rgba -from matplotlib.lines import Line2D -from mpl_toolkits.mplot3d.art3d import Poly3DCollection -import os -import errno -import time -from datetime import datetime -import csv -from core.cylinder import generate_cylinder_n_torch, generate_cylinder_o_torch, snap_to_discrete_values -from core.intersection import center_line_intersections_torch -from core.scoring import cl_score_torch, compute_overlap_ratio_from_cylinder_mask, cl_score_torch_xfr -from imaging.orientation import (azimuth_rotation, analyze_vertebral_tilt_contour, - best_symmetry_plane, best_upper_endplate_plane, - segment_spinous_process, segment_vertebral_body) -from utils.helpers import save_with_unique_name +from core.cylinder import generate_cylinder_n_torch, snap_to_discrete_values +from core.scoring import compute_overlap_ratio_from_cylinder_mask -# Volume absorption 渲染（Beer-Lambert）：每 voxel 不透明度 = 1 - exp(-mu * voxel_width) -# 骨頭核心厚度達 70-90 voxel，沿視線堆疊會使任何 per-voxel alpha 累積成不透明。 -# 因此以「抽稀 (SUBSAMPLE) 降低堆疊數量」+「低 mu 控制每點吸收」兩項共同調出淡薄 X-ray 陰影， -# 同時保留皮質 / 鬆質的吸收入射差異（mu 比值）。 -BONE_MU_CORTICAL = 0.02 # 1/mm → 每 voxel = 1-exp(-0.02*0.5) ~ 0.010 -BONE_MU_TRABECULAR = 0.005 # 1/mm → 每 voxel = 1-exp(-0.005*0.5) ~ 0.0025 -BONE_MARKER_SIZE = 3.0 # 骨骼散點點面積 (pt^2)；略大以補償抽稀後的顆粒感 -BONE_SUBSAMPLE = 1 # 每 10 個骨 voxel 畫 1 個，降低堆疊不透明度（1=全畫） - -def _retry_robust(fn, *args, retries=20, delay=0.5, **kwargs): - """對 ENOENT/EEXIST 重試：NFS 上輸出樹被外部刪除（或多 worker 併發建同一 - output 目錄）會有短暫的 ENOENT 窗口，重試可恢復；其他錯誤直接丟出。""" - for i in range(retries): - try: - return fn(*args, **kwargs) - except OSError as e: - if e.errno not in (errno.ENOENT, errno.EEXIST) or i == retries - 1: - raise - time.sleep(delay) def set_axes_equal_3d(ax): """ Make axes of 3D plot have equal scale so that spheres appear as spheres, - cubes as cubes, etc. + cubes are cubes, etc. """ x_limits = ax.get_xlim3d() y_limits = ax.get_ylim3d() @@ -59,610 +26,31 @@ def set_axes_equal_3d(ax): ax.set_xlim3d([x_middle - plot_radius, x_middle + plot_radius]) ax.set_ylim3d([y_middle - plot_radius, y_middle + plot_radius]) ax.set_zlim3d([z_middle - plot_radius, z_middle + plot_radius]) - + try: ax.set_box_aspect([1, 1, 1]) except AttributeError: pass -def res_plt_2_torch( - spine_tensor: torch.Tensor, - cortical_tensor: torch.Tensor, - image_shape: tuple[int, int, int], - image2_path: str, - base_folder: str, - label_str: str, - diameter_l: float, - length_l: float, - diameter_r: float, - length_r: float, - best_position_l: list[float], - best_position_r: list[float], - swarm_size: int, - max_iter: int, - total_time: float, - spacing: list[float], - CBT: bool, - device: torch.device, - grid=None -) -> None: - """ - Same plotting function as before, but it uses torch-based generation - and then moves data to CPU for matplotlib 3D scatter. - """ - cyl_l = generate_cylinder_n_torch( - diameter_l, - length_l, - best_position_l[0], - best_position_l[1], - best_position_l[2], - best_position_l[3], - best_position_l[4], - image_shape, - spacing, - device, - grid + +def res_plt_2_torch(spine_tensor, cortical_tensor, image_shape, image2_path, + base_folder, label_str, diameter_l, length_l, + diameter_r, length_r, best_position_l, best_position_r, + swarm_size, max_iter, total_time, spacing, CBT, device, + grid): + """相容入口（舊簽名）：實作已合併進 res_bone_figure.render_bone_figure。 + label_str / image_shape 不再使用（volume / level 由 image2_path 反推）。""" + from visualization.res_bone_figure import render_bone_figure + return render_bone_figure( + None, None, spine_tensor, cortical_tensor, base_folder, + spacing=spacing, way='CBT' if CBT else 'TPS', + best_position_l=best_position_l, best_position_r=best_position_r, + diameter_l=diameter_l, length_l=length_l, + diameter_r=diameter_r, length_r=length_r, + image2_path=image2_path, device=device, grid=grid, + swarm_size=swarm_size, max_iter=max_iter, total_time=total_time, ) - cyl_lo = generate_cylinder_o_torch( - diameter_l, - length_l, - best_position_l[0], - best_position_l[1], - best_position_l[2], - best_position_l[3], - best_position_l[4], - image_shape, - spacing, - device, - grid - ) - cyl_r = generate_cylinder_n_torch( - diameter_r, - length_r, - best_position_r[0], - best_position_r[1], - best_position_r[2], - best_position_r[3], - best_position_r[4], - image_shape, - spacing, - device, - grid - ) - cyl_ro = generate_cylinder_o_torch( - diameter_r, - length_r, - best_position_r[0], - best_position_r[1], - best_position_r[2], - best_position_r[3], - best_position_r[4], - image_shape, - spacing, - device, - grid - ) - - intersections_l, line_mask_l = center_line_intersections_torch( - best_position_l[0], - best_position_l[1], - best_position_l[2], - best_position_l[3], - best_position_l[4], - int(length_l), - spine_tensor, - spacing, - device - ) - loss_l = cl_score_torch(cortical_tensor, spine_tensor, cyl_l, cyl_lo, intersections_l) - - intersections_r, line_mask_r = center_line_intersections_torch( - best_position_r[0], - best_position_r[1], - best_position_r[2], - best_position_r[3], - best_position_r[4], - int(length_r), - spine_tensor, - spacing, - device - ) - # loss_r = cl_score_torch(cortical_tensor, spine_tensor, cyl_r, cyl_ro, intersections_r) - # loss_r 放在下方 VBODY mask 計算之後：計入與 PSO 目標函數相同的 VBODY voxel 獎勵 - - if CBT: - azi = float('nan') - alt = float('nan') - else: - azi = azimuth_rotation(image2_path) - res = analyze_vertebral_tilt_contour(image2_path, edge_type='superior', show_plot=False, debug=False) - alt = res['superior']['tilt_angle_deg'] - - # Move data to CPU for plotting - line_mask_l_cpu = line_mask_l.cpu().numpy() - line_mask_r_cpu = line_mask_r.cpu().numpy() - cyl_l_cpu = cyl_l.cpu().numpy() - cyl_lo_cpu = cyl_lo.cpu().numpy() - cyl_r_cpu = cyl_r.cpu().numpy() - cyl_ro_cpu = cyl_ro.cpu().numpy() - spine_cpu = spine_tensor.cpu().numpy() - - z_lin1, y_lin1, x_lin1 = np.where(line_mask_l_cpu == 1) - z_lin2, y_lin2, x_lin2 = np.where(line_mask_r_cpu == 1) - - z_cyl_l1, y_cyl_l1, x_cyl_l1 = np.where(cyl_l_cpu == 1) - z_cyl_l2, y_cyl_l2, x_cyl_l2 = np.where(cyl_lo_cpu == 1) - z_cyl_r1, y_cyl_r1, x_cyl_r1 = np.where(cyl_r_cpu == 1) - z_cyl_r2, y_cyl_r2, x_cyl_r2 = np.where(cyl_ro_cpu == 1) - - # 骨頭 voxel 依「體積吸收」分成兩組：皮質（高不透明度）與鬆質（低不透明度） - cortical_cpu = cortical_tensor.cpu().numpy() - voxel_mm = float(spacing[0]) - alpha_cortical = 1.0 - np.exp(-BONE_MU_CORTICAL * voxel_mm) - alpha_trabecular = 1.0 - np.exp(-BONE_MU_TRABECULAR * voxel_mm) - z_corti, y_corti, x_corti = np.where((spine_cpu == 1) & (cortical_cpu == 1)) - z_trab, y_trab, x_trab = np.where((spine_cpu == 1) & (cortical_cpu == 0)) - - # 中矢狀面：骨頭的最佳鏡稱面，一般平面 a·x + b·y + c·z = d（法線方向任意） - sym = best_symmetry_plane(spine_cpu) - - # 棘突：鏡稱面中線帶（|s|<=w）且在 AP 谷底之後側的骨 voxel，換不同顏色標示 - # 棘突缺如（先前 laminectomy / 棘突切除）時 sp_mask=None，不標示。 - sp_mask, sp_th, sp_info = segment_spinous_process(spine_cpu, sym) - if sp_info['mode'] == 'no_spinous': - top_off = f"{sp_info['top_off']:.1f}" if sp_info.get('top_off') is not None else 'n/a' - print(f"[NO-SP] 中線後側缺如（先前 laminectomy / 棘突切除）: " - f"deficit={sp_info['deficit']:.1f} voxel ({sp_info['deficit'] * 0.5:.1f} mm), " - f"rear3={sp_info['rear3']} voxel, top_off={top_off} voxel " - f"-> 不標示棘突；椎體用放寬後側谷底切分") - sp_corti = sp_trab = None - elif sp_mask is not None and sp_mask.any(): - sp_corti = sp_mask[z_corti, y_corti, x_corti] - sp_trab = sp_mask[z_trab, y_trab, x_trab] - sp_n_bone = max(int(spine_cpu.sum()), 1) - print(f"[SPINOUS] n={sp_info['n_sp']} " - f"({100.0 * sp_info['n_sp'] / sp_n_bone:.1f}% of bone) " - f"band=+/-{sp_info['band_w']:.1f} voxel AP>={sp_info['ap_thresh']:.1f} " - f"mode={sp_info['mode']}") - else: - sp_corti = sp_trab = None - - # 上終板平面：RANSAC 擬合骨頭頂面（前側）的最佳 a·x + b·y + c·z = d - symp = best_upper_endplate_plane(spine_cpu) - - # 椎體：上終板之下（排除跨終板的後側構造）且中線 AP 谷底之前側的骨 voxel， - # 換不同顏色標示（見 segment_vertebral_body；谷底優先取中線帶， - # 中線搜尋 fallback 時退回終板下整體 AP 分佈谷底） - vb_mask, vb_th, vb_info = segment_vertebral_body(spine_cpu, sym, symp, - sp_th, sp_info['mode']) - if vb_mask is not None and vb_mask.any(): - vb_corti = vb_mask[z_corti, y_corti, x_corti] - vb_trab = vb_mask[z_trab, y_trab, x_trab] - print(f"[VBODY] n={vb_info['n_vb']} " - f"({100.0 * vb_info['n_vb'] / max(int(spine_cpu.sum()), 1):.1f}% of bone) " - f"AP<{vb_info['ap_thresh']:.1f} mode={vb_info['mode']}") - if vb_info['mode'] == 'quantile': - print(f"[VBODY] WARNING: 未找到體/弓後側谷底，閾值退回 55 百分位 " - f"（可能切進椎體內），建議人工核對該 level 的椎體邊界") - else: - vb_corti = vb_trab = None - print(f"[VBODY] skipped: {vb_info['mode']}") - - # loss_r 使用與 PSO 目標函數相同的 VBODY 獎勵（回報分數與優化一致） - vbody_tensor = None - if vb_mask is not None and vb_mask.any(): - vbody_tensor = torch.from_numpy(vb_mask.astype(np.uint8)).to(device=device) - loss_r = cl_score_torch_xfr(cortical_tensor, spine_tensor, cyl_r, cyl_ro, intersections_r, - vbody_tensor=vbody_tensor) - - # X-ray 外觀：骨頭合成一個半透明體積吸收點雲（下方）； - # 螺絲（中心線 + 圓柱 + 入口軌跡延长）合成一個全不透明點雲，永遠畫在骨頭之上 - def _rgba_block(n, color, a): - arr = np.empty((n, 4)) - arr[:] = to_rgba(color) - arr[:, 3] = a - return arr - - x_bone = np.concatenate([x_corti, x_trab]) - y_bone = np.concatenate([y_corti, y_trab]) - z_bone = np.concatenate([z_corti, z_trab]) - bone_rgba = np.concatenate([ - _rgba_block(len(x_corti), 'lightblue', float(alpha_cortical)), - _rgba_block(len(x_trab), 'lightblue', float(alpha_trabecular)), - ]) - bone_size = np.full(len(x_bone), BONE_MARKER_SIZE) - - # 抽稀：降低堆疊不透明度以呈現淡薄 X-ray 陰影 - if BONE_SUBSAMPLE > 1: - x_bone = x_bone[::BONE_SUBSAMPLE] - y_bone = y_bone[::BONE_SUBSAMPLE] - z_bone = z_bone[::BONE_SUBSAMPLE] - bone_rgba = bone_rgba[::BONE_SUBSAMPLE] - bone_size = bone_size[::BONE_SUBSAMPLE] - # 平面 patch 的範圍用「完整骨頭」（含 VBODY / SP）算： - # VBODY 前側是整顆骨最前緣，剔除後綠色終板 patch 會縮小 - x_bone_all, y_bone_all, z_bone_all = x_bone, y_bone, z_bone - # VBODY / 棘突拆成獨立上層（在 _fill_ax 內畫）： - # 繪製順序 基底骨(5) < VBODY gold(6) < SP purple(6.5) < 終板(7) < 鏡稱面(8) < 螺絲(10) - # axial 視角（ax3）相機在 +y 前側，椎體 physically 擋在棘突與相機之間， - # 棘突最後畫 -> 紫色不被 gold 遮住 - vb_flag = np.zeros(x_bone.shape, dtype=bool) - sp_flag = np.zeros(x_bone.shape, dtype=bool) - if vb_corti is not None: - f = np.concatenate([vb_corti, vb_trab]).astype(bool) - if BONE_SUBSAMPLE > 1: - f = f[::BONE_SUBSAMPLE] - vb_flag |= f - if sp_corti is not None: - f = np.concatenate([sp_corti, sp_trab]).astype(bool) - if BONE_SUBSAMPLE > 1: - f = f[::BONE_SUBSAMPLE] - sp_flag |= f - - vbody_pts = None - vbody_mask = vb_flag & ~sp_flag - if vbody_mask.any(): - vbody_pts = (x_bone[vbody_mask], y_bone[vbody_mask], z_bone[vbody_mask]) - sp_pts = None - if sp_flag.any(): - sp_pts = (x_bone[sp_flag], y_bone[sp_flag], z_bone[sp_flag]) - - # 基底骨層去掉 VBODY / SP voxel（由上面的專屬圖層畫） - base_mask = ~(vb_flag | sp_flag) - x_bone = x_bone[base_mask] - y_bone = y_bone[base_mask] - z_bone = z_bone[base_mask] - bone_rgba = bone_rgba[base_mask] - bone_size = bone_size[base_mask] - - x_screw = np.concatenate([x_lin1, x_lin2, x_cyl_l1, x_cyl_l2, x_cyl_r1, x_cyl_r2]) - y_screw = np.concatenate([y_lin1, y_lin2, y_cyl_l1, y_cyl_l2, y_cyl_r1, y_cyl_r2]) - z_screw = np.concatenate([z_lin1, z_lin2, z_cyl_l1, z_cyl_l2, z_cyl_r1, z_cyl_r2]) - - _a, _b, _c, _d = sym['plane'] - _n = np.array([_a, _b, _c]) - _u = np.array(sym['u']) - _v = np.array(sym['v']) - _p0 = _d * _n # 平面上最接近原點的點 - xyz_bone = np.stack([x_bone_all - _p0[0], y_bone_all - _p0[1], z_bone_all - _p0[2]], axis=1) - _pu = xyz_bone @ _u - _pv = xyz_bone @ _v - _U_, _V_ = np.meshgrid(np.linspace(_pu.min(), _pu.max(), 8), - np.linspace(_pv.min(), _pv.max(), 8)) - _Xp = _p0[0] + _U_ * _u[0] + _V_ * _v[0] - _Yp = _p0[1] + _U_ * _u[1] + _V_ * _v[1] - _Zp = _p0[2] + _U_ * _u[2] + _V_ * _v[2] - - _EX = _EY = _EZ = None - if symp is not None: - _ea, _eb, _ec, _ed = symp['plane'] - _en = np.array([_ea, _eb, _ec]) - _eu = np.array(symp['u']) - _ev = np.array(symp['v']) - _ep0 = _ed * _en - xz_ep = np.stack([x_bone_all - _ep0[0], y_bone_all - _ep0[1], z_bone_all - _ep0[2]], axis=1) - _pu_ep = xz_ep @ _eu - _pv_ep = xz_ep @ _ev - _EU, _EV = np.meshgrid(np.linspace(_pu_ep.min(), _pu_ep.max(), 8), - np.linspace(_pv_ep.min(), _pv_ep.max(), 8)) - _EX = _ep0[0] + _EU * _eu[0] + _EV * _ev[0] - _EY = _ep0[1] + _EU * _eu[1] + _EV * _ev[1] - _EZ = _ep0[2] + _EU * _eu[2] + _EV * _ev[2] - screw_rgba = np.concatenate([ - _rgba_block(len(x_lin1), 'r', 1.0), - _rgba_block(len(x_lin2), 'r', 1.0), - _rgba_block(len(x_cyl_l1), 'darkcyan', 1.0), - _rgba_block(len(x_cyl_l2), 'pink', 1.0), - _rgba_block(len(x_cyl_r1), 'blue', 1.0), - _rgba_block(len(x_cyl_r2), 'pink', 1.0), - ]) - screw_size = np.concatenate([ - np.full(len(x_lin1), 3), np.full(len(x_lin2), 3), - np.full(len(x_cyl_l1), 36), np.full(len(x_cyl_l2), 36), - np.full(len(x_cyl_r1), 36), np.full(len(x_cyl_r2), 36), - ]) - - fig = plt.figure(figsize=(12, 12)) - - legend_handles = [ - Line2D([], [], marker='o', ls='', ms=6, color='darkcyan', label='Cylinder(L)'), - Line2D([], [], marker='o', ls='', ms=6, color='blue', label='Cylinder(R)'), - ] - if vb_corti is not None: - legend_handles.append( - Line2D([], [], marker='o', ls='', ms=6, color='gold', label='VertebralBody')) - if sp_corti is not None: - legend_handles.append( - Line2D([], [], marker='o', ls='', ms=6, color='purple', label='SpinousProcess')) - - def _fill_ax(ax): - # X-ray 外觀：關閉 mplot3d 依深度自動排序 zorder（否則半透明骨頭會被重繪到 - # 螺絲上方）；改為固定分層： - # 基底骨 zorder=5 < VBODY 6 < SP 6.5 < 終板 7 < 鏡稱面 8 < 螺絲 10 - ax.computed_zorder = False - sc_bone = ax.scatter(x_bone, y_bone, z_bone, c=bone_rgba, s=bone_size, marker='o') - sc_bone.set_zorder(5) - if vbody_pts is not None: - sc_vb = ax.scatter(vbody_pts[0], vbody_pts[1], vbody_pts[2], - c=to_rgba('gold', 0.95), s=BONE_MARKER_SIZE, marker='o') - sc_vb.set_zorder(6) - if sp_pts is not None: - sc_sp = ax.scatter(sp_pts[0], sp_pts[1], sp_pts[2], - c=to_rgba('purple', 0.95), s=BONE_MARKER_SIZE, marker='o') - sc_sp.set_zorder(6.5) - sc_screw = ax.scatter(x_screw, y_screw, z_screw, c=screw_rgba, s=screw_size, marker='o') - sc_screw.set_zorder(10) - # 中矢狀面（理論左右對稱切分面）：半透明橘色平面 x = x_mid - # 平面邊緣畫橘色線，讓 axial / 正視（側看時）也能清楚看到切分線 - plane = ax.plot_surface(_Xp, _Yp, _Zp, color='orange', alpha=0.30, - linewidth=1.0, edgecolor='orange', rstride=1, cstride=1) - plane.set_zorder(8) - # 上終板平面：半透明綠色平面（邊緣綠線） - if _EX is not None: - ep = ax.plot_surface(_EX, _EY, _EZ, color='green', alpha=0.35, - linewidth=1.0, edgecolor='green', rstride=1, cstride=1) - ep.set_zorder(7) - - ax1 = fig.add_subplot(221, projection='3d') - _fill_ax(ax1) - ax1.set_xlabel('X-axis'); ax1.set_ylabel('Y-axis'); ax1.set_zlabel('Z-axis') - set_axes_equal_3d(ax1) - - ax2 = fig.add_subplot(222, projection='3d') - ax2.view_init(elev=90, azim=-90, roll=0) - _fill_ax(ax2) - ax2.set_xlabel('X-axis'); ax2.set_ylabel('Y-axis'); ax2.set_zlabel('Z-axis') - set_axes_equal_3d(ax2) - ax2.legend(handles=legend_handles) - - ax3 = fig.add_subplot(223, projection='3d') - ax3.view_init(elev=0, azim=90, roll=0) - _fill_ax(ax3) - ax3.set_xlabel('X-axis'); ax3.set_ylabel('Y-axis'); ax3.set_zlabel('Z-axis') - set_axes_equal_3d(ax3) - - ax4 = fig.add_subplot(224, projection='3d') - ax4.view_init(elev=0, azim=0, roll=0) - _fill_ax(ax4) - ax4.set_xlabel('X-axis'); ax4.set_ylabel('Y-axis'); ax4.set_zlabel('Z-axis') - set_axes_equal_3d(ax4) - - cyl_points_l = torch.sum(cyl_l).item() - cyl_points_r = torch.sum(cyl_r).item() - - overlap_l = ((cortical_tensor == 1) & (cyl_l == 1)).sum().item() - overlap_r = ((cortical_tensor == 1) & (cyl_r == 1)).sum().item() - overlap_b_l = ((spine_tensor == 1) & (cyl_l == 1)).sum().item() - overlap_b_r = ((spine_tensor == 1) & (cyl_r == 1)).sum().item() - - overlap_cortical_l = (overlap_l / cyl_points_l) * 100 if cyl_points_l else 0.0 - overlap_cortical_r = (overlap_r / cyl_points_r) * 100 if cyl_points_r else 0.0 - overlap_vertebral_l = (overlap_b_l / cyl_points_l) * 100 if cyl_points_l else 0.0 - overlap_vertebral_r = (overlap_b_r / cyl_points_r) * 100 if cyl_points_r else 0.0 - cb_ratio_l = overlap_cortical_l/overlap_vertebral_l if overlap_vertebral_l else 0.0 - cb_ratio_r = overlap_cortical_r/overlap_vertebral_r if overlap_vertebral_r else 0.0 - user_altitude_l = 90 - best_position_l[4] - alt - user_altitude_r = 90 - best_position_r[4] - alt - user_azimuth_l = 90 - best_position_l[3] - azi - user_azimuth_r = 90 - best_position_r[3] - azi - - # 螺絲方向向量（與 generate_cylinder_n_torch 同慣例）： - # d = (cos(az)·sin(alt), sin(az)·sin(alt), cos(alt))，alt = 相对 +z 的極角 - # Azimuth 相对鏡稱面（法線 s，theta_v = atan2(sy, sx)）： - # Azimuth_Lateral = az - theta_v - 90 （面內 AP 軸起的帶號發散角，+ = L 側往外，− = R 側） - # Altitude 相對上終板面（法線 e，朝上）： - # Altitude_Cephalad_Endplate = 90 - (d 與 e 的夾角) - sym_n = np.asarray(sym['normal'], dtype=float) - theta_v = float(np.degrees(np.arctan2(sym_n[1], sym_n[0]))) - if symp is not None: - e_n = np.asarray(symp['normal'], dtype=float) - e_n = e_n / np.linalg.norm(e_n) - tau_y = float(np.degrees(np.arctan2(e_n[1], e_n[2]))) - tau_x = float(np.degrees(np.arctan2(e_n[0], e_n[2]))) - else: - e_n = None - tau_y = float('nan') - tau_x = float('nan') - - def _rel_angles(az_deg, alt_deg): - az_r = np.radians(az_deg) - alt_r = np.radians(alt_deg) - d_v = np.array([np.cos(az_r) * np.sin(alt_r), - np.sin(az_r) * np.sin(alt_r), - np.cos(alt_r)]) - az_lateral = az_deg - theta_v - 90.0 - if e_n is not None: - ang_norm = float(np.degrees(np.arccos(np.clip(d_v @ e_n, -1.0, 1.0)))) - alt_cep = 90.0 - ang_norm - else: - alt_cep = float('nan') - return az_lateral, alt_cep - - azlat_l, acep_l = _rel_angles(best_position_l[3], best_position_l[4]) - azlat_r, acep_r = _rel_angles(best_position_r[3], best_position_r[4]) - - date_str = datetime.now().strftime("%Y%m%d") - # 旋轉後影像存在 /rotated/ 下：parent 是 'rotated'， - # 再上一層才是 volume id（未旋轉路徑不受影響） - img_parent = os.path.dirname(image2_path) - patient_id = os.path.basename(os.path.dirname(img_parent)) \ - if os.path.basename(img_parent) == 'rotated' else os.path.basename(img_parent) - output_folder = os.path.join(base_folder, date_str, patient_id) - _retry_robust(os.makedirs, output_folder, exist_ok=True) - csv_path = os.path.join(output_folder, 'output.csv') - - # 欄位標題 (Header)。CBT 模式下恆為 nan 的 2D 參考欄 - # (Azimuth_Diff / Altitude_Diff / User_Azimuth / User_Altitude) 不寫入 CSV。 - headers = [ - 'Label', 'Side', 'Diameter', 'Length', 'Swarm_Size', 'Max_Iter', - 'Position_XYZ', 'Raw_Azimuth', 'Raw_Altitude', - 'Intersections', 'Best_Loss', 'cyl_points', 'Overlap_Cortical', 'Overlap_Bone', - 'Cortical_Bone_Ratio', - 'Sym_Theta_v_deg', 'Endplate_Tau_y_deg', 'Endplate_Tau_x_deg', - 'Azimuth_Lateral_deg', - 'Altitude_Cephalad_Endplate_deg', - 'Total_Time' - ] - - def _fmt(v): - return '' if not np.isfinite(v) else f"{float(v):.2f}" - - # 檢查檔案是否存在 (決定是否寫入標題)；舊 schema 的檔案按欄位名稱重映射後 - # 改以新 Header 重寫（舊檔多出的欄位捨去、缺的欄位補空白），避免 append 欄位錯位 - file_exists = os.path.isfile(csv_path) - if file_exists: - with _retry_robust(open, csv_path, newline='') as f: - old_rows = [row for row in csv.reader(f) if any(c.strip() for c in row)] - if not old_rows or old_rows[0] != headers: - old_h = old_rows[0] if old_rows else None - with _retry_robust(open, csv_path, 'w', newline='') as f: - w = csv.writer(f) - w.writerow(headers) - for r in (old_rows[1:] if old_rows else []): - if old_h: - d = dict(zip(old_h, r)) - w.writerow([d.get(h, '') for h in headers]) - else: - w.writerow(r + [''] * max(0, len(headers) - len(r))) - - try: - with _retry_robust(open, csv_path, 'a', newline='') as csvfile: - writer = csv.writer(csvfile) - - # 新檔案寫入 Header - if not file_exists: - writer.writerow(headers) - - # 寫入 Left 數據 - writer.writerow([ - label_str, - 'L', - diameter_l, - length_l, - swarm_size, - max_iter, - # f"({best_position_l[0]:.2f}, {best_position_l[1]:.2f}, {best_position_l[2]:.2f})", - f"({best_position_l[2]:.2f}, {best_position_l[1]:.2f}, {best_position_l[0]:.2f})", - f"{best_position_l[3]:.2f}", - f"{best_position_l[4]:.2f}", - intersections_l, - f"{loss_l:.2f}", - cyl_points_l, - f"{overlap_cortical_l:.2f}", - f"{overlap_vertebral_l:.2f}", - f"{(overlap_cortical_l/overlap_vertebral_l if overlap_vertebral_l!=0 else 0):.2f}", - _fmt(theta_v), - _fmt(tau_y), - _fmt(tau_x), - _fmt(azlat_l), - _fmt(acep_l), - f"{total_time:.2f}" - ]) - - # 寫入 Right 數據 - writer.writerow([ - label_str, - 'R', - diameter_r, - length_r, - swarm_size, - max_iter, - # f"({best_position_r[0]:.2f}, {best_position_r[1]:.2f}, {best_position_r[2]:.2f})", - f"({best_position_r[2]:.2f}, {best_position_r[1]:.2f}, {best_position_r[0]:.2f})", - f"{best_position_r[3]:.2f}", - f"{best_position_r[4]:.2f}", - intersections_r, - f"{loss_r:.2f}", - cyl_points_r, - f"{overlap_cortical_r:.2f}", - f"{overlap_vertebral_r:.2f}", - f"{(overlap_cortical_r/overlap_vertebral_r if overlap_vertebral_r!=0 else 0):.2f}", - _fmt(theta_v), - _fmt(tau_y), - _fmt(tau_x), - _fmt(azlat_r), - _fmt(acep_r), - f"{total_time:.2f}" - ]) - print(f"[CSV Saved] {csv_path}") - - except Exception as e: - print(f"[Error] Failed to write CSV: {e}") - - fig.text(0.5, 0.98, f'{label_str} Best Position', ha='center', fontsize=15) - fig.text( - 0.5, 0.44, - f'L: Diameter = {diameter_l} mm, {length_l} mm, ' - f'R: Diameter = {diameter_r} mm, {length_r} mm, ' - f'Swarm size = {swarm_size}, Iteration = {max_iter}, Total time = {total_time:.2f} s', - ha='center', fontsize=12 - ) - - # 角度註記：CBT 沒有 2D 參考面（user az/alt 為 nan），Azimuth/Altitude 直接顯示 - # 最佳化出的原始角（＝CSV 的 Raw_Azimuth / Raw_Altitude）；TPS 沿用 2D 參考之相對角。 - # 另補上相對骨骼的角度（與 CSV 同參數）： - # Azimuth_Lateral = 螺絲在鏡稱面內相對 AP 軸的發散角（+ = L 側往外，− = R 側） - # Altitude_Endplate = 螺絲相對上終板面的仰角 - # 終板面擬合失敗（nan）時該段自動略過。 - def _fig_angle_segs(az, alt, azlat, acep): - segs = [f'Azimuth = {az:.2f}', f'Altitude = {alt:.2f}'] - if np.isfinite(azlat): - segs.append(f'Azimuth_Lateral = {azlat:.2f}') - if np.isfinite(acep): - segs.append(f'Altitude_Endplate = {acep:.2f}') - return ', '.join(segs) - - if CBT: - _ang_l = _fig_angle_segs(float(best_position_l[3]), float(best_position_l[4]), azlat_l, acep_l) - _ang_r = _fig_angle_segs(float(best_position_r[3]), float(best_position_r[4]), azlat_r, acep_r) - else: - _ang_l = _fig_angle_segs(user_azimuth_l, user_altitude_l, azlat_l, acep_l) - _ang_r = _fig_angle_segs(user_azimuth_r, user_altitude_r, azlat_r, acep_r) - - fig.text( - 0.5, 0.03, - f'Left : Position = ({best_position_l[2]:.2f}, {best_position_l[1]:.2f}, {best_position_l[0]:.2f}), ' - f'{_ang_l}, ' - f'Intersection = {intersections_l}, Score = {overlap_cortical_l:.2f} / {overlap_vertebral_l:.2f} / {cb_ratio_l:.2f}', - ha='center', fontsize=8 - ) - fig.text( - 0.5, 0.01, - f'Right : Position = ({best_position_r[2]:.2f}, {best_position_r[1]:.2f}, {best_position_r[0]:.2f}), ' - f'{_ang_r}, ' - f'Intersection = {intersections_r}, Score = {overlap_cortical_r:.2f} / {overlap_vertebral_r:.2f} / {cb_ratio_r:.2f}', - ha='center', fontsize=8 - ) - - fig.tight_layout() - - date_str = datetime.now().strftime("%Y%m%d") - file_name = os.path.basename(image2_path) - level = file_name.split('_')[0] - output_folder = os.path.join(base_folder, date_str, patient_id) - - if CBT == True: - way = 'CBT' - - else: - way = 'TPS' - - # 建目錄 + 存檔一起重試：輸出樹被外部刪除（NFS 刪除競態）時，重試會重建目錄 - path = None - def _save_fig_once(): - nonlocal path - _retry_robust(os.makedirs, output_folder, exist_ok=True) - # 檔名只用 level（volume id 已在資料夾名裡，不重複） - path = save_with_unique_name(output_folder, level, way, - diameter_l, length_l, diameter_r, length_r, - swarm_size, max_iter) - fig.savefig(path, dpi=200, bbox_inches="tight") - - _retry_robust(_save_fig_once) - print("[Saved figure]", path) - plt.close(fig) - def eval_overlap_from_position( pos, @@ -700,6 +88,4 @@ def eval_overlap_from_position( ) overlap = compute_overlap_ratio_from_cylinder_mask(cyl_mask, spine_tensor) - return overlap, d, L - - + return overlap, d, L \ No newline at end of file diff --git a/xfr_preprocess.py b/xfr_preprocess.py index 73e50b6..6db50cd 100644 --- a/xfr_preprocess.py +++ b/xfr_preprocess.py @@ -177,7 +177,8 @@ def _write_rotated_level(vol_dir, level, smd_path, mask_path, roi_path): - _roi.nii.gz - _cortical.nii.gz 旋轉 CT 以骨頭 mask 內 median HU 為門檻 （取代舊的未旋轉 _cortical，定義相同） - 再畫「rotated」平面圖，並在旋轉體上做 VBODY / 棘突分割存 label map + 再畫「rotated」平面圖（bone + 平面 + VBODY [金] / 棘突 [紫] 著色）， + 並在旋轉體上做 VBODY / 棘突分割存 label map （1=VBODY、2=棘突、3=other bone、0=background）。 平面參數經 R 剛性旋轉並換元到輸出 grid 的局部座標。 @@ -316,13 +317,14 @@ def _write_rotated_level(vol_dir, level, smd_path, mask_path, roi_path): else: logger.warning(f'[rotated] {volume_id} {level}: 無旋轉 CT / mask，跳過 _cortical') - # 用旋轉後的平面畫圖（rotated 版 planes）；皮質著色由未旋轉 CT + mask + # 用旋轉後的平面畫圖（rotated 版 planes），並畫出 VBODY / 棘突著色 + #（VBODY=金、棘突=紫，與 label map 對應）；皮質著色由未旋轉 CT + mask # 現算（未旋轉 _cortical 不再存檔） p_fig = os.path.join(rotated_dir, f'{level}_planes.png') if mask_path is not None: fig_cortical = _cortical_from_roi(roi_arr, bin_arr) fig = render_bone_figure(volume_id, level, mask_path, fig_cortical, - planes_only=True, rotation=(R, c_xyz), output_path=p_fig) + planes_only=False, rotation=(R, c_xyz), output_path=p_fig) if fig is not None: logger.info(f'[rotated] saved {fig}') else: @@ -357,7 +359,7 @@ def make_lumbar_post_process(): """每個 volume 處理完後，對其 lumbar level： 1) 畫「骨頭 + 方向平面」圖（不畫螺絲、不做棘突 / 椎體分割）-> /lumbar/ 2) 計算對齊旋轉，存旋轉後的 smd_resampled / binary_sdf / binary_nn / roi - + cortical + 旋轉平面圖 + label map -> /rotated/ + + cortical + 旋轉平面圖（含 VBODY / 棘突著色）+ label map -> /rotated/ post_process 由 process_dataset 呼叫：(volume_dir, processed_labels)。 processed_labels 為該 volume 實際存在的 label id（int），對照 LABEL_MAP。 @@ -397,7 +399,7 @@ def make_lumbar_post_process(): logger.info(f'[lumbar] saved {path}') # 2) 旋轉對齊：rotated/ 的 smd_resampled + binary_sdf + binary_nn - # + roi + cortical + planes 圖 + label + # + roi + cortical + planes 圖（含 VBODY / 棘突著色）+ label _write_rotated_level(vol_dir, level, smd_res_path, mask_path, roi_path) return _post_process @@ -407,7 +409,12 @@ def main(): parser = argparse.ArgumentParser(description='Preprocess CT spine dataset.') parser.add_argument('--max-images', type=int, default=None, dest='max_images', help='Process at most this number of images per dataset (default: all).') + parser.add_argument('--output-dir', type=str, default=None, dest='output_dir', + help='Override the default output dir (e.g. repair run on an ' + 'older generation). Default: the module-level output_dir.') args = parser.parse_args() + # local alias（避免 rebind module-level output_dir 造成 UnboundLocalError） + out_dir = args.output_dir if args.output_dir is not None else output_dir # log 檔（console 與檔案同時輸出） os.makedirs(LOG_DIR, exist_ok=True) @@ -422,6 +429,7 @@ def main(): logger.info(f'Log file: {log_path}') logger.info(f'Command: {sys.executable} {" ".join(sys.argv)}') logger.info(f'Working directory: {os.getcwd()}') + logger.info(f'Output dir: {out_dir}') # metadata db：跳過判定（z spacing / lumbar 層數）命中時免讀影像 / label 檔 metadata_db = ImageMetadataDB() @@ -432,7 +440,7 @@ def main(): for key, value in label_map.items(): data_dir = os.path.join(data_root, key) label_dir = os.path.join(label_root, value) - process_dataset(data_dir, label_dir, output_dir, max_images=args.max_images, + process_dataset(data_dir, label_dir, out_dir, max_images=args.max_images, post_process=post_process, max_z_spacing=MAX_Z_SPACING_MM, allowed_levels=LUMBAR_LEVELS, min_levels=MIN_LUMBAR_LEVELS, metadata_cache=metadata_db) diff --git a/xfr_reprocess_ap.py b/xfr_reprocess_ap.py new file mode 100644 index 0000000..a57fd48 --- /dev/null +++ b/xfr_reprocess_ap.py @@ -0,0 +1,128 @@ +#!/home/xfr/.conda/envs/cbt/bin/python +"""掃 standardize 輸出目錄中「前後（AP）方向翻轉」的個案（prone 伏位掃描， +前側落在 y 小側，例：colon 0003），供重跑修正用。 + +判定：對每個 volume 的 L1~L6 輸出遮罩（優 _binary_sdf、次 _binary_nn、 +再 _binary）個別跑 orientation.anterior_y_side（中線帶椎管兩側質量比較， +見該函式 docstring），多 level 票決： + flip : y_min 票 > y_max 票（前後翻轉，需重跑） + ok : y_max 票 > y_min 票（方向正常） + mixed : 平手（需人工確認） + unknown : 全部無法判定（無明顯前後質量差，如鏡稱面異常案例） + +Usage: + python xfr_reprocess_ap.py # 只回報 + python xfr_reprocess_ap.py --fix # 另刪 flip volume 的輸出 + # 資料夾 + progress.json + # 條目，之後重跑 + # xfr_preprocess.py 即可 +""" +import argparse +import json +import os +import shutil +import sys +import time + +import SimpleITK as sitk + +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) +from imaging.orientation import anterior_y_side + +LEVELS = ('L1', 'L2', 'L3', 'L4', 'L5', 'L6') +MASK_SUFFIXES = ('_binary_sdf.nii.gz', '_binary_nn.nii.gz', '_binary.nii.gz') + + +def volume_decision(vol_dir): + """回傳 (decision, per_level dict)。decision ∈ flip/ok/mixed/unknown/nomask。""" + per = {} + for lvl in LEVELS: + for suf in MASK_SUFFIXES: + p = os.path.join(vol_dir, f'{lvl}{suf}') + if os.path.exists(p): + m = sitk.GetArrayFromImage(sitk.ReadImage(p, sitk.sitkUInt8)) + per[lvl] = anterior_y_side(m) + break + if not per: + return 'nomask', per + votes = [v for v in per.values() if v is not None] + if not votes: + return 'unknown', per + n_min = votes.count('y_min') + n_max = votes.count('y_max') + if n_min > n_max: + return 'flip', per + if n_max > n_min: + return 'ok', per + return 'mixed', per + + +def main(): + parser = argparse.ArgumentParser( + description='Find AP-flipped (prone) volumes in a standardized output dir.') + parser.add_argument('output_dir') + parser.add_argument('--fix', action='store_true', + help='Also delete flipped volumes\' output dirs and ' + 'their progress.json entries') + args = parser.parse_args() + + outdir = args.output_dir + if not os.path.isdir(outdir): + print(f'not a directory: {outdir}') + return + + vols = sorted(d for d in os.listdir(outdir) + if os.path.isdir(os.path.join(outdir, d))) + flip, mixed, unknown, ok, nomask = [], [], [], [], [] + t0 = time.time() + for i, vol in enumerate(vols, 1): + dec, per = volume_decision(os.path.join(outdir, vol)) + tag = {'flip': 'FLIP', 'ok': 'ok ', 'mixed': 'MIXED', + 'unknown': '?!?', 'nomask': '- '}[dec] + detail = ' '.join(f'{k}={v}' for k, v in per.items()) + print(f'[{i}/{len(vols)}] {tag} {vol} {detail}') + {'flip': flip, 'mixed': mixed, 'unknown': unknown, + 'ok': ok, 'nomask': nomask}[dec].append(vol) + if (i % 25) == 0: + print(f' ... {i}/{len(vols)} ({(time.time()-t0)/60:.1f} min)', flush=True) + + print(f'\n=== Summary: {len(vols)} volumes ===') + print(f' ok (normal) : {len(ok)}') + print(f' FLIP (AP-flipped) : {len(flip)}') + for v in flip: + print(f' - {v}') + print(f' mixed (need check) : {len(mixed)}') + for v in mixed: + print(f' - {v}') + print(f' unknown (no vote) : {len(unknown)}') + for v in unknown: + print(f' - {v}') + print(f' no level mask : {len(nomask)}') + + if args.fix and flip: + # 1) progress.json：刪掉 flip 的條目（備份） + prog_path = os.path.join(outdir, 'progress.json') + if os.path.exists(prog_path): + with open(prog_path) as f: + prog = json.load(f) + removed = [v for v in flip if prog.pop(v, None) is not None] + bak = f'{prog_path}.apfix-{time.strftime("%Y%m%d_%H%M%S")}' + shutil.copyfile(prog_path, bak) + with open(prog_path, 'w') as f: + json.dump(prog, f, indent=2) + print(f'\nprogress.json: removed {len(removed)} entr(y/ies) ' + f'[{", ".join(v.split(".")[-1] or v for v in removed)}]; ' + f'backup {bak}') + # 2) 刪輸出資料夾 + for v in flip: + shutil.rmtree(os.path.join(outdir, v)) + print(f'removed {os.path.join(outdir, v)}') + print('\nNext: rerun `python xfr_preprocess.py` — only the removed ' + 'volumes will be reprocessed (with the AP flip applied).') + elif not args.fix and flip: + print('\nRerun with --fix to delete the flipped outputs and progress ' + 'entries, then run `python xfr_preprocess.py`.') + + +if __name__ == '__main__': + main() \ No newline at end of file