From d167c1f7c777c6fd8b24d6a7a7f2704ebcb72eac Mon Sep 17 00:00:00 2001 From: Xiao Furen Date: Wed, 9 Sep 2026 13:39:47 +0800 Subject: [PATCH] feat(imaging): implement coordinate transformation pipeline and directory restructuring Introduce a robust coordinate transformation system to manage the relationship between original CT space and rotated/standardized segmentation spaces. This includes a new directory hierarchy to separate unrotated crops from rotated outputs and utility functions for geometric mapping. Key changes: - Implement `imaging/transforms.py` to handle bounding box metadata, affine standardization, and coordinate mapping between spaces. - Restructure dataset output: unrotated segmentation files (binary, SDF, ROI, etc.) are now stored in a `/crop/` subdirectory to distinguish them from `/rotated/` aligned versions. - Add `level_file_path` utility to abstract file discovery across legacy (top-level) and new (crop-based) directory structures. - Enhance `seg_bone` to capture and export bounding box metadata (`bbox2`, `nn_bbox`, `bbox_orig`) into `transform.json`. - Implement `xfr_cbt_native.py` for mapping screw positions back to original CT space. - Update preprocessing and visualization scripts to support the new directory layout and transformation metadata. - Improve TinyDB metadata migration logic to prevent accidental corruption of existing database structures. --- imaging/affine.py | 5 + imaging/preprocessing.py | 78 +++- imaging/segmentation.py | 38 +- imaging/transforms.py | 304 +++++++++++++++ visualization/res_bone_figure.py | 5 +- xfr_cbt_native.py | 259 +++++++++++++ xfr_check_spinous.py | 6 +- xfr_debug.py | 16 +- xfr_inverse_transform.py | 507 +++++++++++++++++++++++++ xfr_migrate_crop.py | 113 ++++++ xfr_orig_labels.py | 628 +++++++++++++++++++++++++++++++ xfr_plot_level.py | 8 +- xfr_preprocess.py | 78 +++- xfr_reprocess_ap.py | 8 +- 14 files changed, 2006 insertions(+), 47 deletions(-) create mode 100644 imaging/transforms.py create mode 100644 xfr_cbt_native.py create mode 100644 xfr_inverse_transform.py create mode 100644 xfr_migrate_crop.py create mode 100644 xfr_orig_labels.py diff --git a/imaging/affine.py b/imaging/affine.py index c6291df..d6cae16 100644 --- a/imaging/affine.py +++ b/imaging/affine.py @@ -3,6 +3,10 @@ import numpy as np import nibabel as nib def standardize_affine(file_path, output_dir): + """翻轉 affine 對角為負的軸（資料 flip + 平移修正）後重寫同目錄同名檔。 + + 回傳：實際翻轉的軸 index list（未翻為 []）。呼叫端（transform 記錄） + 需要這份 ground truth：反 index 鏈路用「N-1-i」鏡射補回翻軸。""" img = nib.load(file_path) data = img.get_fdata() @@ -39,4 +43,5 @@ def standardize_affine(file_path, output_dir): os.makedirs(os.path.dirname(output_path), exist_ok=True) nib.save(standardized_img, output_path) + return flip_axes diff --git a/imaging/preprocessing.py b/imaging/preprocessing.py index 97a4499..964a373 100644 --- a/imaging/preprocessing.py +++ b/imaging/preprocessing.py @@ -9,6 +9,8 @@ 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 +from imaging.transforms import (build_volume_meta, img_geom, margined_box, + save_transform, std_flip_axes_for_direction) def flip_y_sitk(img): @@ -33,10 +35,13 @@ def save_progress(progress, PROGRESS_FILE): with open(PROGRESS_FILE, "w") as f: json.dump(progress, f, indent=2) -def process_single_image(image_path, label_path, output_dir_base=None, max_z_spacing=None, allowed_levels=None, min_levels=None, metadata_cache=None): +def process_single_image(image_path, label_path, output_dir_base=None, max_z_spacing=None, allowed_levels=None, min_levels=None, metadata_cache=None, crop_subdir=False): """metadata_cache：可選，需提供 get(name) -> dict|None 與 put(name, dict)。dict 可含 spacing=[x,y,z]、labels=[label id]。 - 兩者都在 db 裡時整支跳過判定不需讀影像 / label 檔。""" + 兩者都在 db 裡時整支跳過判定不需讀影像 / label 檔。 + crop_subdir=True：各 level 的未旋轉輸出（_binary / _smd / + _smd_resampled / _binary_sdf / _binary_nn / _roi）改写到 + /crop/ 子資料夾（False＝舊佈局、/ 頂層）。""" file_name = os.path.basename(image_path) name = file_name.replace(".nii.gz", "") @@ -191,6 +196,13 @@ def process_single_image(image_path, label_path, output_dir_base=None, max_z_spa name = file_name.replace(".nii.gz", "") output_dir = os.path.join(output_dir_base, name) os.makedirs(output_dir, exist_ok=True) + # 各 level 的未旋轉輸出寫入位置：crop_subdir=True 時在 /crop/ + # 子資料夾（旋轉版在 /rotated/，xfr_preprocess 產出） + if crop_subdir: + seg_dir = os.path.join(output_dir, 'crop') + os.makedirs(seg_dir, exist_ok=True) + else: + seg_dir = output_dir # 存現有 label 到 txt txt_path = os.path.join(output_dir, f"{name}_labels.txt") @@ -202,6 +214,8 @@ def process_single_image(image_path, label_path, output_dir_base=None, max_z_spa # 或 seg_bone 報錯）只跳過該 label，不中斷整檔處理 processed = [] skipped = [] + level_entries = {} # level name -> transform.json 的 level 條目 + std_flips = None # standardize_affine 實際翻的軸（整卷一致；記錄用） for n in existing_labels: if n not in LABEL_MAP: print(f"Label {n} not found in label_map, skipping this label (file continues).") @@ -213,7 +227,7 @@ def process_single_image(image_path, label_path, output_dir_base=None, max_z_spa skipped.append(n) continue try: - res = seg_bone(n, name, resampled_sitk_img, resampled_sitk_lbl, output_dir, + res = seg_bone(n, name, resampled_sitk_img, resampled_sitk_lbl, seg_dir, label_map=LABEL_MAP, original_label=label) if res is None: print(f"Label {n}: empty after largest-CC extraction, skipping this label.") @@ -221,24 +235,73 @@ def process_single_image(image_path, label_path, output_dir_base=None, max_z_spa continue (roi_path, binary_path, roi2_path, cortical_path, binary_nn_path, binary_linear_path, smd_path, resampled_path, binary_sdf_path, - binary_erode_path) = res + binary_erode_path, seg_meta) = res for path in [roi_path, binary_path, roi2_path, cortical_path, binary_nn_path, binary_linear_path, smd_path, resampled_path, binary_sdf_path, binary_erode_path]: if path is not None: - standardize_affine(path, output_dir) + flips = standardize_affine(path, seg_dir) + if flips: + if std_flips is None: + std_flips = sorted(flips) + elif sorted(flips) != std_flips: + print(f"WARNING: standardize_affine 翻軸不一致 " + f"({name}): {std_flips} vs {sorted(flips)}") + # 該 level 的裁切 box（x,y,z 序；imaging/transforms.py 定義） + boxes = {} + if seg_meta['has_smd']: + boxes['smd_resampled'] = seg_meta['bbox2'] + boxes['binary_sdf'] = seg_meta['bbox2'] + boxes['roi'] = seg_meta['bbox2'] + else: + boxes['roi'] = seg_meta['bbox2'] + if seg_meta['nn_bbox'] is not None: + boxes['binary_nn'] = seg_meta['nn_bbox'] + if seg_meta['bbox_orig'] is not None: + boxes['binary'] = seg_meta['bbox_orig'] + boxes['smd'] = margined_box(seg_meta['bbox_orig'], + image.GetSize(), margin=4) + level_entries[LABEL_MAP[n]] = { + 'label': int(n), + 'std_flip_axes': list(std_flips) if std_flips is not None else [], + 'boxes': boxes, + } processed.append(n) except RuntimeError as e: print(f"Label {n} could not be processed, skipping. Error: {e}") skipped.append(n) + # transform.json：每 level 記錄 原始 CT <-> 標準化 grid 的完整座標鏈 + # （box / flip 軸 / ap_flip / 原 CT 幾何），rotated/ 的 R/center/start 由 + # _write_rotated_level 補寫（xfr_preprocess）；供分割後 mask 映回原始座標 + if level_entries: + try: + tmeta = build_volume_meta(name, image_path, img_geom(image), + img_geom(resampled_sitk_img), ap_flip, + level_entries) + save_transform(output_dir, tmeta) + print(f"Transform metadata: {os.path.join(output_dir, 'transform.json')} " + f"({len(level_entries)} level(s), std_flips={tmeta['levels'][list(level_entries)[0]]['std_flip_axes']})") + try: + expected = std_flip_axes_for_direction(image.GetDirection()) + if std_flips is not None and list(std_flips) != expected: + print(f"WARNING: {name} std_flips {std_flips} != 預期 " + f"{expected}（依原 CT direction 推斷）") + elif std_flips is None and expected: + print(f"WARNING: {name} 預期 std_flips {expected} 但未翻任何軸") + except ValueError: + print(f"WARNING: {name} 原 CT direction 非對角 ±1；" + f"transform.json 的反向映射將無法使用該卷") + except Exception as e: + print(f"WARNING: 寫 transform.json 失敗（{name}）：{e}") + return { "processed_labels": processed, "missing_labels": skipped } -def process_dataset(image_dir, label_dir, output_dir, labels_to_process=None, max_images=None, post_process=None, max_z_spacing=None, allowed_levels=None, min_levels=None, metadata_cache=None): +def process_dataset(image_dir, label_dir, output_dir, labels_to_process=None, max_images=None, post_process=None, max_z_spacing=None, allowed_levels=None, min_levels=None, metadata_cache=None, crop_subdir=False): image_files = sorted(glob.glob(os.path.join(image_dir, "*.nii.gz"))) total_files = len(image_files) print(f"Total files: {total_files}") @@ -286,7 +349,8 @@ def process_dataset(image_dir, label_dir, output_dir, labels_to_process=None, ma max_z_spacing=max_z_spacing, allowed_levels=allowed_levels, min_levels=min_levels, - metadata_cache=metadata_cache) + metadata_cache=metadata_cache, + crop_subdir=crop_subdir) # print(result) # exit() except Exception as e: diff --git a/imaging/segmentation.py b/imaging/segmentation.py index 936ca54..3f337d1 100644 --- a/imaging/segmentation.py +++ b/imaging/segmentation.py @@ -48,6 +48,7 @@ def seg_bone(n, name, resampled_sitk_img, resampled_sitk_lbl, output_base=None, # 1. 提取標籤 n 的二值遮罩 (將標籤 n 設為 1，其餘為 0)，最大連通區域 smd_path = resampled_path = binary_sdf_path = binary_erode_path = None binary_linear_path = binary_nn_path = None + nn_bbox = bbox_orig = None # 裁切 box（x,y,z 序），供 transform.json 記錄 if original_label is not None: bin_orig = sitk.BinaryThreshold(original_label, n, n, 1, 0) cc_orig = _largest_cc_bbox(bin_orig) @@ -142,20 +143,37 @@ def seg_bone(n, name, resampled_sitk_img, resampled_sitk_lbl, output_base=None, roi_path = os.path.join(output_base, f"{label_name}_roi.nii.gz") sitk.WriteImage(roi, roi_path) + # 第 11 項：該 level 的裁切 box（x,y,z 序 [x0,y0,z0,xs,ys,zs]）與來源資訊， + # 供 transform.json（imaging/transforms.py）記錄原始 <-> 標準化座標鏈： + # bbox2 : 0.5mm 最大連通 bbox（_smd_resampled / _binary_sdf / _roi 共用） + # nn_bbox : _binary_nn 自己的 bbox（_smd_resampled 路徑才產生） + # bbox_orig : 原解析度最大連通 bbox（_binary；_smd = 其 +4 margin） return roi_path, binary_path, None, None, binary_nn_path, \ - binary_linear_path, smd_path, resampled_path, binary_sdf_path, binary_erode_path + binary_linear_path, smd_path, resampled_path, binary_sdf_path, \ + binary_erode_path, { + 'label': int(n), + 'level': label_name, + 'has_smd': smd_path is not None, + 'bbox2': [int(v) for v in bbox2], + 'nn_bbox': [int(v) for v in nn_bbox] if nn_bbox is not None else None, + 'bbox_orig': [int(v) for v in bbox_orig] if bbox_orig is not None else None, + } """ Dataset/ └── standardized/ └── subject001/ - ├── L1_binary.nii.gz # 原解析度【未插值】遮罩（最大连通區域、裁物件 bbox） - ├── L1_smd.nii.gz # SignedMaurerDistanceMap（內負/外正，原始 index 單位； - │ # bbox 外扩 4 voxel 背景輪，供重取樣插值用） - ├── L1_smd_resampled.nii.gz # _smd 經线性插值重取樣到 0.5mm（浮點，裁 0.5mm bbox） - ├── L1_binary_sdf.nii.gz # _smd_resampled 於 0.5 閾值 -> 0.5mm 平滑 mask（0/1，裁同 bbox） - ├── L1_binary_nn.nii.gz # 最近邻版 0/1（對比用，各自 bbox） - ├── L1_roi.nii.gz - ├── L2_binary.nii.gz - ... + ├── crop/ # 未旋转裁切片（旋转版见 /rotated/） + │ ├── L1_binary.nii.gz # 原解析度【未插值】遮罩（最大连通区域、裁物件 bbox） + │ ├── L1_smd.nii.gz # SignedMaurerDistanceMap（内负/外正，原始 index 单位； + │ │ # bbox 外扩 4 voxel 背景轮，供重取样插值用） + │ ├── L1_smd_resampled.nii.gz # _smd 经线性插值重取样到 0.5mm（浮点，裁 0.5mm bbox） + │ ├── L1_binary_sdf.nii.gz # _smd_resampled 于 0.5 阈值 -> 0.5mm 平滑 mask（0/1，裁同 bbox） + │ ├── L1_binary_nn.nii.gz # 最近邻版 0/1（对比用，各自 bbox） + │ ├── L1_roi.nii.gz + │ ├── L2_binary.nii.gz + │ ... + │ └── L1_planes.png # 未旋转 planes 图（xfr_preprocess 后处理） + ├── transform.json + └── rotated/ # 旋转对齐版（xfr_preprocess 后处理） """ \ No newline at end of file diff --git a/imaging/transforms.py b/imaging/transforms.py new file mode 100644 index 0000000..0fa1ef9 --- /dev/null +++ b/imaging/transforms.py @@ -0,0 +1,304 @@ +"""Per-level 座標變換記錄（transform.json）與反向映射：把標準化 +（0.5mm / rotated）grid 上的 segmentation mask 映回 data_root 原始 CT grid。 + +座標鏈（純 index 空間、(x,y,z) 向量慣例、array = (z,y,x)；不依賴 standardized +輸出的物理 header —— standardize_affine 的 origin 處理不可靠，2026-09-08 已 +驗證，xfr_cbt_native / xfr_orig_labels 同結論）： + + original index o（data_root CT） + -> 真 0.5mm full grid: t = o * (so / 0.5) + -> 工作系（ap_flip 後）: ap_flip 時 g_y = N05_y - 1 - t_y + -> 未旋轉檔 pre-standardize: p = g - box（box = 該檔在 full grid 的裁切） + -> disk（standardize 後）: std_flip_axes 的軸: disk_i = N_i - 1 - p_i + -> rotated disk: a = R (p_tpl - center) + center - start + +transform.json 記錄每一步所需的參數（寫檔時即為 ground truth，反向映射不需 +重新估算任何平面 / 幾何）： + + { + "version": 1, + "name": "", + "original": {"path", "size", "spacing", "direction", "origin"}, + "resampled05": {"size", "spacing", "direction", "origin"}, // 0.5mm full grid + "ap_flip": false, + "levels": { + "L1": { + "label": 20, + "std_flip_axes": [0, 1], // standardize_affine 實際翻的軸 + "boxes": { // [x0,y0,z0,xs,ys,zs]（x,y,z 序）； + "smd_resampled": [...], // 0.5mm 檔: 0.5mm full grid（ap_flip 後） + "binary_sdf": [...], // 原解析度檔 (binary/smd): 原 index + "roi": [...], // （ap_flip 後） + "binary_nn": [...], + "binary": [...], + "smd": [...] + }, + "rotated": { // _write_rotated_level 補寫 + "template": "smd_resampled", // 或 "binary_nn"（fallback） + "R": [[..]], // 作用於 (x,y,z)；forward dest = R(src-c)+c + "center": [cx, cy, cz], // 未旋轉檔 disk（standardize 後）index 系 + "start": [sx, sy, sz], // rotated 檔 origin（模板 index 系） + "size": [nx, ny, nz] + } + } + } + } + +original_to_source 回傳 (M (3,3), t (3,))：source disk index = M @ original +index + t（連續座標，x,y,z 系）。mask_to_original 對每個 original voxel 以 +order=0（NN）在 source mask 上採樣，得原始 grid 的 mask。 +""" + +import json +import os + +import numpy as np +import SimpleITK as sitk +from scipy.ndimage import map_coordinates + +import logging + +logger = logging.getLogger('imaging.transforms') + +TRANSFORM_FILENAME = 'transform.json' + +# 0.5mm 標準化 grid 的 source（boxes 記在 0.5mm full grid、ap_flip 後） +SOURCES_05 = ('smd_resampled', 'binary_sdf', 'roi', 'binary_nn') +# 原解析度 grid 的 source（boxes 記在原始 index、ap_flip 後） +SOURCES_ORIG = ('binary', 'smd') +SOURCES = SOURCES_05 + SOURCES_ORIG + ('rotated',) + + +def transform_path(vol_dir): + return os.path.join(vol_dir, TRANSFORM_FILENAME) + + +def level_file_path(vol_dir, level, source): + """_.nii.gz 的位置（未旋轉、裁切 level bbox 的檔）： + 新世代（xfr-3 起）落在 /crop/ 子資料夾（旋轉版在 /rotated/）， + 舊世代（xfr-2 等）落在 / 頂層。回傳先存在者（crop/ 優先）； + 兩者皆無時回 crop/ 路徑（呼叫端以 os.path.exists 判定）.""" + p_crop = os.path.join(vol_dir, 'crop', f'{level}_{source}.nii.gz') + if os.path.exists(p_crop): + return p_crop + p_top = os.path.join(vol_dir, f'{level}_{source}.nii.gz') + if os.path.exists(p_top): + return p_top + return p_crop + + +def load_transform(vol_dir): + p = transform_path(vol_dir) + if not os.path.exists(p): + raise FileNotFoundError( + f'{p} 不存在（需先跑 xfr_preprocess pipeline，' + f'或 xfr_inverse_transform.py --rebuild）') + with open(p) as f: + return json.load(f) + + +def save_transform(vol_dir, meta): + p = transform_path(vol_dir) + tmp = f'{p}.tmp' + with open(tmp, 'w') as f: + json.dump(meta, f, indent=2, ensure_ascii=False) + os.replace(tmp, p) + + +def merge_rotated_into_transform(vol_dir, level, sec): + """讀入既有 transform.json，補 / 覆寫 levels[level]['rotated'] 後寫回 + （_write_rotated_level 用；base 部分由 process_single_image 寫入）。""" + meta = load_transform(vol_dir) + lv = meta['levels'].setdefault(level, {}) + lv['rotated'] = sec + save_transform(vol_dir, meta) + + +def img_geom(img): + """sitk img -> {"size"(x,y,z), "spacing", "direction"(row-major 9), "origin"}。""" + return { + 'size': [int(v) for v in img.GetSize()], + 'spacing': [float(v) for v in img.GetSpacing()], + 'direction': [float(v) for v in img.GetDirection()], + 'origin': [float(v) for v in img.GetOrigin()], + } + + +def require_diagonal(direction, what='direction'): + d = np.asarray(direction, dtype=float).reshape(3, 3) + diag = np.diag(d) + if not (np.allclose(d, np.diag(diag), atol=1e-6) + and np.allclose(np.abs(diag), 1.0, atol=1e-6)): + raise ValueError(f'{what} 非對角 ±1（本座標鏈只支援對角 direction）: ' + f'{d.tolist()}') + return diag + + +def std_flip_axes_for_direction(direction): + """standardize_affine 會翻的軸：儲存 NIfTI affine（RAS）對角 < 0 者。 + RAS = LPS 的 x/y 取反：ras_diag = (-d_x, -d_y, +d_z)。 + direction 非對角 ±1 時 raise ValueError。""" + d = require_diagonal(direction, 'direction') + ras = np.array([-d[0], -d[1], d[2]]) + return [int(i) for i in range(3) if ras[i] < 0] + + +def resampled05_geom_from_original(orig_geom): + """resample_img 的 0.5mm full grid 幾何（origin/direction 同原、尺寸 ceil）。""" + n = np.asarray(orig_geom['size'], dtype=float) + s = np.asarray(orig_geom['spacing'], dtype=float) + return { + 'size': [max(1, int(vv)) for vv in np.ceil(n * s / 0.5 - 1e-6)], + 'spacing': [0.5, 0.5, 0.5], + 'direction': [float(v) for v in orig_geom['direction']], + 'origin': [float(v) for v in orig_geom['origin']], + } + + +def margined_box(box, n_xyz, margin=4): + """seg_bone 的 _bbox_roi 同式：box [x0,y0,z0,xs,ys,zs] 對稱外擴 margin、 + clamp 到 n_xyz（x,y,z 序）。""" + n = np.asarray(n_xyz, dtype=float) + idx = [max(0, int(box[i]) - margin) for i in range(3)] + size = [min(int(n[i]) - idx[i], int(box[i + 3]) + 2 * margin) for i in range(3)] + return idx + size + + +def build_volume_meta(name, orig_path, orig_geom, r05_geom, ap_flip, level_entries): + return { + 'version': 1, + 'name': name, + 'original': {'path': orig_path, **orig_geom}, + 'resampled05': r05_geom, + 'ap_flip': bool(ap_flip), + 'levels': level_entries, + } + + +# ---------------------------------------------------------------- 座標鏈 +def _compose(f1, f2): + """f1(f2(x))：先 f2 後 f1。f = (M, t)：x -> M @ x + t。""" + M1, t1 = f1 + M2, t2 = f2 + return M1 @ M2, M1 @ t2 + t1 + + +def original_to_source(meta, level, source): + """原 CT index (x,y,z) -> source disk index (x,y,z) 的連續仿射 + (M (3,3), t (3,))：s = M @ o + t（ NN 前座標）。 + + source ∈ SOURCES。rotated 時 chain 末端加 R 旋轉 + start 平移 + （template = 未旋轉 0.5mm 檔，flip / box 用 template 的）。 + 所有參數來自 meta（transform.json 記錄值，非從 header 估算）。""" + if source not in SOURCES: + raise ValueError(f'unknown source {source!r}; 預期 {SOURCES}') + lv = meta['levels'].get(level) + if lv is None: + raise ValueError(f'level {level!r} 不在 transform.json ' + f'({sorted(meta["levels"])})') + if source == 'rotated' and 'rotated' not in lv: + raise ValueError(f'{level} 沒有 rotated 記錄（该 level 未走 _write_rotated_level）') + + so = np.asarray(meta['original']['spacing'], dtype=float) + n_orig = np.asarray(meta['original']['size'], dtype=float) + n05 = np.asarray(meta['resampled05']['size'], dtype=float) + ap_flip = bool(meta.get('ap_flip', False)) + + if source == 'rotated': + rot = lv['rotated'] + tpl = rot['template'] + Mrot = np.asarray(rot['R'], dtype=float) + c = np.asarray(rot['center'], dtype=float) + start = np.asarray(rot['start'], dtype=float) + else: + tpl = source + Mrot = c = start = None + if tpl not in lv['boxes']: + raise ValueError(f'{level} 缺 boxes[{tpl!r}]（transform.json 未完）') + box = np.asarray(lv['boxes'][tpl], dtype=float) + + res = '05' if tpl in SOURCES_05 else 'orig' + flips = [int(i) for i in lv.get('std_flip_axes', [])] + n_file = box[3:6] # 該檔尺寸（= flip 時的 N） + + f = (np.eye(3), np.zeros(3)) + # 1) original index -> 真 0.5mm full grid index（0.5mm source 才需要） + if res == '05': + f = _compose((np.diag(so / 0.5), np.zeros(3)), f) + # 2) -> 工作系（ap_flip 後的 0.5mm / 原解析度 full grid index） + if ap_flip: + extent = n05[1] if res == '05' else n_orig[1] + M = np.eye(3) + M[1, 1] = -1.0 + f = _compose((M, np.array([0.0, float(extent) - 1.0, 0.0])), f) + # 3) -> 未旋轉檔 pre-standardize index（裁切 box 的局部系） + f = _compose((np.eye(3), -box[:3]), f) + # 4) -> disk（standardize_affine 翻軸後的存檔 index） + for i in flips: + if i >= 3: + raise ValueError(f'std_flip_axes 含非法軸 {i}') + M = np.eye(3) + M[i, i] = -1.0 + t = np.zeros(3) + t[i] = float(n_file[i]) - 1.0 + f = _compose((M, t), f) + # 5) (rotated) 未旋轉 template disk -> rotated disk： + # forward: W = R^T (V - c) + c, V = a + start => a = R(W - c) + c - start + if Mrot is not None: + f = _compose((Mrot, c - Mrot @ c - start), f) + return f + + +def expected_source_size(meta, level, source): + """該 source 檔的 (x,y,z) 尺寸（rotated 用 rotated.size，其餘用 boxes）。""" + lv = meta['levels'][level] + if source == 'rotated': + return [int(v) for v in lv['rotated']['size']] + box = lv['boxes'][source] + return [int(box[i + 3]) for i in range(3)] + + +def mask_to_original(mask_img, meta, level, source, chunk=16): + """把 source grid 上的 mask（nifti 或 sitk image）NN 反向映射回原始 CT + grid。回傳 (z,y,x) ndarray（dtype = 輸入 dtype）。 + + chunk：每個 z 區塊的 slice 數（原始 CT z 可能上百；分塊控制峰值記憶體）。 + mask 的 size 與該 source 檔不一致時 raise（grid 不對應，映射必錯）。""" + M, t = original_to_source(meta, level, source) + arr = mask_img if isinstance(mask_img, np.ndarray) \ + else sitk.GetArrayFromImage(mask_img) + if isinstance(mask_img, sitk.Image): + actual = np.asarray(mask_img.GetSize()) + else: + actual = np.array(arr.shape[::-1]) + expect = np.asarray(expected_source_size(meta, level, source)) + if not np.array_equal(actual, expect): + raise ValueError( + f'mask 尺寸 {tuple(map(int, actual))} 與 {level}/{source} 記錄尺寸 ' + f'{tuple(map(int, expect))} 不符（grid 不對應，無法映射）') + if arr.size == 0: + raise ValueError('輸入 mask 為空') + + nx, ny, nz = (int(v) for v in meta['original']['size']) + out = np.zeros((nz, ny, nx), dtype=np.float64) + for z0 in range(0, nz, int(chunk)): + z1 = min(nz, z0 + int(chunk)) + oz, oy, ox = np.indices((z1 - z0, ny, nx)) + X = ox.astype(np.float64) + Y = oy.astype(np.float64) + Z = (oz + z0).astype(np.float64) + P = np.stack([X, Y, Z], axis=0) # (3, cz, ny, nx) (x,y,z) + S = np.tensordot(M, P, axes=([1], [0])) + t[:, np.newaxis, np.newaxis, + np.newaxis] + out[z0:z1, ...] = map_coordinates(arr, [S[2], S[1], S[0]], + order=0, cval=0.0, mode='constant') + return out.astype(arr.dtype, copy=False) + + +def save_as_original(arr, orig_img, out_path): + """arr (z,y,x) -> nifti，幾何取 orig_img（保留 arr 自身 pixel type）。""" + out = sitk.GetImageFromArray(arr) + out.SetSpacing(orig_img.GetSpacing()) + out.SetDirection(orig_img.GetDirection()) + out.SetOrigin(orig_img.GetOrigin()) + sitk.WriteImage(out, out_path) \ No newline at end of file diff --git a/visualization/res_bone_figure.py b/visualization/res_bone_figure.py index fe490a6..4e8ddfd 100644 --- a/visualization/res_bone_figure.py +++ b/visualization/res_bone_figure.py @@ -640,11 +640,12 @@ def render_bone_figure(volume_id, level, binary_path, cortical_path, sc_bone = ax.scatter(x_base, y_base, z_base, c=rgba_base, s=size_base, marker="o") sc_bone.set_zorder(5) if x_vb.size: - sc_vb = ax.scatter(x_vb, y_vb, z_vb, c=to_rgba("gold", 0.95), + # (1,4) 單行 2D 陣列：整組點共用同一色，避免 *c* 被視為數值映射 + sc_vb = ax.scatter(x_vb, y_vb, z_vb, c=np.array([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), + sc_sp = ax.scatter(x_sp, y_sp, z_sp, c=np.array([to_rgba("purple", 0.95)]), s=BONE_MARKER_SIZE, marker="o") sc_sp.set_zorder(6.5) if _EX is not None: diff --git a/xfr_cbt_native.py b/xfr_cbt_native.py new file mode 100644 index 0000000..d28d472 --- /dev/null +++ b/xfr_cbt_native.py @@ -0,0 +1,259 @@ +#!/home/xfr/.conda/envs/cbt/bin/python +"""把 PSO 找到的 CBT 螺絲位置（旋轉 0.5mm 標準系）映回原始 +（未旋轉、未重取樣）native CT 空間，整支 volume 各 level 的螺絲 +（L1-L5 x L/R，最多 10 支）存成單一 label 體積： + + Output_dir///cbt.nii.gz + +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) + +其中每 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 上）。 + +螺絲參數化（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）。 +""" + +import argparse +import json +import logging +import os +import sys + +import numpy as np +import SimpleITK as sitk +from tinydb import TinyDB, Query + +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 visualization.res_bone_figure import compute_normalizing_rotation + +standardized_dir = '/mnt/1248/open2/cyrou/CBT/Seg/Resample/standardized-xfr-3' +data_root = '/mnt/1220/Public/dataset/Spine/CTSpine1K/data/' +label_root = '/mnt/1220/Public/dataset/Spine/CTSpine1K/label/' +SUBDIRS = { + 'colon': 'conlon', + 'COVID-19': 'COVID-19', + 'HNSCC-3DCT-RT_neck': 'HNSCC-3DCT-RT_neck', + 'liver': 'Liver', +} +Output_dir = '/mnt/1248/open/cyrou/Output' + +_PROJ_DIR = os.path.dirname(os.path.abspath(__file__)) +LOG_DIR = os.path.join(_PROJ_DIR, 'logs') +SIDE_RESULT_DIR = os.path.join(LOG_DIR, 'side_results') +META_DB = os.path.join(_PROJ_DIR, 'xfr_image_metadata.json') + +LEVELS = ('L1', 'L2', 'L3', 'L4', 'L5') +LEVEL_LABEL_VAL = {v: int(k) for k, v in LABEL_MAP.items() if v in LEVELS} # {'L1': 20, ...} + +logger = logging.getLogger('xfr_cbt_native') + + +def find_native_paths(volume_id): + """回傳 (ct_path, label_path, ap_flip)；找不到 native 檔對時 raise。""" + ap_flip = False + try: + meta = TinyDB(META_DB, access_mode='r').table('images') \ + .get(Query().name == volume_id) + if meta is None: + logger.warning(f'{volume_id}: not in metadata db; ap_flip defaults to False') + else: + ap_flip = bool(meta.get('ap_flip', False)) + except Exception as e: + logger.warning(f'{volume_id}: metadata db unavailable ({e}); ap_flip defaults to False') + for sub, lab_sub in SUBDIRS.items(): + ct = f'{data_root}{sub}/{volume_id}.nii.gz' + lb = f'{label_root}{lab_sub}/{volume_id}_seg.nii.gz' + if os.path.isfile(ct) and os.path.isfile(lb): + return ct, lb, ap_flip + raise FileNotFoundError(f'{volume_id}: no native CT/label pair under {data_root}') + + +def level_geometry(volume_id, level, lb_img, lb_arr): + """每 level 的 (R, c, fstart, bbox2s, wx, wy)；與預處理/驗證腳本同定義。""" + vol_dir = os.path.join(standardized_dir, volume_id) + # 未旋轉檔：新世代在 crop/ 子資料夾、舊世代在頂層 + smd_img = sitk.ReadImage(level_file_path(vol_dir, level, 'smd_resampled')) + smd = sitk.GetArrayFromImage(smd_img).astype(np.float32) + if not (np.isfinite(smd).all() and smd.min() < 0 and smd.max() > 0): + b = (sitk.GetArrayFromImage( + sitk.ReadImage(level_file_path(vol_dir, level, 'binary_nn'))) > 0).astype(np.uint8) + else: + b = (smd < 0.5).astype(np.uint8) + sym = best_symmetry_plane(b) + symp = best_upper_endplate_plane(b) + R, c = compute_normalizing_rotation(b, sym, symp) + rot_img = sitk.ReadImage(f'{vol_dir}/rotated/{level}_binary_sdf.nii.gz') + fstart = np.round(np.array( + smd_img.TransformPhysicalPointToIndex(rot_img.GetOrigin()))).astype(float) + + lv = LEVEL_LABEL_VAL[level] + bin_img = sitk.GetImageFromArray((lb_arr == lv).astype(np.uint8)) + bin_img.CopyInformation(lb_img) + bin_lin = resample_img(sitk.Cast(bin_img, sitk.sitkFloat32)) + m_full = sitk.GetImageFromArray((sitk.GetArrayFromImage(bin_lin) > 0.5).astype(np.uint8)) + m_full.CopyInformation(bin_lin) + cc = _largest_cc_bbox(m_full) + if cc is None: + raise ValueError(f'{volume_id} {level}: empty native level mask (label {lv})') + _, bbox2 = cc + wx, wy = smd_img.GetSize()[0], smd_img.GetSize()[1] + return R, c, fstart, np.array(bbox2[:3], float), wx, wy + + +def rotated_to_native(r_xyz, geom, N05, sn, ap_flip): + """r_xyz: (N,3) rotated disk (x,y,z) index -> (N,3) native 連續 index。""" + R, c, fstart, bbox2s, wx, wy = geom + t = (np.asarray(r_xyz, float) + fstart - c) @ R + c + g0 = np.empty_like(t) + g0[:, 0] = bbox2s[0] + (wx - 1 - t[:, 0]) + g0[:, 1] = bbox2s[1] + (wy - 1 - t[:, 1]) + g0[:, 2] = bbox2s[2] + t[:, 2] + if ap_flip: + g0[:, 1] = N05[1] - 1 - g0[:, 1] + return g0 * (0.5 / sn) + + +def screw_voxel_xyz(pos): + """pos = [z, y, x, az°, alt°, d, L]（rotated 系 0.5mm index）。 + 回傳柱體內 (x,y,z) float voxel 座標（N,3）。""" + z, y, x, az, alt, d, L = (float(v) for v in pos[:7]) + azr, altr = np.radians(az), np.radians(alt) + ca, sa = np.cos(azr), np.sin(azr) + ct_, st_ = np.cos(altr), np.sin(altr) + dv = np.array([ca * st_, sa * st_, ct_]) + e1 = np.array([ca * ct_, sa * ct_, -st_]) + e2 = np.array([-sa, ca, 0.0]) + p0 = np.array([x, y, z]) + p1 = p0 + (L / 0.5) * dv + rad = d / 0.5 + 3.0 + lo = np.floor(np.minimum(p0, p1)) - rad + hi = np.ceil(np.maximum(p0, p1)) + rad + xs, ys, zs = np.meshgrid(np.arange(lo[0], hi[0] + 1), + np.arange(lo[1], hi[1] + 1), + np.arange(lo[2], hi[2] + 1), indexing='xy') + P = np.stack([xs.ravel(), ys.ravel(), zs.ravel()], 1) + dP = P - p0 + xr, yr, zr = dP @ e1, dP @ e2, dP @ dv + rrad = d # (d/2 mm 半徑) / 0.5mm 每 voxel = d 個 voxel + m = (xr ** 2 + yr ** 2 <= rrad ** 2) & (zr >= 0) & (zr <= L / 0.5) + return P[m] + + +def write_volume_cbt(volume_id, run_id, output_root=Output_dir, date=None): + """把一個 volume 的 side_results 螺絲全部映回 native 空間寫 cbt.nii.gz。 + 回傳 (path, n_screws)；該 run 無此 volume 時回傳 (None, 0)。""" + date = date or run_id[:8] + side_vol_dir = os.path.join(SIDE_RESULT_DIR, run_id, volume_id) + if not os.path.isdir(side_vol_dir): + logger.warning(f'{volume_id}: no side results under {side_vol_dir}; cbt skipped') + return None, 0 + if not any(os.path.isfile(os.path.join(side_vol_dir, f'{l}_{s}.json')) + for l in LEVELS for s in ('L', 'R')): + logger.warning(f'{volume_id}: no _.json in {side_vol_dir}; cbt skipped') + return None, 0 + + 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) + sn = np.array(ct.GetSpacing(), float) + N05 = np.maximum(1, np.ceil(Nn * sn / 0.5 - 1e-6).astype(int)) + + out_arr = np.zeros(lb_arr.shape, np.uint8) + geom_cache = {} + n_screws, skipped = 0, [] + for li, level in enumerate(LEVELS): + for side in ('L', 'R'): + jp = os.path.join(side_vol_dir, f'{level}_{side}.json') + if not os.path.isfile(jp): + 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) + ni = np.rint(nat).astype(int) + valid = (ni >= 0).all(1) & (ni < np.array(lb_arr.shape[::-1])).all(1) + idx = ni[valid] + val = li * 2 + (1 if side == 'R' else 0) + out_arr[idx[:, 2], idx[:, 1], idx[:, 0]] = val + n_screws += 1 + except Exception as e: + logger.error(f'{volume_id} {level}_{side}: {e}') + skipped.append(f'{level}_{side} ({e})') + + if n_screws == 0: + raise ValueError(f'{volume_id}: no screw could be written (skipped: {skipped or "none read"})') + out_dir = os.path.join(output_root, date, volume_id) + os.makedirs(out_dir, exist_ok=True) + out_img = sitk.GetImageFromArray(out_arr) + out_img.CopyInformation(ct) + out_path = os.path.join(out_dir, 'cbt.nii.gz') + 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})') + return out_path, n_screws + + +def main(): + parser = argparse.ArgumentParser( + description='Write ///cbt.nii.gz (screws mapped to native space) ' + 'from a side_results run') + parser.add_argument('run_id', help='e.g. 20260907_230652') + parser.add_argument('volumes', nargs='*', + help='volume ids (full or trailing digits); default: all in the run') + parser.add_argument('--date', default=None, help='output date dir (default: run_id[:8])') + args = parser.parse_args() + logging.basicConfig(level=logging.INFO, + format='%(asctime)s [%(levelname)s] %(name)s: %(message)s', + datefmt='%Y-%m-%d %H:%M:%S') + run_dir = os.path.join(SIDE_RESULT_DIR, args.run_id) + if not os.path.isdir(run_dir): + sys.exit(f'side results dir not found: {run_dir}') + run_vols = sorted(d for d in os.listdir(run_dir) if os.path.isdir(os.path.join(run_dir, d))) + if args.volumes: + vols = [v for v in run_vols + if v in args.volumes or v.rsplit('.', 1)[-1] in args.volumes] + if not vols: + sys.exit(f'none of {args.volumes} found in run {args.run_id}') + else: + vols = run_vols + ok, fail = 0, 0 + for vid in vols: + try: + write_volume_cbt(vid, args.run_id, date=args.date) + ok += 1 + except Exception as e: + fail += 1 + logger.error(f'{vid}: {e}') + print(f'done: {ok} written, {fail} failed / {len(vols)} volume(s)') + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/xfr_check_spinous.py b/xfr_check_spinous.py index 02b8775..878b584 100644 --- a/xfr_check_spinous.py +++ b/xfr_check_spinous.py @@ -42,6 +42,7 @@ import SimpleITK as sitk from imaging.orientation import (best_symmetry_plane, segment_spinous_process, best_upper_endplate_plane, segment_vertebral_body, diagnose_spinous_process) +from imaging.transforms import level_file_path standardized_dir = '/mnt/1248/open2/cyrou/CBT/Seg/Resample/standardized-xfr/' @@ -141,9 +142,10 @@ def full_check(m): # ===================================================================== def load_binary(volume_dir, level): # _binary.nii.gz 現為原解析度（未插值）；0.5mm 分析優 _binary_sdf（SDF 平滑遮罩） - p = os.path.join(volume_dir, f'{level}_binary_sdf.nii.gz') + # 未旋轉檔：新世代在 crop/ 子資料夾、舊世代在頂層 + p = level_file_path(volume_dir, level, 'binary_sdf') if not os.path.exists(p): - p = os.path.join(volume_dir, f'{level}_binary.nii.gz') + p = level_file_path(volume_dir, level, 'binary') if not os.path.exists(p): return None m = sitk.GetArrayFromImage(sitk.ReadImage(p, sitk.sitkUInt8)) > 0 diff --git a/xfr_debug.py b/xfr_debug.py index 9e04e4f..8999504 100644 --- a/xfr_debug.py +++ b/xfr_debug.py @@ -17,6 +17,9 @@ from core.cylinder import create_coordinate_grid from core.objective import set_global_context from core.optimizer import run_pso_torch, run_de_torch, run_nm_torch, run_pso_torch_xfr from imaging.orientation import azimuth_rotation, analyze_vertebral_tilt_contour +from imaging.transforms import level_file_path + +import xfr_cbt_native standardized_dir = '/mnt/1248/open2/cyrou/CBT/Seg/Resample/standardized-xfr-3/' @@ -165,9 +168,10 @@ def debug_orientation(volume_id, level): volume_dir = os.path.join(standardized_dir, volume_id) # _binary.nii.gz 現為原解析度；0.5mm 用 _binary_sdf（SDF 平滑遮罩） - sdf = os.path.join(volume_dir, f'{level}_binary_sdf.nii.gz') + # 未旋轉檔：新世代在 crop/ 子資料夾、舊世代在頂層 + sdf = level_file_path(volume_dir, level, 'binary_sdf') binary_path = sdf if os.path.exists(sdf) \ - else os.path.join(volume_dir, f'{level}_binary.nii.gz') + else level_file_path(volume_dir, level, 'binary') # azi = azimuth_rotation(binary_path) # res = analyze_vertebral_tilt_contour(binary_path, edge_type='superior', show_plot=False, debug=False) @@ -449,6 +453,14 @@ def main(): for _, g, vid, level, side, _, err in fail: logger.error(f'[GPU {g}] {vid} {level} {side}: {err}') + # 收尾：螺絲位置映回原 CT 空間 -> Output_dir///cbt.nii.gz + # （label 1-10 = L1L L1R L2L L2R ... L5L L5R；無 side 結果的 volume 跳過） + for vid in volumes: + try: + xfr_cbt_native.write_volume_cbt(vid, run_id) + except Exception as e: + logger.error(f'[CBT-NATIVE] {vid}: {e}') + if __name__ == '__main__': main() \ No newline at end of file diff --git a/xfr_inverse_transform.py b/xfr_inverse_transform.py new file mode 100644 index 0000000..9ae30ba --- /dev/null +++ b/xfr_inverse_transform.py @@ -0,0 +1,507 @@ +#!/home/xfr/.conda/envs/cbt/bin/python +"""把標準化 grid 上的 segmentation mask 反向映射回 data_root 原始 CT grid。 + +座標參數來自 xfr_preprocess 記錄的 /transform.json（每 level： +裁切 box / standardize 翻軸 / ap_flip / 原 CT 幾何 / rotated 的 R、center、 +start）。純 index 空間映射（不依賴 standardized 輸出的物理 header），NN +（order=0），輸出幾何 = 原始 CT，pixel type = 輸入 mask。 + +支援的 source grid（--source，預設由檔名自動判定）： + rotated /rotated/_*.nii.gz 上的任何 mask（含 label map、 + 自己跑的分割結果；整格同幾何） + smd_resampled [/crop]/_smd_resampled.nii.gz（0.5mm 裁切） + binary_sdf [/crop]/_binary_sdf.nii.gz + roi [/crop]/_roi.nii.gz + binary_nn [/crop]/_binary_nn.nii.gz + binary [/crop]/_binary.nii.gz（原解析度） + smd [/crop]/_smd.nii.gz（原解析度 SMD，裁切 +4 margin） + （[/crop]：新世代未旋轉檔在 crop/ 子資料夾、舊世代在頂層，自動判定） +mask 尺寸須與該 source 檔一致（同 grid）；level 取檔名前綴（L1..L5/T#..）。 + +用法： + xfr_inverse_transform.py VOL_DIR MASK [--source auto] [--level auto] + [--out PATH] [--orig CT] [--rebuild] [--chunk 16] + + VOL_DIR 標準化 volume 目錄（含 transform.json 或可 --rebuild） + MASK mask nifti（label / mask / 浮點分數皆可，dtype 保留） + --rebuild transform.json 缺失時從 disk + data_root + metadata db 重建 + （box 由 disk 檔 origin 回復（精確；缺檔時 fallback native + label 重算）、rotated R/center 重算並做正向 identity 驗證、 + fstart 由 disk 幾何回復、ap_flip 以骨頭/native label 對位 IoU + 檢測（db 值僅先驗——舊世代與新世代 disk 資料的 ap_flip 可能 + 不同）；結果寫回 transform.json）。舊世代（xfr-2 等）用這個。 + +輸出預設：_to_original.nii.gz（與 MASK 同目錄）。 +""" + +import argparse +import logging +import os +import re +import sys +import time + +import numpy as np +import SimpleITK as sitk + +from config.constant import LABEL_MAP +from imaging.orientation import best_symmetry_plane, best_upper_endplate_plane +from imaging.resample import resample_img +from imaging.segmentation import _largest_cc_bbox +from imaging.transforms import (SOURCES, build_volume_meta, img_geom, + level_file_path, load_transform, + mask_to_original, margined_box, + resampled05_geom_from_original, save_transform, + std_flip_axes_for_direction, transform_path) +from visualization.res_bone_figure import compute_normalizing_rotation + +from xfr_orig_labels import (_fstart_from_geometries, _forward_identity_check, + _load_unrotated_template, _bone_mask, + find_original_ct, load_ap_flip) + +DATA_ROOT = '/mnt/1220/Public/dataset/Spine/CTSpine1K/data/' +LABEL_ROOT = '/mnt/1220/Public/dataset/Spine/CTSpine1K/label/' +LABEL_SUBDIRS = ('conlon', 'COVID-19', 'HNSCC-3DCT-RT_neck', 'Liver') + +LEVEL_NAMES = set(LABEL_MAP.values()) +_LABEL_ID = {v: int(k) for k, v in LABEL_MAP.items()} + +LOG_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'logs') + +logger = logging.getLogger('xfr_inverse_transform') + + +def setup_tee(log_path): + """console 與 log 檔同時輸出（append、line-buffered）。""" + class _Tee: + def __init__(self, console, fh): + self.console, self.fh, self.buf = console, fh, '' + + def write(self, data): + if not data: + return + self.buf += data + while True: + idx_n, idx_r = self.buf.find('\n'), self.buf.find('\r') + idx = min([i for i in (idx_n, idx_r) if i != -1], default=-1) + if idx == -1: + break + line, self.buf = self.buf[:idx], self.buf[idx + 1:] + self.console.write(line + '\n') + self.fh.write(line + '\n') + + def flush(self): + if self.buf: + line, self.buf = self.buf, '' + self.console.write(line + '\n') + self.fh.write(line + '\n') + self.console.flush() + self.fh.flush() + + fh = open(log_path, 'a', buffering=1) + sys.stdout = _Tee(sys.stdout, fh) + sys.stderr = _Tee(sys.stderr, fh) + + +def _level_of(stem): + """檔名 stem（無 .nii.gz）的前綴若為 LABEL_MAP 的 level 名則回傳，否則 None。""" + prefix = stem.split('_')[0] + return prefix if prefix in LEVEL_NAMES else None + + +def detect_source_level(mask_path, vol_dir): + """由 MASK 相對於 vol_dir 的路徑推 (source, level)；推不出時 raise。""" + rel = os.path.relpath(os.path.abspath(mask_path), os.path.abspath(vol_dir)) + parts = rel.split(os.sep) + if len(parts) < 1 or not parts[-1].endswith('.nii.gz'): + raise ValueError(f'{mask_path}：需為 .nii.gz') + stem = parts[-1][:-len('.nii.gz')] + level = _level_of(stem) + if level is None: + raise ValueError(f'{mask_path}：檔名前綴推不出 level（預期 /C#>_...）') + if len(parts) >= 2 and parts[-2] == 'rotated': + return 'rotated', level + suffix = stem.split('_', 1)[1] if '_' in stem else '' + if suffix in SOURCES and suffix != 'rotated': + return suffix, level + raise ValueError(f'{mask_path}：推不出 source（目錄需為 .../rotated/ 或檔名需為 ' + f'_{{{"|".join(SOURCES)}}}.nii.gz）') + + +def discover_levels(vol_dir): + """volume 目錄（含 crop/、rotated/）出現的所有 level，依 label id 排序。 + 未旋轉檔：新世代在 crop/ 子資料夾、舊世代在頂層，兩者都掃。""" + lv = set() + for d in (vol_dir, os.path.join(vol_dir, 'crop'), os.path.join(vol_dir, 'rotated')): + if not os.path.isdir(d): + continue + for f in os.listdir(d): + if f.endswith('.nii.gz'): + l = _level_of(f[:-len('.nii.gz')]) + if l is not None: + lv.add(l) + return sorted(lv, key=lambda l: _LABEL_ID[l]) + + +def find_native_label(name): + """LABEL_ROOT 各子目錄找 _seg.nii.gz；回傳 path 或 None。""" + for sub in LABEL_SUBDIRS: + p = os.path.join(LABEL_ROOT, sub, f'{name}_seg.nii.gz') + if os.path.isfile(p): + return p + return None + + +def _cc_bbox05(mask_img): + """最大 26-CC bbox：mask_img（0/1）-> [x0,y0,z0,xs,ys,zs]。""" + cc = _largest_cc_bbox(mask_img) + if cc is None: + raise RuntimeError('empty level mask after largest-CC extraction') + return cc[1] + + +def _disk_box(orig_img, disk_img, flips, spacing): + """由 disk 檔 origin 回復該檔的裁切 box [x0,y0,z0,xs,ys,zs]（x,y,z 序）。 + + standardize_affine 翻軸的檔：disk origin (LPS) = -pre-standardize + origin（資料被鏡射到世界原點鏡射面）；未翻軸：disk origin = + pre-standardize origin = O_full + box0·(d·h)。spacing 為該檔每軸 mm + （0.5mm 檔 = [0.5,0.5,0.5]；原解析度檔 = 原 CT spacing）。 + 非整數（>1e-3）時 raise（幾何模型不符）。""" + o_full = np.array(orig_img.GetOrigin(), dtype=float) + o_disk = np.array(disk_img.GetOrigin(), dtype=float) + d = np.diag(np.array(orig_img.GetDirection(), dtype=float).reshape(3, 3)) + s = np.asarray(spacing, dtype=float) + if s.ndim == 0: + s = np.array([float(s)] * 3) + n = np.array(disk_img.GetSize(), dtype=float) + x0 = [((-o_disk[i] - o_full[i]) if i in flips else (o_disk[i] - o_full[i])) + / (d[i] * s[i]) for i in range(3)] + for i, v in enumerate(x0): + if abs(v - round(v)) > 1e-3: + raise ValueError(f'{disk_img} 翻軸軸 {i} 隱含 box 非整數 ({v:.4f})，' + f'幾何模型不符') + return [int(round(v)) for v in x0] + [int(v) for v in n] + + +def _recompute_box05(lb_arr, lab_id, lb_img, ct_img): + """fallback：native label 重算（未翻轉慣例）的最大 26-CC 0.5mm bbox + （與 seg_bone 的 bbox2 同定義）。舊世代 ap_flip=True 個案的 disk 裁切 + 位置與此可能有 ≤1 voxel 差（舊版 flip 路徑），故優先用 _disk_box。""" + m0_img = sitk.GetImageFromArray((lb_arr == lab_id).astype(np.uint8)) + m0_img.CopyInformation(lb_img) + bin_lin = resample_img(sitk.Cast(m0_img, sitk.sitkFloat32)) + m_full = (sitk.GetArrayFromImage(bin_lin) > 0.5).astype(np.uint8) + mf_img = sitk.GetImageFromArray(m_full) + mf_img.CopyInformation(bin_lin) + return [int(v) for v in _cc_bbox05(mf_img)] + + +def rebuild_level_boxes(vol_dir, level, lab_id, lb_img, lb_arr, orig_img, + flips): + """回復 level 各檔的裁切 box：優先由 disk 檔 origin 回復（精確、含舊世代 + ap_flip 差異）；該檔缺失時 fallback native label 重算（未翻轉慣例）並 + 警告。0.5mm 檔（smd_resampled/binary_sdf/roi 共用 box2、binary_nn 自己 + 的 nn box）+ 原解析度檔（binary、smd）都在此。""" + name = os.path.basename(os.path.abspath(vol_dir)) + boxes = {} + shared = None + for sname in ('smd_resampled', 'binary_sdf', 'roi'): + p = level_file_path(vol_dir, level, sname) + if os.path.exists(p): + try: + b = _disk_box(orig_img, sitk.ReadImage(p), flips, 0.5) + if shared is None: + shared = b + elif b != shared and max(abs(a - c) for a, c in zip(b, shared)) > 0: + logger.warning(f'{name} {level}: {sname} disk box ' + f'{b} != smd_resampled {shared}（應相同）') + except ValueError as e: + logger.warning(f'{name} {level}: {sname} disk box 回復失敗' + f'（{e}）') + if shared is None: + shared = _recompute_box05(lb_arr, lab_id, lb_img, orig_img) + shared = list(shared) + logger.warning(f'{name} {level}: 無 0.5mm 檔案可回復 disk box，' + f'fallback native label 重算 {shared}（ap_flip=True 舊' + f'世代時可能 ≤1 voxel 偏差）') + boxes['smd_resampled'] = list(shared) + boxes['binary_sdf'] = list(shared) + boxes['roi'] = list(shared) + + nn_path = level_file_path(vol_dir, level, 'binary_nn') + if os.path.exists(nn_path): + boxes['binary_nn'] = _disk_box(orig_img, sitk.ReadImage(nn_path), + flips, 0.5) + else: + nn_img = resample_img(lb_img, is_label=True) + mn = (sitk.GetArrayFromImage(nn_img) == lab_id).astype(np.uint8) + mn_img = sitk.GetImageFromArray(mn) + mn_img.CopyInformation(nn_img) + nn_res = _largest_cc_bbox(mn_img) + if nn_res is not None: + boxes['binary_nn'] = [int(v) for v in nn_res[1]] + else: + logger.warning(f'{name} {level}: NN fallback 無連通區域') + + for sname, fb in (('binary', 'm0'), ('smd', 'smd4')): + p = level_file_path(vol_dir, level, sname) + try: + boxes[sname] = _disk_box(orig_img, sitk.ReadImage(p), flips, + orig_img.GetSpacing()) + except (ValueError, OSError) as e: + m0_img = sitk.GetImageFromArray((lb_arr == lab_id).astype(np.uint8)) + m0_img.CopyInformation(lb_img) + bbox_orig = [int(v) for v in _cc_bbox05(m0_img)] + boxes[sname] = bbox_orig if fb == 'm0' else \ + margined_box(bbox_orig, orig_img.GetSize(), margin=4) + logger.warning(f'{name} {level}: {sname} disk box 回復失敗' + f'（{e}），fallback native label') + return boxes + + +def rebuild_rotated_section(vol_dir, level): + """重算 rotated 座標：template 選擇與 R/center 重算（與 _write_rotated_level + 同路徑，確定性）、fstart 由 rotated/target 檔的 disk 幾何回復，並做正向 + identity 驗證（mismatch > 1e-4 raise）。回傳 section dict 或 + (None, reason)。""" + tpl = _load_unrotated_template(vol_dir, level) + if tpl is None: + return None, 'no unrotated template' + t_img, t_arr, is_smd = tpl + rot_dir = os.path.join(vol_dir, 'rotated') + rot_path = None + for s in ('smd_resampled', 'binary_sdf', 'label'): + p = os.path.join(rot_dir, f'{level}_{s}.nii.gz') + if os.path.exists(p): + rot_path = p + break + if rot_path is None: + return None, 'no rotated file' + rot_img = sitk.ReadImage(rot_path) + fstart = _fstart_from_geometries(rot_img, t_img) + if fstart is None: + return None, 'rotated/unrotated geometry mismatch' + bin_arr = _bone_mask(t_arr, is_smd).astype(np.uint8) + sym = best_symmetry_plane(bin_arr) + symp = best_upper_endplate_plane(bin_arr) + R, c = compute_normalizing_rotation(bin_arr, sym, symp) + fwd = _forward_identity_check(vol_dir, level, R, c, t_img, t_arr, is_smd) + if fwd is not None and fwd > 1e-4: + raise RuntimeError(f'{level}: 正向重旋轉 mismatch {fwd:.3e} > 1e-4' + f'（重算 R/center 與生成參數不符；拒絕寫入）') + if fwd is None: + logger.warning(f'{level}: 無可比對的 rotated 檔，R/center 未經驗證') + return { + 'template': 'smd_resampled' if is_smd else 'binary_nn', + 'R': [[float(v) for v in row] for row in R], + 'center': [float(v) for v in c], + 'start': [int(v) for v in fstart], + 'size': [int(v) for v in rot_img.GetSize()], + }, ('forward_mismatch=%.1e' % fwd) if fwd is not None else 'unverified' + + +def _detect_ap_flip(vol_dir, meta, lb_arr, ap_db, max_levels=3): + """disk 資料是否含 ap_flip 的對位檢測：以「template 骨頭 -> 原始 grid」 + 與 native label 的 bone IoU 比對兩種假設，取高者。 + + metadata db 的 ap_flip 只反映【新 pipeline 世代】（09-07 起）的行為； + 舊世代（如 standardized-xfr-2，09-06 產出）disk 上沒有 y 翻轉，db 值 + （可能為 True）不可直接用於舊世代資料 —— 故重建時一律以對位結果為準， + db 值僅作為先驗（試驗順序 / 高 IoU 時短路）。""" + import copy + lvls = [l for l in meta['levels'] if 'rotated' in meta['levels'][l]] + name = os.path.basename(os.path.abspath(vol_dir)) + if not lvls: + chosen = bool(ap_db) if ap_db is not None else False + logger.warning(f'{name}: 無 rotated level 可驗證 ap_flip；' + f'{"用 db 值" if ap_db is not None else "缺 db 值"} {chosen}') + return chosen + order = [bool(ap_db), not bool(ap_db)] if ap_db is not None else [False, True] + best_hyp, best_iou = None, -1.0 + for hyp in order: + m = copy.deepcopy(meta) + m['ap_flip'] = hyp + ious = [] + for level in lvls[:max_levels]: + tpl = m['levels'][level]['rotated']['template'] + a = sitk.GetArrayFromImage( + sitk.ReadImage(level_file_path(vol_dir, level, tpl))) + bone = ((a < 0.5) if tpl == 'smd_resampled' else (a > 0)).astype(np.uint8) + if not bone.any(): + continue + out = mask_to_original(bone, m, level, tpl) + lab = (lb_arr == _LABEL_ID[level]) + on = out > 0 + union = int((on | lab).sum()) + ious.append(int((on & lab).sum()) / union if union else 0.0) + iou = float(np.mean(ious)) if ious else 0.0 + logger.info(f'{name}: ap_flip={hyp} -> 骨頭 native IoU={iou:.4f}') + if iou > best_iou: + best_iou, best_hyp = iou, hyp + if ap_db is not None and hyp == bool(ap_db) and iou >= 0.95: + break # 先驗假設已充分對位，免試另一半 + if ap_db is not None and best_hyp != bool(ap_db): + logger.warning(f'{name}: ap_flip 對位檢測={best_hyp} 與 metadata db=' + f'{bool(ap_db)} 不同（db 值屬新 pipeline 世代；disk 資料' + f'應為另一世代），採對位檢測結果') + if best_iou < 0.9: + logger.warning(f'{name}: ap_flip 對位檢測最佳 IoU={best_iou:.4f} < 0.9，' + f'座標鏈可能有誤，請人工核對') + return best_hyp + + +def rebuild_transform(vol_dir, orig_path=None): + """舊 volume（無 transform.json）：從 disk 輸出 + data_root 原始 CT + + native label + metadata db 重建完整 transform.json 並寫回。 + 回傳 meta。所有推導與 pipeline 記錄同式（box 由 disk 檔 origin 回復、 + 缺檔 fallback native label 重算、flip 軸由原 CT direction 推斷、 + rotated R/center 重算 + 正向 identity 驗證、ap_flip 以對位 IoU 檢測 + （db 值僅先驗，見 _detect_ap_flip））。""" + name = os.path.basename(os.path.abspath(vol_dir)) + if orig_path is None: + orig_path, _ = find_original_ct(name) + if orig_path is None: + raise FileNotFoundError(f'{DATA_ROOT}*/{name}.nii.gz 找不到') + orig_img = sitk.ReadImage(orig_path) + flip_axes = std_flip_axes_for_direction(orig_img.GetDirection()) + lb_path = find_native_label(name) + if lb_path is None: + raise FileNotFoundError(f'{LABEL_ROOT}*/{name}_seg.nii.gz 找不到') + lb_img = sitk.ReadImage(lb_path) + lb_arr = sitk.GetArrayFromImage(lb_img) + + levels = discover_levels(vol_dir) + if not levels: + raise ValueError(f'{vol_dir}: 找不到任何 _*.nii.gz') + level_entries = {} + for level in levels: + boxes = rebuild_level_boxes(vol_dir, level, _LABEL_ID[level], lb_img, + lb_arr, orig_img, flip_axes) + lv = {'label': _LABEL_ID[level], 'std_flip_axes': flip_axes, + 'boxes': boxes} + sec, note = rebuild_rotated_section(vol_dir, level) + if sec is not None: + lv['rotated'] = sec + logger.info(f'{name} {level}: rebuilt (rotated verified: {note})') + else: + logger.info(f'{name} {level}: rebuilt (no rotated section: {note})') + level_entries[level] = lv + orig_geom = img_geom(orig_img) + meta = build_volume_meta(name, orig_path, orig_geom, + resampled05_geom_from_original(orig_geom), + False, level_entries) + ap_db = load_ap_flip().get(name) + meta['ap_flip'] = _detect_ap_flip(vol_dir, meta, lb_arr, ap_db) + save_transform(vol_dir, meta) + logger.info(f'{name}: rebuild -> {transform_path(vol_dir)} ' + f'({len(level_entries)} level(s), ap_flip={meta["ap_flip"]}, ' + f'std_flip_axes={flip_axes})') + return meta + + +def _check_orig_geometry(meta, orig_img, src='recorded'): + rec = meta['original'] + if (list(orig_img.GetSize()) != rec['size'] + or any(abs(float(a) - b) > 1e-6 + for a, b in zip(orig_img.GetSpacing(), rec['spacing']))): + sys.exit(f'原始 CT 幾何與記錄不符（{src}）：size {orig_img.GetSize()} ' + f'spacing {orig_img.GetSpacing()} vs 記錄 {rec["size"]} ' + f'{rec["spacing"]}；座標鏈以記錄值為準，無法映射') + + +def main(): + parser = argparse.ArgumentParser( + description='Inverse-map a standardized-grid segmentation mask to the ' + 'original CT coordinates (data_root) via transform.json.') + parser.add_argument('vol_dir', help='standardized volume dir (含 transform.json)') + parser.add_argument('mask', help='mask nifti（在 source grid 上）') + parser.add_argument('--source', default=None, + choices=['auto'] + list(SOURCES), + help='source grid（default: 由檔名自動判定）') + parser.add_argument('--level', default=None, + help='level（default: 由檔名前綴自動判定）') + parser.add_argument('--out', default=None, + help='輸出 nifti（default: _to_original.nii.gz）') + parser.add_argument('--orig', default=None, + help='原始 CT path（default: transform.json 記錄值）') + parser.add_argument('--rebuild', action='store_true', + help='transform.json 缺失時重建（舊世代 volume）') + parser.add_argument('--chunk', type=int, default=16, + help='反向映射的 z 區塊 slice 數（記憶體控制, default 16）') + args = parser.parse_args() + + if not os.path.isdir(args.vol_dir): + sys.exit(f'vol_dir 不存在：{args.vol_dir}') + if not os.path.isfile(args.mask): + sys.exit(f'mask 不存在：{args.mask}') + + os.makedirs(LOG_DIR, exist_ok=True) + log_path = os.path.join(LOG_DIR, + f'xfr_inverse_transform_{time.strftime("%Y%m%d_%H%M%S")}.log') + setup_tee(log_path) + logging.basicConfig(level=logging.INFO, + format='%(asctime)s [%(levelname)s] %(name)s: %(message)s', + datefmt='%Y-%m-%d %H:%M:%S') + logger.info(f'Log file: {log_path}') + logger.info(f'Command: {sys.executable} {" ".join(sys.argv)}') + + tp = transform_path(args.vol_dir) + if os.path.exists(tp): + meta = load_transform(args.vol_dir) + elif args.rebuild: + meta = rebuild_transform(args.vol_dir, args.orig) + else: + sys.exit(f'{tp} 不存在；重跑 pipeline 或在既有 volume 上加 --rebuild') + + source = args.source if (args.source and args.source != 'auto') else None + level = args.level + if source is None or level is None: + a_source, a_level = detect_source_level(args.mask, args.vol_dir) + source = source if source is not None else a_source + level = level if level is not None else a_level + + name = meta.get('name') or os.path.basename(os.path.abspath(args.vol_dir)) + if level not in meta['levels']: + sys.exit(f'{level} 不在 transform.json（{sorted(meta["levels"])}）') + if source == 'rotated' and 'rotated' not in meta['levels'][level]: + sys.exit(f'{level} 無 rotated 記錄（未跑 _write_rotated_level 的 level ' + f'沒有 rotated grid 座標鏈）') + if source in ('smd_resampled', 'binary_sdf', 'roi', 'binary_nn', + 'binary', 'smd') and source not in meta['levels'][level]['boxes']: + sys.exit(f'{level} 無 {source} 記錄（該檔未產出）') + + if args.orig: + orig_img = sitk.ReadImage(args.orig) + _check_orig_geometry(meta, orig_img, '--orig') + else: + p = meta['original']['path'] + if not os.path.exists(p): + sys.exit(f'記錄的原始 CT 不存在：{p}（可用 --orig 指定）') + orig_img = sitk.ReadImage(p) + + mask_img = sitk.ReadImage(args.mask) + t0 = time.time() + arr = mask_to_original(mask_img, meta, level, source, chunk=args.chunk) + dt = time.time() - t0 + + n_on = int((arr != 0).sum()) + ct = sitk.GetArrayFromImage(orig_img) + if n_on: + hu = ct[arr != 0] + frac_body = float((hu > -100).mean()) + med_hu = float(np.median(hu)) + logger.info(f'alignment (info): nonzero={n_on}, in-body(HU>-100)=' + f'{frac_body:.3f}, median HU={med_hu:.0f}') + else: + logger.warning('反向映射結果全零（mask 可能全空或 grid 錯位）') + + out_path = args.out or (args.mask[:-len('.nii.gz')] + '_to_original.nii.gz') + from imaging.transforms import save_as_original + save_as_original(arr, orig_img, out_path) + logger.info(f'saved {out_path} ({name} {level} {source} -> original, ' + f'nonzero={n_on}, {dt:.1f}s)') + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/xfr_migrate_crop.py b/xfr_migrate_crop.py new file mode 100644 index 0000000..5b7ce12 --- /dev/null +++ b/xfr_migrate_crop.py @@ -0,0 +1,113 @@ +#!/home/xfr/.conda/envs/cbt/bin/python +"""一次性的 layout 遷移：把已產出 volume 的「未旋轉 檔」移入 crop/ 子資料夾。 + +新 layout（2026-09-08 起 xfr_preprocess 直接產出，無須再跑本腳本）： + /crop/_binary|_smd|_smd_resampled|_binary_sdf|_binary_nn|_roi.nii.gz + /crop/_planes.png （原 /lumbar/） + /rotated/_*.nii.gz （不變） + /transform.json、_labels.txt（頂層，不變） + +本腳本供改動前已產出、尚未重新處理的 volume 用（例：standardized-xfr-3）： + 1) / 頂層的 _*.nii.gz（L = LABEL_MAP level 名）移到 /crop/ + 2) /lumbar/*.png 移到 /crop/，lumbar/ 清空後刪除 +Idempotent：已遷就的 volume 重跑無作用；crop/ 已有同名檔時不覆蓋、跳過該檔 +（頂層檔保留，人工核對）。rotated/ 與其他頂層檔不動。 + +讀端（xfr_orig_labels / xfr_inverse_transform / xfr_cbt_native 等）以 +imaging.transforms.level_file_path 同時認 crop/ 與頂層，舊世代 +（standardized-xfr-2 等）不遷移也能讀；如需一致化可對該 root 跑本腳本。 + +Usage: + python xfr_migrate_crop.py [ROOT] [--dry-run] +""" +import argparse +import os +import shutil +import sys + +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) +from config.constant import LABEL_MAP + +DEFAULT_ROOT = '/mnt/1248/open2/cyrou/CBT/Seg/Resample/standardized-xfr-3' +LEVEL_SET = set(LABEL_MAP.values()) + + +def migrate_volume(vol_dir, dry_run=False): + """回傳 (moved, kept) 兩個清單（人讀的說明文字）。""" + crop_dir = os.path.join(vol_dir, 'crop') + moved, kept = [], [] + + for f in sorted(os.listdir(vol_dir)): + if not f.endswith('.nii.gz'): + continue + if f.split('_', 1)[0] not in LEVEL_SET: + kept.append(f) + continue + dst = os.path.join(crop_dir, f) + if os.path.exists(dst): + kept.append(f'{f}（crop/ 已有同名檔，未動）') + continue + src = os.path.join(vol_dir, f) + if dry_run: + moved.append(f'{f} (dry-run)') + else: + os.makedirs(crop_dir, exist_ok=True) + shutil.move(src, dst) + moved.append(f) + + lumbar_dir = os.path.join(vol_dir, 'lumbar') + if os.path.isdir(lumbar_dir): + for f in sorted(os.listdir(lumbar_dir)): + dst = os.path.join(crop_dir, f) + if os.path.exists(dst): + kept.append(f'lumbar/{f}（crop/ 已有同名檔，未動）') + continue + if dry_run: + moved.append(f'lumbar/{f} (dry-run)') + else: + os.makedirs(crop_dir, exist_ok=True) + shutil.move(os.path.join(lumbar_dir, f), dst) + moved.append(f'lumbar/{f}') + if not dry_run: + try: + os.rmdir(lumbar_dir) + except OSError: + kept.append('lumbar/（仍有其他檔，未刪）') + return moved, kept + + +def main(): + parser = argparse.ArgumentParser( + description='Migrate per-level pre-rotated files from / top level ' + 'to /crop/ (plus lumbar/*.png -> crop/).') + parser.add_argument('root', nargs='?', default=DEFAULT_ROOT, + help=f'standardized generation dir (default: {DEFAULT_ROOT})') + parser.add_argument('--dry-run', action='store_true', + help='Only report what would be moved.') + args = parser.parse_args() + + if not os.path.isdir(args.root): + sys.exit(f'not a directory: {args.root}') + vols = sorted(d for d in os.listdir(args.root) + if os.path.isdir(os.path.join(args.root, d))) + tag = '[dry-run] ' if args.dry_run else '' + n_done = n_empty = 0 + for i, vol in enumerate(vols, 1): + vol_dir = os.path.join(args.root, vol) + moved, kept = migrate_volume(vol_dir, dry_run=args.dry_run) + if not moved: + n_empty += 1 + print(f'{tag}[{i}/{len(vols)}] - {vol} (no file to move)') + continue + n_done += 1 + print(f'{tag}[{i}/{len(vols)}] ok {vol} moved {len(moved)} file(s)') + for m in moved: + print(f' -> {m}') + for k in kept: + print(f' !! kept: {k}') + print(f'\n{tag}done: {n_done} volume(s) migrated, {n_empty} untouched ' + f'out of {len(vols)}') + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/xfr_orig_labels.py b/xfr_orig_labels.py new file mode 100644 index 0000000..5e77361 --- /dev/null +++ b/xfr_orig_labels.py @@ -0,0 +1,628 @@ +#!/home/xfr/.conda/envs/cbt/bin/python + +"""把每個 lumbar level 的 rotated/{L}_label.nii.gz（旋轉座標， +1=VBODY 椎體、2=SP 棘突、3=other bone）反旋轉回原始座標， +把所有 level 合併成「每 patient 一個」.nii.gz，存到 CTSpine1K +原始 dataset 下面（座標系與 data/ 裡的原始影像檔一致）。 + +做法： + 1) 依 level 重算對齊旋轉 (R, c_xyz)（與 _write_rotated_level 完全 + 同路徑：_smd_resampled 有效時以 <0.5 取骨頭 mask，否則 _binary_nn + 以 >0），再從旋轉 label 檔與未旋轉模板檔的 origin 差回復裁切 + 起點 fstart（(O_rot - O_tpl) 換元、必須為整數）。 + 2) NN 反向映射把 label 映回未旋轉 0.5mm 模板 grid。旋轉是 index 空間 + 剛性變換（物理座標系統只平移 fstart），所以換元後只是單純的 + R 旋轉，不需要重新估算任何平面。 + 3) std（0.5mm 標準化空間）-> 原始 CT 物理座標的對角仿射 mapping： + standardize_affine 每翻一個軸就把整份資料平移到世界原點的鏡射面 + （b_i = -(2O_i + d_i·0.5·(L_i-1))，翻軸條件 x: d_x=+1、y: d_y=+1、 + z: d_z=-1），ap_flip（metadata db）再令 y 軸鏡射。全部參數都可從 + 原始 CT 幾何 + ap_flip 算出，見 _std_to_orig_affine。 +4) NN 重取樣到原始 CT grid（/mnt/1220/Public/dataset/Spine/CTSpine1K/ + data/<子目錄>/.nii.gz 的幾何），依 level 重編號後合併： + Lx 的 VBODY = 2x-1、SP = 2x（L1..L6 -> 1..12），0 = 背景； + rotated label 的 3（other bone）不納入。 + 每 level 驗證重算的 (R, c)：主要用正向 identity —— 以 _write_rotated_level + 完全同參數重旋轉未旋轉 check 檔（_smd_resampled 三線性、缺則 + _binary_nn NN），與存檔旋轉檔逐體素比對，mismatch > 1e-4 跳過 + （精確驗證 (R,c) 與該輸出生成版本一致，不受薄結構影響）；無可比檔 + 時 fallback 用 round-trip 包含率（旋轉 check 檔反向旋轉後與未旋轉 + bone mask 互相包含於 2-voxel 膨脹內，min 方向 < --min-contain + 跳過該 level）。 + 寫檔前對位檢查：合併 label 落在原始 CT 體內（HU > -100）比例 < 0.98 + 或 median HU < 100 時判定錯位、跳過。 + +輸出：/<子目錄>/.nii.gz（uint8），預設 dest = +/mnt/1220/Public/dataset/Spine/CTSpine1K/label_vbody_sp/， +檔名與原始影像檔相同。已存在的輸出預設跳過（--force 覆蓋）。 +""" + +import argparse +import json +import logging +import os +import sys +import time + +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) +from imaging.transforms import level_file_path +from visualization.res_bone_figure import (compute_normalizing_rotation, + rotated_grid, rotate_volume_to) + +DATA_ROOT = '/mnt/1220/Public/dataset/Spine/CTSpine1K/data/' +DATASET_ROOT = '/mnt/1220/Public/dataset/Spine/CTSpine1K/' +DEFAULT_DEST = os.path.join(DATASET_ROOT, 'label_vbody_sp') +# 預設依序掃的世代資料夾（先新後舊；同名 volume 以先處理者為準） +DEFAULT_ROOTS = [ + '/mnt/1248/open2/cyrou/CBT/Seg/Resample/standardized-xfr-3', + '/mnt/1248/open2/cyrou/CBT/Seg/Resample/standardized-xfr-2', +] +LUMBAR_LEVELS = ('L1', 'L2', 'L3', 'L4', 'L5', 'L6') + +LOG_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'logs') +METADATA_DB = os.path.join(os.path.dirname(os.path.abspath(__file__)), + 'xfr_image_metadata.json') + +logger = logging.getLogger('xfr_orig_labels') + + +def load_ap_flip(): + """metadata db（TinyDB 純 JSON：{'images': {doc_id: doc}}）的 + name -> ap_flip（bool）；db 缺失 / 解析失敗時回 {}。""" + try: + with open(METADATA_DB) as f: + data = json.load(f) + except (json.JSONDecodeError, OSError) as e: + logger.warning(f'metadata db {METADATA_DB} 讀取失敗：{e}') + return {} + tab = data.get('images') if isinstance(data, dict) else None + if not isinstance(tab, dict): + return {} + return {v['name']: bool(v.get('ap_flip', False)) + for v in tab.values() if isinstance(v, dict) and 'name' in v} + + +# ---------------------------------------------------------------- log 設定 +def setup_tee(log_path): + """console 與 log 檔同時輸出（append、line-buffered）。""" + class _Tee: + def __init__(self, console, fh): + self.console, self.fh, self.buf = console, fh, '' + + def write(self, data): + if not data: + return + self.buf += data + while True: + idx_n, idx_r = self.buf.find('\n'), self.buf.find('\r') + idx = min([i for i in (idx_n, idx_r) if i != -1], default=-1) + if idx == -1: + break + line, self.buf = self.buf[:idx], self.buf[idx + 1:] + self.console.write(line + '\n') + self.fh.write(line + '\n') + + def flush(self): + if self.buf: + line, self.buf = self.buf, '' + self.console.write(line + '\n') + self.fh.write(line + '\n') + self.console.flush() + self.fh.flush() + + fh = open(log_path, 'a', buffering=1) + sys.stdout = _Tee(sys.stdout, fh) + sys.stderr = _Tee(sys.stderr, fh) + + +# ---------------------------------------------------------------- 幾何工具 +def _axis_vectors(img): + """direction 為純對角（±1）時回傳 (signs (x,y,z), spacing)，否則 None。""" + d = np.asarray(img.GetDirection(), dtype=float).reshape(3, 3) + diag = np.diag(d) + if not (np.allclose(d, np.diag(diag), atol=1e-6) + and np.allclose(np.abs(diag), 1.0, atol=1e-6)): + return None + return diag, np.asarray(img.GetSpacing(), dtype=float) + + +def _fstart_from_geometries(rot_img, t_img): + """從旋轉檔與未旋轉模板檔的 origin 差回復 fstart（模板 index 系， + (x,y,z) 整數）；幾何不一致（非對角 direction、spacing / direction + 不同）或差不是整數時回傳 None。""" + ax_r, ax_t = _axis_vectors(rot_img), _axis_vectors(t_img) + if ax_r is None or ax_t is None: + return None + (d_r, s_r), (d_t, s_t) = ax_r, ax_t + if not (np.allclose(d_r, d_t, atol=1e-6) and np.allclose(s_r, s_t, atol=1e-6)): + return None + o_r = np.asarray(rot_img.GetOrigin(), dtype=float) + o_t = np.asarray(t_img.GetOrigin(), dtype=float) + fs = (o_r - o_t) / (d_t * s_t) + if not np.allclose(fs, np.round(fs), atol=1e-6): + return None + return np.round(fs).astype(int) + + +def _std_to_orig_affine(orig_img, t_img, ap_flip): + """std（0.5mm 裁切檔、standardize_affine 後、ap_flip 後）-> 原始 CT + 物理座標（LPS）：P_orig = A·P_std + b。回傳 (A (3,) 對角, b (3,))； + 原始 direction 非純對角時回傳 None。t_img = 該 level 的未旋轉標準化 + 裁切檔（disk 上的後 standardize 幾何）。 + + standardize_affine 是【逐檔（每 level 裁切）】處理：翻軸（nibabel： + affine 對角與平移同反號 + 資料 flip）在 worlds 中是整份資料平移 + t_i = 2·O_post_i - d_i·h·(L_i-1)（O_post/L 為 disk 上該檔的 + origin/尺寸、d 為原始 direction、h=0.5；翻軸條件 RAS 對角 < 0 + = x: d_x=+1、y: d_y=+1、z: d_z=-1），故 b_i = -t_i（翻軸的軸）。 + ap_flip 是 y 軸純資料鏡射（幾何不變，作用在整張 0.5mm full grid）： + y_ap = C - y_orig，C = 2·O_full_y + d_y·h·(L_full_y-1)（full grid + 幾何 = 原始 CT），發生在 standardize 之前。""" + d_full = np.asarray(orig_img.GetDirection(), dtype=float).reshape(3, 3) + d = np.diag(d_full) + if not (np.allclose(d_full, np.diag(d), atol=1e-6) + and np.allclose(np.abs(d), 1.0, atol=1e-6)): + return None + h = 0.5 + O_full = np.asarray(orig_img.GetOrigin(), dtype=float) + s = np.asarray(orig_img.GetSpacing(), dtype=float) + n = np.asarray(orig_img.GetSize(), dtype=float) # (x, y, z) + L_full = np.ceil(n * s / h - 1e-6) # 0.5mm full grid + O_post = np.asarray(t_img.GetOrigin(), dtype=float) + L = np.asarray(t_img.GetSize(), dtype=float) # 裁切檔尺寸 + flips = (d[0] == 1.0, d[1] == 1.0, d[2] == -1.0) + A = np.ones(3) + b = np.zeros(3) + for i in range(3): + if flips[i]: + b[i] = -(2.0 * O_post[i] - d[i] * h * (L[i] - 1.0)) + if ap_flip: + # y0 --AP 鏡射--> y1 = C - y0 --std 平移--> y_std = y1 + t_y + # => y0 = -y_std + (C + t_y) + A[1] = -1.0 + C_full = 2.0 * O_full[1] + d[1] * h * (L_full[1] - 1.0) + if flips[1]: + t_y = 2.0 * O_post[1] - d[1] * h * (L[1] - 1.0) + b[1] = C_full + t_y + else: + b[1] = C_full + return A, b + + +def _affine_sitk(A, b): + """對角 A + 平移 b 的 sitk 仿射（LPS 物理 -> LPS 物理）。 + p' = A·p + b；SetMatrix 為 row-major 3x3、SetTranslation 獨立。 + 自檢三軸 + 原點。""" + tr = sitk.AffineTransform(3) + tr.SetMatrix([float(A[0]), 0.0, 0.0, + 0.0, float(A[1]), 0.0, + 0.0, 0.0, float(A[2])]) + tr.SetTranslation([float(b[0]), float(b[1]), float(b[2])]) + for i in range(3): + p = [0.0, 0.0, 0.0] + p[i] = 123.0 + got = np.array(tr.TransformPoint(tuple(p))) + exp = np.array(b, dtype=float) + exp[i] = A[i] * 123.0 + b[i] + if not np.allclose(got, exp, atol=1e-6): + raise RuntimeError(f'affine self-check failed: axis {i} ' + f'{got} != {exp}') + return tr + + +def _rotate_back(arr, R, c_xyz, fstart, shape_zyx, order=0, cval=0.0): + """rotated grid 上的 field 反向旋轉到未旋轉模板 grid（shape_zyx）。 + rotated 檔 array 座標 a = 模板座標 V - fstart；反向映射 + V = R(W - c) + c，W 為模板 grid 上每格（x,y,z 向量，同 + compute_normalizing_rotation 慣例）。order / cval 同 map_coordinates。""" + nz, ny, nx = shape_zyx + idx = np.indices((nz, ny, nx), dtype=np.float64) + W = np.stack([idx[2], idx[1], idx[0]], axis=0) # (3, z, y, x) + c = np.asarray(c_xyz, dtype=np.float64)[:, None, None, None] + Rf = np.asarray(R, dtype=np.float64) + V = np.tensordot(Rf, W - c, axes=([1], [0])) + c + fs = np.asarray(fstart, dtype=np.float64)[:, None, None, None] + a = V - fs # rotated 檔 array 座標 + return map_coordinates(arr, [a[2], a[1], a[0]], + order=order, cval=cval, mode='constant') + + +# ---------------------------------------------------------------- level 處理 +def _load_unrotated_template(vol_dir, level): + """_write_rotated_level 的旋轉來源：_smd_resampled（有效時）否則 + _binary_nn。回傳 (img, arr, is_smd)；皆無則 None。 + 未旋轉檔位置經 level_file_path（新世代 crop/、舊世代頂層）。""" + smd_path = level_file_path(vol_dir, level, 'smd_resampled') + if os.path.exists(smd_path): + sa = sitk.GetArrayFromImage(sitk.ReadImage(smd_path)).astype(np.float32) + if sa.size > 0 and np.isfinite(sa).all() and sa.min() < 0.0 and sa.max() > 0.0: + return sitk.ReadImage(smd_path), sa, True + nn_path = level_file_path(vol_dir, level, 'binary_nn') + if os.path.exists(nn_path): + na = sitk.GetArrayFromImage(sitk.ReadImage(nn_path)) + if na.size > 0 and int((na > 0).sum()) > 0: + return sitk.ReadImage(nn_path), na, False + return None + + +def _bone_mask(t_arr, is_smd): + return (t_arr < 0.5) if is_smd else (t_arr > 0) + + +def _forward_identity_check(vol_dir, level, R, c_xyz, t_img, t_arr, is_smd): + """以 _write_rotated_level 完全同參數（fstart/fsize 由未旋轉 _binary_nn + 旋轉後 bbox±4、order / cval 同 writer）重旋轉、與存檔的旋轉檔逐體素 + 比對，回傳 mismatch 比例（(R,c) 相同時 ≈0，精確驗證、不受物件厚度 + 影響）；無可比檔（缺檔 / 形狀不符）時回 None。""" + rotated_dir = os.path.join(vol_dir, 'rotated') + if is_smd: + target = os.path.join(rotated_dir, f'{level}_smd_resampled.nii.gz') + src_arr, order, cval = t_arr, 1, float(t_arr.max()) + else: + target = os.path.join(rotated_dir, f'{level}_binary_nn.nii.gz') + src_arr, order, cval = (t_arr > 0).astype(np.uint8), 0, 0.0 + if not os.path.exists(target): + return None + start, size = rotated_grid(t_arr.shape, R, c_xyz, margin=4) + fstart = (int(start[0]), int(start[1]), int(start[2])) + fsize = (int(size[0]), int(size[1]), int(size[2])) + nn_path = level_file_path(vol_dir, level, 'binary_nn') + if os.path.exists(nn_path): + nn_img = sitk.ReadImage(nn_path) + nn_arr = (sitk.GetArrayFromImage( + sitk.Resample(nn_img, t_img, interpolator=sitk.sitkNearestNeighbor, + defaultPixelValue=0)) > 0).astype(np.uint8) + rot_nn_full = (rotate_volume_to(nn_arr, R, c_xyz, fstart, fsize, + order=0, cval=0.0) > 0.5).astype(np.uint8) + if int(rot_nn_full.sum()) > 0: + zz, yy, xx = np.where(rot_nn_full > 0) + m = 4 + z0 = max(0, int(zz.min()) - m) + z1 = min(rot_nn_full.shape[0] - 1, int(zz.max()) + m) + y0 = max(0, int(yy.min()) - m) + y1 = min(rot_nn_full.shape[1] - 1, int(yy.max()) + m) + x0 = max(0, int(xx.min()) - m) + x1 = min(rot_nn_full.shape[2] - 1, int(xx.max()) + m) + fstart = (fstart[0] + x0, fstart[1] + y0, fstart[2] + z0) + fsize = (x1 - x0 + 1, y1 - y0 + 1, z1 - z0 + 1) + out = rotate_volume_to(src_arr, R, c_xyz, fstart, fsize, order=order, + cval=cval) + saved = sitk.GetArrayFromImage(sitk.ReadImage(target)) + if out.shape != saved.shape: + return None + if order == 0: + return float(((out > 0.5) != (saved > 0.5)).mean()) + return float((np.abs(out.astype(np.float64) + - saved.astype(np.float64)) > 0.02).mean()) + + +def _containment_roundtrip(vol_dir, level, R, c_xyz, fstart, t_arr, is_smd, + it=2): + """fallback 驗證：旋轉 check 檔反向旋轉後與未旋轉 bone mask 互相包含於 + it-voxel 膨脹內（薄結構的 round-trip IoU 天生低，用包含率代替）。 + 回傳 min 方向包含率；無 check 檔時回 None。""" + from scipy.ndimage import binary_dilation + rotated_dir = os.path.join(vol_dir, 'rotated') + candidates = [ + (f'{level}_smd_resampled.nii.gz', 1), + (f'{level}_binary_sdf.nii.gz', 0), + (f'{level}_binary_nn.nii.gz', 0), + ] + if not is_smd: + candidates = candidates[1:] + B0 = _bone_mask(t_arr, is_smd) + for fname, order in candidates: + p = os.path.join(rotated_dir, fname) + if not os.path.exists(p): + continue + a = sitk.GetArrayFromImage(sitk.ReadImage(p)).astype(np.float64) + if a.size == 0: + continue + cval = float(a.max()) if fname.endswith('_smd_resampled.nii.gz') else 0.0 + b = _rotate_back(a, R, c_xyz, fstart, t_arr.shape, order=order, + cval=cval) + B1 = b < 0.5 if fname.endswith('_smd_resampled.nii.gz') else b > 0.5 + if not B1.any() or not B0.any(): + continue + da = binary_dilation(B1, iterations=it) + db = binary_dilation(B0, iterations=it) + return float(min((B1 & db).sum() / B1.sum(), (B0 & da).sum() / B0.sum())) + return None + + +def process_level(vol_dir, level, min_contain): + """單一 level：重算 (R, c)、回復 fstart、驗證、反向旋轉 label。 + 回傳 (unrotated_label_arr (z,y,x) uint8 0/1/2/3, t_img, info)；失敗回 None。""" + name = os.path.basename(vol_dir) + rot_lbl_path = os.path.join(vol_dir, 'rotated', f'{level}_label.nii.gz') + if not os.path.exists(rot_lbl_path): + return None + rot_img = sitk.ReadImage(rot_lbl_path) + rot_arr = sitk.GetArrayFromImage(rot_img).astype(np.uint8) + + tpl = _load_unrotated_template(vol_dir, level) + if tpl is None: + logger.warning(f'{name} {level}: 無可用未旋轉模板 ' + f'（_smd_resampled / _binary_nn），跳過') + return None + t_img, t_arr, is_smd = tpl + + fstart = _fstart_from_geometries(rot_img, t_img) + if fstart is None: + logger.warning(f'{name} {level}: 旋轉 / 未旋轉幾何不一致，跳過') + return None + + bin_arr = _bone_mask(t_arr, is_smd).astype(np.uint8) + sym = best_symmetry_plane(bin_arr) + symp = best_upper_endplate_plane(bin_arr) + R, c_xyz = compute_normalizing_rotation(bin_arr, sym, symp) + + # 驗證重算的 (R, c)：主要用 forward identity（與生成時同參數重旋轉、 + # 逐體素比對存檔；(R,c) 相同時 ≈0）；無可比檔時 fallback 用 + # round-trip 包含率（min_contain 門檻）。 + fwd = _forward_identity_check(vol_dir, level, R, c_xyz, t_img, t_arr, is_smd) + contain = None + if fwd is not None: + if fwd > 1e-4: + logger.warning(f'{name} {level}: 正向重旋轉 mismatch {fwd:.4f} ' + f'> 1e-4，(R,c) 與生成參數不符，跳過') + return None + else: + contain = _containment_roundtrip(vol_dir, level, R, c_xyz, fstart, + t_arr, is_smd) + if contain is None: + logger.warning(f'{name} {level}: 無 check 檔可驗證，跳過') + return None + if contain < min_contain: + logger.warning(f'{name} {level}: round-trip 包含率 {contain:.4f} ' + f'< {min_contain}，跳過該 level') + return None + + unrot = _rotate_back(rot_arr.astype(np.float64), R, c_xyz, fstart, + t_arr.shape).astype(np.uint8) + n_lbl = int((unrot > 0).sum()) + B0 = _bone_mask(t_arr, is_smd) + bone_frac = (int(((unrot > 0) & B0).sum()) / n_lbl) if n_lbl else 0.0 + info = { + 'forward_mismatch': fwd, + 'containment': contain, + 'fstart': [int(v) for v in fstart], + 'bone_frac': bone_frac, + 'rotated': {'vbody': int((rot_arr == 1).sum()), + 'sp': int((rot_arr == 2).sum()), + 'other': int((rot_arr == 3).sum())}, + 'unrotated': {'vbody': int((unrot == 1).sum()), + 'sp': int((unrot == 2).sum()), + 'other': int((unrot == 3).sum())}, + } + return unrot, t_img, info + + +def _resample_nearest(ref_img, src_img, transform=None): + r = sitk.ResampleImageFilter() + r.SetReferenceImage(ref_img) + r.SetInterpolator(sitk.sitkNearestNeighbor) + r.SetDefaultPixelValue(0) + if transform is not None: + r.SetTransform(transform) + return r.Execute(src_img) + + +def process_volume(vol_dir, orig_img, key, dest_dir, min_contain, force, + ap_flip): + """所有 level 反旋轉 + 合併 + 重取樣到原始 CT grid，存「一個」.nii.gz。 + 回傳 (out_path, status, {level: info}, vol_info)；失敗時 out_path=None、 + status 為原因。""" + name = os.path.basename(vol_dir) + nxi, nyi, nzi = orig_img.GetSize() # sitk GetSize = (x, y, z) + merged = np.zeros((nzi, nyi, nxi), dtype=np.uint8) # array = (z, y, x) + level_info = {} + maps = [] + + for li, level in enumerate(LUMBAR_LEVELS): + res = process_level(vol_dir, level, min_contain) + if res is None: + continue + unrot, t_img, info = res + # std->orig mapping 依該 level 裁切檔幾何逐 level 算 + aff = _std_to_orig_affine(orig_img, t_img, bool(ap_flip)) + if aff is None: + logger.warning(f'{name} {level}: 原始 CT direction 非純對角，' + f'跳過') + continue + A, b = aff + # Resample 的 transform 方向是 dest(原始) physical -> src(std) + # physical：P_std = A·(P_orig - b) = A·P_orig - A*b（A 對角 ±1） + tr = _affine_sitk(A, -A * b) + # 重編號：Lx VBODY=2x-1、SP=2x（x = L1..L6 的 1-based）； + # other(3) 不納入 + new = np.zeros_like(unrot) + new[unrot == 1] = 2 * li + 1 + new[unrot == 2] = 2 * li + 2 + maps.append({'A': [float(v) for v in A], + 'b': [round(float(v), 3) for v in b]}) + if not new.any(): + logger.info(f'{name} {level}: 無 VBODY / SP voxel，不計入') + continue + lvl_img = sitk.GetImageFromArray(new) + lvl_img.CopyInformation(t_img) + arr = sitk.GetArrayFromImage(_resample_nearest(orig_img, lvl_img, tr)) + merged = np.where(arr > 0, arr, merged) + info['orig'] = {'vbody': int((arr == 2 * li + 1).sum()), + 'sp': int((arr == 2 * li + 2).sum())} + level_info[level] = info + check = (f'fwd_mismatch={info["forward_mismatch"]:.1e}' + if info.get('forward_mismatch') is not None + else f'containment={info["containment"]:.4f}') + logger.info(f'{name} {level}: {check} ' + f'bone_frac={info["bone_frac"]:.4f} ' + f'std->orig A={np.round(A, 3)} b={np.round(b, 2)} ' + f'orig vbody={info["orig"]["vbody"]} ' + f'sp={info["orig"]["sp"]} (labels {2 * li + 1}/{2 * li + 2})') + + if not level_info: + return None, 'no usable level', {}, {} + if int(merged.sum()) == 0: + return None, 'merged label empty after resample', level_info, {} + + # 對位檢查：label 應落在原始 CT 體內（HU > -100）且主要是骨頭 + #（median HU >= 100；鬆質骨 HU 可 < 100，故不用 bone 比例作硬門檻） + ct = sitk.GetArrayFromImage(orig_img) + hu = ct[merged > 0] + frac_body = float((hu > -100).mean()) + med_hu = float(np.median(hu)) + frac_bone = float((hu > 100).mean()) + vol_info = {'in_body_frac': frac_body, 'bone_frac': frac_bone, + 'median_hu': med_hu, 'std_to_orig': maps, 'ap_flip': bool(ap_flip)} + if frac_body < 0.98 or med_hu < 100: + logger.warning(f'{name}: 對位檢查 in_body={frac_body:.4f} ' + f'median HU={med_hu:.0f}，疑似錯位，跳過') + return None, f'alignment check failed ' \ + f'(in_body={frac_body:.3f}, median_hu={med_hu:.0f})', \ + level_info, vol_info + logger.info(f'{name}: alignment OK in_body={frac_body:.4f} ' + f'bone_frac={frac_bone:.4f} median HU={med_hu:.0f}') + + out_path = os.path.join(dest_dir, key, f'{name}.nii.gz') + if os.path.exists(out_path) and not force: + return out_path, 'exists (skipped)', level_info, vol_info + os.makedirs(os.path.dirname(out_path), exist_ok=True) + out_img = sitk.GetImageFromArray(merged) + out_img.CopyInformation(orig_img) + sitk.WriteImage(out_img, out_path) + logger.info(f'{name}: saved {out_path} (levels={list(level_info)})') + return out_path, 'written', level_info, vol_info + + +def find_original_ct(name): + """在 DATA_ROOT 子目錄裡找 .nii.gz；回傳 (path, 子目錄) 或 + (None, None)。""" + if not os.path.isdir(DATA_ROOT): + return None, None + for key in sorted(os.listdir(DATA_ROOT)): + p = os.path.join(DATA_ROOT, key, f'{name}.nii.gz') + if os.path.isfile(p): + return p, key + return None, None + + +def list_volumes_with_labels(root): + """root 下有 rotated/*_label.nii.gz 的 volume 目錄（依名稱排序）。""" + vols = [] + if not os.path.isdir(root): + return vols + for d in sorted(os.listdir(root)): + rot = os.path.join(root, d, 'rotated') + if os.path.isdir(rot) and any( + f.endswith('_label.nii.gz') for f in os.listdir(rot)): + vols.append(os.path.join(root, d)) + return vols + + +def main(): + parser = argparse.ArgumentParser( + description='Merge unrotated lumbar VBODY/SP labels into one ' + 'original-CT coordinate .nii.gz per patient.') + parser.add_argument('--roots', nargs='*', default=DEFAULT_ROOTS, + help='Output generation dir(s) to scan (newest first).') + parser.add_argument('--dest', default=DEFAULT_DEST, + help=f'Destination dir (default: {DEFAULT_DEST})') + parser.add_argument('--force', action='store_true', + help='Overwrite existing outputs.') + parser.add_argument('--min-contain', type=float, default=0.85, + help='Fallback per-level check: skip a level whose ' + 'unrotation round-trip containment is below ' + 'this (default: 0.85).') + parser.add_argument('--name', default=None, + help='Process only this volume name (debug).') + args = parser.parse_args() + + os.makedirs(LOG_DIR, exist_ok=True) + log_path = os.path.join(LOG_DIR, f'xfr_orig_labels_{time.strftime("%Y%m%d_%H%M%S")}.log') + setup_tee(log_path) + logging.basicConfig(level=logging.INFO, + format='%(asctime)s [%(levelname)s] %(name)s: %(message)s', + datefmt='%Y-%m-%d %H:%M:%S') + logger.info(f'Log file: {log_path}') + logger.info(f'Command: {sys.executable} {" ".join(sys.argv)}') + logger.info(f'Roots: {args.roots}') + logger.info(f'Dest: {args.dest}') + logger.info('Label scheme: Lx VBODY=2x-1, SP=2x (L1:1/2, L2:3/4, ' + 'L3:5/6, L4:7/8, L5:9/10, L6:11/12); other bone excluded; ' + 'uint8; geometry = original CT') + + os.makedirs(args.dest, exist_ok=True) + ap_flip_db = load_ap_flip() + logger.info(f'ap_flip metadata: {len(ap_flip_db)} entries ' + f'({sum(1 for v in ap_flip_db.values() if v)} flipped)') + done, summary, no_ap = set(), {}, set() + summary_path = os.path.join(args.dest, 'summary.json') + if os.path.exists(summary_path): + try: + with open(summary_path) as f: + summary = json.load(f) + except (json.JSONDecodeError, OSError): + summary = {} + + for root in args.roots: + vols = list_volumes_with_labels(root) + if args.name is not None: + vols = [v for v in vols if os.path.basename(v) == args.name] + logger.info(f'Scanning {root}: {len(vols)} volume(s) with rotated labels') + for vol_dir in vols: + name = os.path.basename(vol_dir) + if name in done: + logger.info(f'[skip] {name}: 已由較新世代處理') + continue + done.add(name) + + orig_path, key = find_original_ct(name) + if orig_path is None: + logger.warning(f'[skip] {name}: 找不到原始 CT ' + f'({DATA_ROOT}*/{name}.nii.gz)') + continue + ap_flip = ap_flip_db.get(name) + if ap_flip is None and name not in no_ap: + no_ap.add(name) + logger.warning(f'{name}: metadata db 缺 ap_flip，假設 False ' + f'（若為 AP 翻轉個案將對位失敗並被抓出）') + ap_flip = False + out_path = os.path.join(args.dest, key, f'{name}.nii.gz') + if os.path.exists(out_path) and not args.force: + logger.info(f'[skip] {name}: {out_path} 已存在') + continue + + try: + orig_img = sitk.ReadImage(orig_path) + t0 = time.time() + out, status, level_info, vol_info = process_volume( + vol_dir, orig_img, key, args.dest, args.min_iou, args.force, + ap_flip) + dt = time.time() - t0 + common = {'status': status, 'source_root': root, + 'vol_info': vol_info, + 'updated': time.strftime('%Y-%m-%d %H:%M:%S')} + if out is None: + logger.warning(f'[skip] {name}: {status} ({dt:.1f}s)') + summary[name] = {**common, 'dest_root': DATA_ROOT, + 'levels': level_info} + else: + logger.info(f'[done] {name}: {status} in {dt:.1f}s') + summary[name] = {**common, 'path': out, 'key': key, + 'levels': level_info} + except Exception as e: + logger.warning(f'[error] {name}: {e}') + summary[name] = {'status': f'error: {e}', 'source_root': root, + 'updated': time.strftime('%Y-%m-%d %H:%M:%S')} + + with open(summary_path, 'w') as f: + json.dump(summary, f, indent=2, ensure_ascii=False) + logger.info('All done.') + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/xfr_plot_level.py b/xfr_plot_level.py index 588268f..8a45517 100644 --- a/xfr_plot_level.py +++ b/xfr_plot_level.py @@ -17,6 +17,7 @@ import argparse import os import sys +from imaging.transforms import level_file_path from visualization.res_bone_figure import render_bone_figure standardized_dir = '/mnt/1248/open2/cyrou/CBT/Seg/Resample/standardized-xfr/' @@ -72,10 +73,11 @@ def main(): for i, (vid, lvl) in enumerate(tasks, 1): vol_dir = os.path.join(args.dir, vid) # _binary.nii.gz 現為原解析度；0.5mm 用 _binary_sdf（SDF 平滑遮罩） - sdf_path = os.path.join(vol_dir, f'{lvl}_binary_sdf.nii.gz') + # 未旋轉檔：新世代在 crop/ 子資料夾、舊世代在頂層 + sdf_path = level_file_path(vol_dir, lvl, 'binary_sdf') binary_path = sdf_path if os.path.exists(sdf_path) \ - else os.path.join(vol_dir, f'{lvl}_binary.nii.gz') - cortical_path = os.path.join(vol_dir, f'{lvl}_cortical.nii.gz') + else level_file_path(vol_dir, lvl, 'binary') + cortical_path = level_file_path(vol_dir, lvl, 'cortical') path = render_bone_figure(vid, lvl, binary_path, cortical_path, base_folder=args.output) if path is None: diff --git a/xfr_preprocess.py b/xfr_preprocess.py index 6db50cd..20ef746 100644 --- a/xfr_preprocess.py +++ b/xfr_preprocess.py @@ -22,6 +22,7 @@ from visualization.res_bone_figure import (render_bone_figure, rotated_sitk_image_at, _rotate_plane_params, _shift_plane_params) +from imaging.transforms import merge_rotated_into_transform, level_file_path data_root = '/mnt/1220/Public/dataset/Spine/CTSpine1K/data/' label_root = '/mnt/1220/Public/dataset/Spine/CTSpine1K/label/' @@ -58,9 +59,20 @@ def _upsert_by_name(table, name, meta): table.update(entry, q) +def _looks_like_tinydb_table(v): + """TinyDB table 資料格式：{數字字串 doc_id: doc}。""" + return (isinstance(v, dict) and bool(v) + and all(isinstance(k, str) and k.isdigit() for k in v)) + + def _migrate_legacy_metadata(path): - """舊版純 JSON 檔（{name: {...}}、無 TinyDB 的 _default 結構）： - 改名成 .legacy- 備份，內容匯入新 TinyDB。""" + """舊版純 JSON 檔（{影像 name: {...}}、無 TinyDB table 結構）： + 改名成 .legacy- 備份，內容匯入新 TinyDB。 + + 注意：TinyDB 4.x 落檔不含頂層 _default 鍵（只有 {table: {id: doc}}）， + 舊判定「無 _default 即 legacy」會把正常 TinyDB 檔整體当成單一 legacy + 條目重寫（反覆包層、全檔損毀）。改以結構判定：頂層有任一 + {數字 id: doc} 形式的 table 即視為 TinyDB 檔、不做遷移。""" if not os.path.exists(path): return try: @@ -68,9 +80,12 @@ def _migrate_legacy_metadata(path): data = json.load(f) except (json.JSONDecodeError, OSError): return - legacy = {k: v for k, v in data.items() - if isinstance(v, dict)} if isinstance(data, dict) else {} - if not (isinstance(data, dict) and '_default' not in data and legacy): + if not isinstance(data, dict) or not data: + return + if any(_looks_like_tinydb_table(v) for v in data.values()): + return + legacy = {k: v for k, v in data.items() if isinstance(v, dict)} + if not legacy: return bak = f'{path}.legacy-{time.strftime("%Y%m%d_%H%M%S")}' os.replace(path, bak) @@ -182,8 +197,14 @@ def _write_rotated_level(vol_dir, level, smd_path, mask_path, roi_path): （1=VBODY、2=棘突、3=other bone、0=background）。 平面參數經 R 剛性旋轉並換元到輸出 grid 的局部座標。 + 完成後把本 level 的旋轉座標（template / R / center / start / size）補寫 + 進 /transform.json（base 部分由 process_single_image 寫入）： + 供分割後把 rotated mask 反向映射回 data_root 原始 CT 座標 + （imaging/transforms.py）。 + 旋轉來源 + 幾何模板：patient 未旋轉的 0.5mm _smd_resampled.nii.gz - （float、裁物件 bbox，SMD < 0.5 = 內部，與 _binary_sdf 同幾何）； + （新世代在 /crop/、舊世代在 / 頂層；float、裁物件 bbox， + SMD < 0.5 = 內部，與 _binary_sdf 同幾何）； 缺失或格式不符時回退 _binary_nn.nii.gz（0.5mm 最近邻 0/1）， 此情形只產出 _binary_nn（+ roi / cortical / 圖 / label）。 mask_path（未旋轉 0.5mm 0/1 遮罩，如 _binary_sdf）只供畫圖 @@ -203,7 +224,7 @@ def _write_rotated_level(vol_dir, level, smd_path, mask_path, roi_path): logger.warning(f'[rotated] {volume_id} {level}: _smd_resampled 不可用 ' f'(value)，改用 _binary_nn') if template is None: - nn_in_path = os.path.join(vol_dir, f'{level}_binary_nn.nii.gz') + nn_in_path = level_file_path(vol_dir, level, 'binary_nn') if not os.path.exists(nn_in_path): logger.warning(f'[rotated] {volume_id} {level}: _smd_resampled / _binary_nn 皆不可用, skip') return @@ -220,7 +241,7 @@ def _write_rotated_level(vol_dir, level, smd_path, mask_path, roi_path): # 圖用遮罩須與模板同幾何（SDF 與 SMD 皆裁 bbox2；NN fallback 只有 # NN 遮罩一致）；0.5mm 0/1，float 檔會被 uint8 讀取成全 0 - nn_in_path = os.path.join(vol_dir, f'{level}_binary_nn.nii.gz') + nn_in_path = level_file_path(vol_dir, level, 'binary_nn') if smd_arr is None: mask_path = nn_in_path if os.path.exists(nn_in_path) else mask_path if mask_path is None or not os.path.exists(mask_path): @@ -354,10 +375,29 @@ def _write_rotated_level(vol_dir, level, smd_path, mask_path, roi_path): f'other={int((label_arr == 3).sum())}, ' f'sp_mode={sp_info["mode"]}, vb_mode={vb_info["mode"]})') + # 記錄 rotated 座標鏈（template 選擇 / R / center / start），補進 + # transform.json；rotated disk a -> 未旋轉 template disk: + # W = R^T (a + start - center) + center（template disk 系，standardize 後） + try: + merge_rotated_into_transform(vol_dir, level, { + 'template': 'smd_resampled' if smd_arr is not None else 'binary_nn', + 'R': [[float(v) for v in row] for row in R], + 'center': [float(v) for v in c_xyz], + 'start': [int(v) for v in fstart], + 'size': [int(v) for v in fsize], + }) + logger.info(f'[rotated] transform.json: recorded {level} rotated section ' + f'(template={ "smd_resampled" if smd_arr is not None else "binary_nn" }, ' + f'start={list(fstart)})') + except Exception as e: + logger.warning(f'[rotated] {volume_id} {level}: transform.json 記錄失敗（{e}）；' + f'該 level 需 --rebuild 才能反向映射') + def make_lumbar_post_process(): """每個 volume 處理完後，對其 lumbar level： - 1) 畫「骨頭 + 方向平面」圖（不畫螺絲、不做棘突 / 椎體分割）-> /lumbar/ + 1) 畫「骨頭 + 方向平面」圖（不畫螺絲、不做棘突 / 椎體分割） + -> /crop/_planes.png（未旋轉檔同在 crop/） 2) 計算對齊旋轉，存旋轉後的 smd_resampled / binary_sdf / binary_nn / roi + cortical + 旋轉平面圖（含 VBODY / 棘突著色）+ label map -> /rotated/ @@ -367,23 +407,25 @@ def make_lumbar_post_process(): def _post_process(vol_dir, processed_labels): volume_id = os.path.basename(vol_dir) - lumbar_dir = os.path.join(vol_dir, 'lumbar') + crop_dir = os.path.join(vol_dir, 'crop') for n in processed_labels: level = LABEL_MAP.get(int(n)) if level not in LUMBAR_LEVELS: continue - smd_res_path = os.path.join(vol_dir, f'{level}_smd_resampled.nii.gz') - nn_path = os.path.join(vol_dir, f'{level}_binary_nn.nii.gz') + # 未旋轉輸入：新世代在 crop/、舊世代在 / 頂層 + smd_res_path = level_file_path(vol_dir, level, 'smd_resampled') + nn_path = level_file_path(vol_dir, level, 'binary_nn') if not os.path.exists(smd_res_path) and not os.path.exists(nn_path): continue # 畫圖用 0/1 遮罩：優 _binary_sdf（0.5mm 平滑），缺則 _binary_nn - sdf_path = os.path.join(vol_dir, f'{level}_binary_sdf.nii.gz') + sdf_path = level_file_path(vol_dir, level, 'binary_sdf') mask_path = sdf_path if os.path.exists(sdf_path) else nn_path - roi_path = os.path.join(vol_dir, f'{level}_roi.nii.gz') + roi_path = level_file_path(vol_dir, level, 'roi') # 1) 原始（未旋轉）planes 圖；皮質著色由未旋轉 CT + mask 現算 - #（未旋轉 _cortical 不再存檔） - output_path = os.path.join(lumbar_dir, f'{level}_planes.png') + #（未旋轉 _cortical 不再存檔）；圖與未旋轉檔同在 crop/ + os.makedirs(crop_dir, exist_ok=True) + output_path = os.path.join(crop_dir, f'{level}_planes.png') if os.path.exists(mask_path): fig_cortical = None if os.path.exists(roi_path): @@ -396,7 +438,7 @@ def make_lumbar_post_process(): path = render_bone_figure(volume_id, level, mask_path, fig_cortical, planes_only=True, output_path=output_path) if path is not None: - logger.info(f'[lumbar] saved {path}') + logger.info(f'[crop] saved {path}') # 2) 旋轉對齊：rotated/ 的 smd_resampled + binary_sdf + binary_nn # + roi + cortical + planes 圖（含 VBODY / 棘突著色）+ label @@ -443,7 +485,7 @@ def main(): 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) + metadata_cache=metadata_db, crop_subdir=True) if __name__ == '__main__': main() \ No newline at end of file diff --git a/xfr_reprocess_ap.py b/xfr_reprocess_ap.py index a57fd48..4b2b4cc 100644 --- a/xfr_reprocess_ap.py +++ b/xfr_reprocess_ap.py @@ -28,17 +28,19 @@ import SimpleITK as sitk sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) from imaging.orientation import anterior_y_side +from imaging.transforms import level_file_path LEVELS = ('L1', 'L2', 'L3', 'L4', 'L5', 'L6') -MASK_SUFFIXES = ('_binary_sdf.nii.gz', '_binary_nn.nii.gz', '_binary.nii.gz') +MASK_SOURCES = ('binary_sdf', 'binary_nn', 'binary') 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}') + # 未旋轉遮罩：新世代在 crop/ 子資料夾、舊世代在頂層 + for src in MASK_SOURCES: + p = level_file_path(vol_dir, lvl, src) if os.path.exists(p): m = sitk.GetArrayFromImage(sitk.ReadImage(p, sitk.sitkUInt8)) per[lvl] = anterior_y_side(m)