From 0a928d8f8e8c9b804734e47bcdeda63f446b8e08 Mon Sep 17 00:00:00 2001 From: Xiao Furen Date: Wed, 9 Sep 2026 19:08:52 +0800 Subject: [PATCH] refactor(imaging): update coordinate transformation logic to use transform.json Update the coordinate transformation pipeline to prioritize the affine transformation chain recorded in `transform.json` (original_to_source) instead of relying on re-calculated geometric parameters. The previous method relied on re-estimating rotation, center, and bounding boxes from the bone mask, which led to inaccuracies in screw pillar mapping (e.g., dropping from 99% to 30-94% in-label accuracy). The new approach uses the precise inverse affine mapping `o = M^-1 (s - t)` from the transformation metadata. A fallback mechanism is maintained for legacy volumes lacking `transform.json`, which continues to use the re-calculation method. - Implement `load_transform` and `original_to_source` integration - Update documentation to reflect the new primary/fallback coordinate chains - Improve precision of screw parameter mapping to native index space --- xfr_cbt_native.py | 98 +++++++++++++++++++++++++++++++++-------------- 1 file changed, 69 insertions(+), 29 deletions(-) diff --git a/xfr_cbt_native.py b/xfr_cbt_native.py index d28d472..e9076ba 100644 --- a/xfr_cbt_native.py +++ b/xfr_cbt_native.py @@ -8,27 +8,34 @@ label 值：L1L=1 L1R=2 L2L=3 L2R=4 L3L=5 L3R=6 L4L=7 L4R=8 L5L=9 L5R=10 （0 = 背景）。 -座標鏈（純 index 空間；不依賴 standardized 輸出的物理 header —— -standardize_affine 的 origin 處理不可靠，2026-09-08 已驗證）： - rotated disk r（(x,y,z) index） - -> template 0.5mm: t = R^T (r + fstart - c) + c [行向量: (r+fstart-c) @ R + c] - -> 記憶體 0.5mm 全域: g0 = bbox2s + (wx-1-tx, wy-1-ty, tz) - -> ap_flip 時: g0y = N05y - 1 - g0y - -> native index: rint(g0 * 0.5 / sn) +座標鏈首選：直接用預處理寫檔時記錄在 /transform.json 的正向鏈 +（imaging.transforms.original_to_source：boxes / std_flip_axes / ap_flip / +rotated R、center、start / 原 CT 幾何；純 index 空間、不依賴 standardized +輸出的物理 header）。該鏈為仿射（original index <-> rotated disk index）， +螺絲點經仿射逆向 o = M^-1 (s - t) 精確映射——不重新估算任何平面 / 幾何、 +不重取樣 native label。 -其中每 level 的 R/c 由 template 骨頭 mask（smd_resampled < 0.5，與 -_write_rotated_level 的輸入同定義；退化時退回 _binary_nn）重算； -fstart = rotated 檔 origin 反映到 template index 的整數； -bbox2 = native label 線性重取樣 0.5mm（>0.5）的最大 26-連通區域 bbox， -與 seg_bone 同定義（template/roi/binary_sdf 都裁在這個 bbox2 上）。 +fallback（僅舊世代 volume 沒有 transform.json、或該 level 沒有 rotated +記錄時）：純 index 空間重算 R/center（template 骨頭 mask 的 +best_symmetry_plane / best_upper_endplate_plane）、fstart（rotated 檔 +origin 反映到 template index 的整數）、bbox2（native label 0.5mm 線性 +重取樣 >0.5 的最大 26-連通區域 bbox）： + rotated disk r -> t = R^T (r + fstart - c) + c + -> g0 = bbox2s + (wx-1-tx, wy-1-ty, tz)；ap_flip 時 g0y = N05y-1-g0y + -> native: rint(g0 * 0.5 / sn) +重算結果可能與生成時參數不完全一致：2026-09-09 驗證 liver_100 螺絲柱 +in-label 由 93-99%（transform.json）掉到 30-94%（重算）；新世代資料 +（transform.json 在）一律走主路徑。 螺絲參數化（rotated 系、0.5mm index）：pos = [z, y, x, az°, alt°, d mm, L mm]； 方向 d_v = (cos az sin alt, sin az sin alt, cos alt)；末端 = p0 + L/0.5 * d_v； 柱半徑 = d/2 mm。 -驗證（2026-09-08）：整顆骨頭 15/15 (volume, level) 100% 落在 native label -（/tmp/kilo/validate_final.py）；螺絲柱體全部點 20/20 level-side -in-label >= 96%（2-voxel 膨脹）（/tmp/kilo/validate_screws2.py）。 +驗證（2026-09-09，transform.json 路徑）：20260909_141448 run — +0001（10/10）、liver_100（10/10）、covid19（L1/L2 六側）螺絲柱體 +in-label >= 92%、2-voxel 膨脹 >= 98%。covid19 L3 兩側除外：該 volume 的 +L3 旋轉 bone mask 是 4.5k voxel / 5 slice 的退化區域（L4/L5 預處理缺檔、 +先前棘突切除），屬資料問題非映射問題（L3 bone 碎片本身逆向 in-label 97.8%）。 """ import argparse @@ -45,7 +52,8 @@ from config.constant import LABEL_MAP from imaging.resample import resample_img from imaging.segmentation import _largest_cc_bbox from imaging.orientation import best_symmetry_plane, best_upper_endplate_plane -from imaging.transforms import level_file_path +from imaging.transforms import (level_file_path, load_transform, + original_to_source, transform_path) from visualization.res_bone_figure import compute_normalizing_rotation standardized_dir = '/mnt/1248/open2/cyrou/CBT/Seg/Resample/standardized-xfr-3' @@ -176,15 +184,32 @@ def write_volume_cbt(volume_id, run_id, output_root=Output_dir, date=None): ct_path, lb_path, ap_flip = find_native_paths(volume_id) ct = sitk.ReadImage(ct_path) - lb_img = sitk.ReadImage(lb_path) - lb_arr = sitk.GetArrayFromImage(lb_img) - Nn = np.array(ct.GetSize(), float) + n_xyz = np.array(ct.GetSize(), int) # (x,y,z) + n_zyx = n_xyz[::-1].copy() # 輸出 (z,y,x) sn = np.array(ct.GetSpacing(), float) - N05 = np.maximum(1, np.ceil(Nn * sn / 0.5 - 1e-6).astype(int)) + N05 = np.maximum(1, np.ceil(n_xyz.astype(float) * sn / 0.5 - 1e-6).astype(int)) - out_arr = np.zeros(lb_arr.shape, np.uint8) - geom_cache = {} + # 主路徑：transform.json 記錄的正向鏈（仿射），逆向 o = M^-1 (s - t) + meta = None + try: + meta = load_transform(os.path.join(standardized_dir, volume_id)) + except FileNotFoundError: + logger.warning(f'{volume_id}: no transform.json ' + f'({transform_path(os.path.join(standardized_dir, volume_id))}); ' + f'全部 level 退回重算路徑（近似）') + if meta is not None and not np.array_equal( + np.asarray(meta['original']['size'], int), n_xyz): + logger.warning(f'{volume_id}: transform.json original size ' + f'{meta["original"]["size"]} != CT {list(n_xyz)}; ' + f'忽略 transform.json，退回重算路徑') + meta = None + + out_arr = np.zeros(n_zyx, np.uint8) + lb_img = lb_arr = None # fallback 才需要讀 native label + geom_cache = {} # fallback：每 level 的 (R, c, fstart, bbox2s, wx, wy) + aff_cache = {} # 主路徑：每 level 的 (M^-1, t) n_screws, skipped = 0, [] + n_tf = n_fb = 0 for li, level in enumerate(LEVELS): for side in ('L', 'R'): jp = os.path.join(side_vol_dir, f'{level}_{side}.json') @@ -192,15 +217,29 @@ def write_volume_cbt(volume_id, run_id, output_root=Output_dir, date=None): continue pos = json.load(open(jp))['position'] try: - if level not in geom_cache: - geom_cache[level] = level_geometry(volume_id, level, lb_img, lb_arr) - geom = geom_cache[level] cyl = screw_voxel_xyz(pos) - nat = rotated_to_native(cyl, geom, N05, sn, ap_flip) + lv = meta['levels'].get(level) if meta is not None else None + if lv is not None and 'rotated' in lv: + if level not in aff_cache: + M, t = original_to_source(meta, level, 'rotated') + aff_cache[level] = (np.linalg.inv(M), t) + Mi, t = aff_cache[level] + nat = (cyl - t) @ Mi.T # (N,3) 原 CT 連續 index + n_tf += 1 + else: + if lb_arr is None: + lb_img = sitk.ReadImage(lb_path) + lb_arr = sitk.GetArrayFromImage(lb_img) + if level not in geom_cache: + geom_cache[level] = level_geometry(volume_id, level, lb_img, lb_arr) + nat = rotated_to_native(cyl, geom_cache[level], N05, sn, ap_flip) + n_fb += 1 ni = np.rint(nat).astype(int) - valid = (ni >= 0).all(1) & (ni < np.array(lb_arr.shape[::-1])).all(1) + valid = (ni >= 0).all(1) & (ni < n_xyz).all(1) idx = ni[valid] - val = li * 2 + (1 if side == 'R' else 0) + if idx.shape[0] == 0: + raise ValueError('all screw voxels out of original CT bounds') + val = li * 2 + 1 + (1 if side == 'R' else 0) out_arr[idx[:, 2], idx[:, 1], idx[:, 0]] = val n_screws += 1 except Exception as e: @@ -217,7 +256,8 @@ def write_volume_cbt(volume_id, run_id, output_root=Output_dir, date=None): sitk.WriteImage(out_img, out_path) if skipped: logger.warning(f'{volume_id}: failed sides: {", ".join(skipped)}') - logger.info(f'{volume_id}: {n_screws}/10 screws -> {out_path} (ap_flip={ap_flip})') + logger.info(f'{volume_id}: {n_screws}/10 screws -> {out_path} ' + f'(transform.json={n_tf}, re-est={n_fb}, ap_flip={ap_flip})') return out_path, n_screws