{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.12.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceType":"competition","sourceId":117682,"databundleVersionId":15062069,"isSourceIdPinned":false},{"sourceType":"datasetVersion","sourceId":14614416,"datasetId":9335004,"databundleVersionId":15451161},{"sourceType":"datasetVersion","sourceId":14290662,"datasetId":9121922,"databundleVersionId":15093407},{"sourceType":"datasetVersion","sourceId":15173847,"datasetId":9716518,"databundleVersionId":16066249},{"sourceType":"datasetVersion","sourceId":14483320,"datasetId":9250691,"databundleVersionId":15307051},{"sourceType":"datasetVersion","sourceId":15173876,"datasetId":9716537,"databundleVersionId":16066279},{"sourceType":"datasetVersion","sourceId":14611845,"datasetId":9333366,"databundleVersionId":15448370},{"sourceType":"modelInstanceVersion","sourceId":723518,"databundleVersionId":15364551,"modelInstanceId":550541,"modelId":563168,"isSourceIdPinned":false},{"sourceType":"modelInstanceVersion","sourceId":719284,"databundleVersionId":15315523,"modelInstanceId":547016,"modelId":559836,"isSourceIdPinned":false},{"sourceType":"modelInstanceVersion","sourceId":733533,"databundleVersionId":15484876,"modelInstanceId":559011,"modelId":571599,"isSourceIdPinned":false},{"sourceType":"modelInstanceVersion","sourceId":728047,"databundleVersionId":15416012,"modelInstanceId":554345,"modelId":566907,"isSourceIdPinned":false},{"sourceType":"modelInstanceVersion","sourceId":729781,"databundleVersionId":15439544,"modelInstanceId":555833,"modelId":568387,"isSourceIdPinned":false},{"sourceType":"modelInstanceVersion","sourceId":730460,"databundleVersionId":15447480,"modelInstanceId":556409,"modelId":568973,"isSourceIdPinned":false},{"sourceType":"modelInstanceVersion","sourceId":785389,"databundleVersionId":16066177,"modelInstanceId":599290,"modelId":611549,"isSourceIdPinned":false},{"sourceType":"modelInstanceVersion","sourceId":750427,"databundleVersionId":15689394,"modelInstanceId":573182,"modelId":585527,"isSourceIdPinned":false},{"sourceType":"modelInstanceVersion","sourceId":745440,"databundleVersionId":15628003,"modelInstanceId":569096,"modelId":581419,"isSourceIdPinned":false},{"sourceType":"modelInstanceVersion","sourceId":719283,"databundleVersionId":15315506,"modelInstanceId":547015,"modelId":559835,"isSourceIdPinned":false},{"sourceType":"modelInstanceVersion","sourceId":732636,"databundleVersionId":15474274,"modelInstanceId":558311,"modelId":570878,"isSourceIdPinned":false}],"dockerImageVersionId":31236,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# !pip download nnunetv2 -d /kaggle/working/nnunetv2_offline","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Install Packages","metadata":{}},{"cell_type":"code","source":"import os\nimport shutil\n\n# 1. 設定路徑\ninput_path = '/kaggle/input/nnu-net-install/nnunetv2_offline.tar.gz.raw' # 改成你的 Dataset 路徑\nworking_path = '/kaggle/working/nnunetv2_offline.tar.gz'\n\n# 2. 將檔案從唯讀的 input 複製到可讀寫的 working\n# 注意：這會佔用 4GB 的 working 空間\nshutil.copy(input_path, working_path)\n\n# 3. 解壓縮\n# -x: 解壓, -z: gzip, -f: 指定檔案\n# -C: 指定解壓到的目錄\n!mkdir -p /kaggle/working/nnunet_install\n!tar -xzf {working_path} -C /kaggle/working/nnunet_install\n\n# 4. 刪除剛才複製過來的壓縮檔以節省空間\nos.remove(working_path)\n\nprint(\"解壓完成！現在檔案位於 /kaggle/working/nnunet_install，且具有讀寫權限。\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-13T15:29:43.479113Z","iopub.execute_input":"2026-01-13T15:29:43.479317Z","iopub.status.idle":"2026-01-13T15:31:44.153594Z","shell.execute_reply.started":"2026-01-13T15:29:43.479290Z","shell.execute_reply":"2026-01-13T15:31:44.152632Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!pip install --no-index --find-links=/kaggle/working/nnunet_install/nnunetv2_offline nnunetv2","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-13T15:32:03.529137Z","iopub.execute_input":"2026-01-13T15:32:03.530182Z","iopub.status.idle":"2026-01-13T15:32:25.152002Z","shell.execute_reply.started":"2026-01-13T15:32:03.530141Z","shell.execute_reply":"2026-01-13T15:32:25.151149Z"},"collapsed":true,"jupyter":{"outputs_hidden":true}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!pip install --no-index --find-links=/kaggle/input/edt-install edt","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# !python --version","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-13T15:33:44.205719Z","iopub.execute_input":"2026-01-13T15:33:44.206568Z","iopub.status.idle":"2026-01-13T15:33:44.325854Z","shell.execute_reply.started":"2026-01-13T15:33:44.206529Z","shell.execute_reply":"2026-01-13T15:33:44.325262Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# !pip list","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-13T15:33:55.477444Z","iopub.execute_input":"2026-01-13T15:33:55.477731Z","iopub.status.idle":"2026-01-13T15:33:58.551664Z","shell.execute_reply.started":"2026-01-13T15:33:55.477701Z","shell.execute_reply":"2026-01-13T15:33:58.550758Z"},"collapsed":true,"jupyter":{"outputs_hidden":true}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Set variables","metadata":{}},{"cell_type":"code","source":"import os\n\nos.environ['nnUNet_results'] = \"/kaggle/input/models/ggayoayogg/villa-musgd-msr-3d-160-group32-gelu-4000/pytorch/default/1/Villa_MuSGD_MSR_3d_160_Group32_GELU_4000/nnUNet_results\"\nos.environ['nnUNet_raw'] = \"/kaggle/working/nnUNet_raw\"\nos.environ['nnUNet_preprocessed'] = \"/kaggle/working/nnUNet_preprocessed\"\n\n# 檢查是否設定成功\nprint(os.getenv('nnUNet_raw'))","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2026-01-13T15:34:24.197693Z","iopub.execute_input":"2026-01-13T15:34:24.198492Z","iopub.status.idle":"2026-01-13T15:34:24.202904Z","shell.execute_reply.started":"2026-01-13T15:34:24.198452Z","shell.execute_reply":"2026-01-13T15:34:24.202235Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nimport shutil\n\n# Kaggle Python 環境變數設定\nnnUNet_raw = os.environ['nnUNet_raw']\n\n# 建立目錄並複製所有 test images\nimages_ts_dir = os.path.join(nnUNet_raw, 'Dataset501_Vesuvius3D', 'imagesTs')\nos.makedirs(images_ts_dir, exist_ok=True)\n\ntest_dir = '/kaggle/input/vesuvius-challenge-surface-detection/test_images'\nfor filename in os.listdir(test_dir):\n    if filename.endswith('.tif'):\n        src = os.path.join(test_dir, filename)\n        base = os.path.splitext(filename)[0]\n        dst = os.path.join(images_ts_dir, f'{base}_0000.tif')\n        shutil.copy(src, dst)\n        print(f'Copied: {base}_0000.tif')\n\nprint('All test images copied successfully!')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile gpu_worker.py\nimport sys\nimport os\nimport torch\nimport numpy as np\n\n# 1. 暴力清除：如果記憶體中有任何 nnunetv2 相關的模組，全部刪除\nmodules_to_kill = [k for k in sys.modules if k.startswith('nnunetv2')]\nfor k in modules_to_kill:\n    del sys.modules[k]\n\n# ==========================================\n# 1. 確保子進程能找到 nnUNet 與自定義路徑\n# ==========================================\ncustom_path = \"/kaggle/input/datasets/ggayoayogg/villa-nnunet-custom-ayogg\"\nif custom_path not in sys.path:\n    sys.path.insert(0, custom_path)\n\nimport nnunetv2\nfrom nnunetv2.inference.predict_from_raw_data import nnUNetPredictor\nfrom nnunetv2.imageio.simpleitk_reader_writer import SimpleITKIO\nfrom batchgenerators.utilities.file_and_folder_operations import join\n\n# ==========================================\n# 2. 定義路徑配置 (必須與主程式一致)\n# ==========================================\nINPUT_FOLDER = \"/kaggle/working/nnUNet_raw/Dataset501_Vesuvius3D/imagesTs\"\nOUTPUT_FOLDER = \"/kaggle/working/nnUNet_upload/3d_fullres_TTA\"\nNNUNET_RESULTS_FOLDER = \"/kaggle/input/models/ggayoayogg/villa-musgd-msr-3d-160-group32-gelu-4000/pytorch/default/1/Villa_MuSGD_MSR_3d_160_Group32_GELU_4000/nnUNet_results\"\nDATASET_NAME = \"Dataset501_Vesuvius3D\"\nTRAINER_NAME = \"MuSGD_MSR\"\nPLANS_NAME = \"nnUNetResEncUNetMPlans\"\nCONFIG_NAME = \"3d_fullres_160_Group_GELU\"\nFOLD_ID = 1\n\nMODEL_FOLDER = join(\n    NNUNET_RESULTS_FOLDER, \n    DATASET_NAME, \n    f\"{TRAINER_NAME}__{PLANS_NAME}__{CONFIG_NAME}\"\n)\n\n# ==========================================\n# 3. 定義工作函數 (Worker)\n# ==========================================\ndef inference_worker(file_list, gpu_id):\n    if not file_list:\n        return\n\n    print(f\"[GPU {gpu_id}] 啟動... 準備處理 {len(file_list)} 個檔案\")\n    \n    # 初始化預測器\n    predictor = nnUNetPredictor(\n        tile_step_size=0.4,\n        use_gaussian=True,\n        use_mirroring=True,\n        perform_everything_on_device=True,\n        device=torch.device('cuda', gpu_id), \n        verbose=False,\n        verbose_preprocessing=False,\n        allow_tqdm=True\n    )\n\n    predictor.initialize_from_trained_model_folder(\n        MODEL_FOLDER,\n        use_folds=(FOLD_ID,),\n        checkpoint_name='checkpoint_final.pth',\n    )\n    \n    rw = SimpleITKIO()\n\n    for filename in file_list:\n        full_input_path = join(INPUT_FOLDER, filename)\n        \n        # 檔名處理: 去除 _0000\n        name, ext = os.path.splitext(filename)\n        clean_name = name.replace(\"_0000\", \"\") \n        output_filename = clean_name + ext\n        full_output_path = join(OUTPUT_FOLDER, output_filename)\n        \n        print(f\"[GPU {gpu_id}] 正在推理: {clean_name}\")\n\n        try:\n            # 修改處：save_probabilities 設為 False，直接獲取分割結果\n            result = predictor.predict_from_files(\n                list_of_lists_or_source_folder=[[full_input_path]],\n                output_folder_or_list_of_truncated_output_files=None,\n                save_probabilities=False\n            )\n            \n            # 修改處：直接取用預測出的 Label (通常在 index 0)\n            final_seg = result[0].astype(np.uint8)\n\n            # 確保維度正確 (如果是 2D 則補上通道)\n            if final_seg.ndim == 2:\n                final_seg = final_seg[None, ...]\n                \n            dummy_props = {\n                'sitk_stuff': {\n                    'spacing': (1.0, 1.0, 1.0), \n                    'origin': (0.0, 0.0, 0.0),\n                    'direction': (1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)\n                },\n                'spacing': [1.0, 1.0, 1.0]\n            }\n\n            rw.write_seg(final_seg, full_output_path, dummy_props)\n            print(f\"[GPU {gpu_id}] 完成儲存: {output_filename}\")\n            \n        except Exception as e:\n            print(f\"[GPU {gpu_id}] 發生錯誤 ({filename}): {e}\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import sys\nimport os\nimport torch.multiprocessing as mp\nfrom batchgenerators.utilities.file_and_folder_operations import subfiles, maybe_mkdir_p\n\n# 1. 將 Kaggle 工作目錄加入系統路徑，確保能 import 剛剛寫的檔案\nsys.path.append(\"/kaggle/working\")\n\n# 2. 清除舊的緩存 (如果有的話)，避免修改後沒更新\nif 'gpu_worker' in sys.modules:\n    del sys.modules['gpu_worker']\n\n# 3. 從剛剛建立的 gpu_worker.py 匯入函數與變數\ntry:\n    from gpu_worker import inference_worker, INPUT_FOLDER, OUTPUT_FOLDER, MODEL_FOLDER\nexcept ImportError as e:\n    print(\"錯誤：找不到 gpu_worker.py。請確認您已執行上一個 '%%writefile' 的 Cell。\")\n    raise e\n\ndef run_inference_manager():\n    if not os.path.exists(MODEL_FOLDER):\n        raise FileNotFoundError(f\"找不到模型資料夾: {MODEL_FOLDER}\")\n\n    # 搜尋檔案\n    if not os.path.exists(INPUT_FOLDER):\n        print(f\"警告：輸入資料夾不存在 {INPUT_FOLDER}\")\n        return\n        \n    files = subfiles(INPUT_FOLDER, suffix='.tif', join=False)\n    files.sort()\n    maybe_mkdir_p(OUTPUT_FOLDER)\n\n    num_files = len(files)\n    print(f\"=== 任務開始 ===\")\n    print(f\"發現檔案數: {num_files}\")\n\n    if num_files == 0:\n        print(\"無檔案需處理。\")\n        return\n\n    # === 核心判斷邏輯 ===\n    if num_files == 1:\n        print(\"模式: 單一檔案 -> 使用單卡 (GPU 0)\")\n        # 直接呼叫，不需開新 Process\n        inference_worker(files, 0)\n    else:\n        print(f\"模式: 多檔案 -> 啟動雙卡並行加速 (GPU 0 & GPU 1)\")\n        \n        # 設定 spawn 模式 (CUDA 必須)\n        try:\n            mp.set_start_method('spawn', force=True)\n        except RuntimeError:\n            pass\n\n        # 拆分檔案\n        mid = num_files // 2\n        files_gpu0 = files[:mid]\n        files_gpu1 = files[mid:]\n\n        # 建立並行任務\n        p1 = mp.Process(target=inference_worker, args=(files_gpu0, 0))\n        p2 = mp.Process(target=inference_worker, args=(files_gpu1, 1))\n\n        p1.start()\n        p2.start()\n        \n        p1.join()\n        p2.join()\n        print(\"所有 GPU 任務已完成。\")\n\nif __name__ == \"__main__\":\n    run_inference_manager()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Post-Processing Function","metadata":{}},{"cell_type":"code","source":"import cv2\nimport math\nfrom numba import njit\nfrom skimage.morphology import skeletonize, medial_axis, reconstruction, disk\nfrom skimage.measure import label\nfrom skimage.draw import line as skimage_line\nfrom scipy.ndimage import convolve, distance_transform_edt, binary_closing, binary_dilation, generate_binary_structure, \\\n    gaussian_filter, binary_opening, binary_hit_or_miss, binary_erosion\nimport numpy as np\nfrom scipy import ndimage as ndi\nfrom scipy.ndimage import binary_fill_holes\nfrom skimage.morphology import remove_small_objects, remove_small_holes\nfrom skimage.draw import line\nfrom skimage.measure import label as measure_label, regionprops\nfrom skimage.segmentation import watershed, find_boundaries\nfrom scipy.spatial.distance import cdist\nfrom typing import List, Tuple\nfrom skimage.graph import route_through_array\n\n# ==========================================\n# 1. 幾何與邏輯計算核心 (從原補線代碼移植)\n# ==========================================\n# 嘗試導入路徑計算模組\ntry:\n    from skimage.graph import route_through_array\nexcept ImportError:\n    from skimage.graph import MCP_Geometric\n\n\n    def route_through_array(cost, start, end):\n        mcp = MCP_Geometric(cost)\n        cumulative_costs, traceback = mcp.find_costs([start], [end])\n        return mcp.traceback(end), 0\n\n\n# =========================================================================\n#  Numba 加速核心 (Geodesic, PCA, Ray-Casting)\n# =========================================================================\n\n@njit(cache=True)\ndef compute_geodesic_distance_numba(mask, start_point):\n    \"\"\"計算 Geodesic Distance Transform\"\"\"\n    d, h, w = mask.shape\n    dist_map = np.full((d, h, w), np.inf, dtype=np.float32)\n    sz, sy, sx = start_point\n\n    # 起點修正 (若起點在背景，搜尋鄰域)\n    if mask[sz, sy, sx] == 0:\n        found = False\n        for dz in range(-1, 2):\n            for dy in range(-1, 2):\n                for dx in range(-1, 2):\n                    nz, ny, nx = sz + dz, sy + dy, sx + dx\n                    if 0 <= nz < d and 0 <= ny < h and 0 <= nx < w:\n                        if mask[nz, ny, nx] != 0:\n                            sz, sy, sx = nz, ny, nx\n                            found = True\n                            break\n                if found: break\n            if found: break\n        if not found: return np.zeros((d, h, w), dtype=np.float32)\n\n    # BFS 初始化\n    max_queue_len = d * h * w\n    queue_z = np.empty(max_queue_len, dtype=np.int32)\n    queue_y = np.empty(max_queue_len, dtype=np.int32)\n    queue_x = np.empty(max_queue_len, dtype=np.int32)\n    head, tail = 0, 0\n\n    queue_z[tail], queue_y[tail], queue_x[tail] = sz, sy, sx\n    tail += 1\n    dist_map[sz, sy, sx] = 0.0\n\n    dz_off = np.array([0, 0, 0, 0, 1, -1], dtype=np.int8)\n    dy_off = np.array([0, 0, 1, -1, 0, 0], dtype=np.int8)\n    dx_off = np.array([1, -1, 0, 0, 0, 0], dtype=np.int8)\n\n    while head < tail:\n        cz, cy, cx = queue_z[head], queue_y[head], queue_x[head]\n        head += 1\n        next_dist = dist_map[cz, cy, cx] + 1.0\n\n        for i in range(6):\n            nz, ny, nx = cz + int(dz_off[i]), cy + int(dy_off[i]), cx + int(dx_off[i])\n            if 0 <= nz < d and 0 <= ny < h and 0 <= nx < w:\n                if mask[nz, ny, nx] != 0 and dist_map[nz, ny, nx] == np.inf:\n                    dist_map[nz, ny, nx] = next_dist\n                    queue_z[tail], queue_y[tail], queue_x[tail] = nz, ny, nx\n                    tail += 1\n    return dist_map\n\n\n@njit(fastmath=True)\ndef compute_pca_field_numba(labels, radius, step=1):\n    h, w = labels.shape\n    norm_y = np.zeros((h, w), dtype=np.float32)\n    norm_x = np.zeros((h, w), dtype=np.float32)\n\n    for r in range(0, h, step):\n        for c in range(0, w, step):\n            current_id = labels[r, c]\n            if current_id == 0: continue\n            sum_xx, sum_yy, sum_xy, count = 0.0, 0.0, 0.0, 0.0\n            r_min, r_max = max(0, r - radius), min(h, r + radius + 1)\n            c_min, c_max = max(0, c - radius), min(w, c + radius + 1)\n            for nr in range(r_min, r_max):\n                for nc in range(c_min, c_max):\n                    if labels[nr, nc] == current_id:\n                        dy, dx = nr - r, nc - c\n                        sum_xx += dx * dx\n                        sum_yy += dy * dy\n                        sum_xy += dx * dy\n                        count += 1.0\n            if count < 3: continue\n            A, B, C = sum_yy, sum_xy, sum_xx\n            tangent_angle = 0.5 * np.arctan2(2 * B, C - A)\n            normal_angle = tangent_angle + np.pi / 2.0\n            norm_y[r, c] = np.sin(normal_angle)\n            norm_x[r, c] = np.cos(normal_angle)\n    return norm_y, norm_x\n\n\n@njit(fastmath=True)\ndef find_first_collision_pair(y_coords, x_coords, pca_ny, pca_nx, mask, labels, ray_len):\n    h, w = mask.shape\n    num_points = len(y_coords)\n    for i in range(num_points):\n        r, c = y_coords[i], x_coords[i]\n        ny, nx = pca_ny[r, c], pca_nx[r, c]\n        if ny == 0.0 and nx == 0.0: continue\n        source_id = labels[r, c]\n        for sign in (1, -1):\n            has_left_self = False\n            for t in range(1, ray_len):\n                cy = int(r + sign * t * ny + 0.5)\n                cx = int(c + sign * t * nx + 0.5)\n                if cy < 0 or cy >= h or cx < 0 or cx >= w: break\n                pixel_val = mask[cy, cx]\n                target_id = labels[cy, cx]\n                if not has_left_self:\n                    if pixel_val == 0:\n                        has_left_self = True\n                    elif target_id != source_id:\n                        return True, r, c, cy, cx\n                else:\n                    if pixel_val > 0:\n                        if target_id != source_id:\n                            return True, r, c, cy, cx\n                        else:\n                            break\n    return False, 0, 0, 0, 0\n\n\n# =========================================================================\n#  Python 輔助 (Seeds, Path)\n# =========================================================================\n\ndef get_split_seeds(mask_2d, ray_len=64, pca_radius=8, pca_step=3):\n    mask_uint8 = (mask_2d > 0).astype(np.uint8)\n    y_coords, x_coords = np.where(mask_uint8 > 0)\n    if len(y_coords) < 10: return None\n    valid_indices = (y_coords % pca_step == 0) & (x_coords % pca_step == 0)\n    y_sampled = y_coords[valid_indices]\n    x_sampled = x_coords[valid_indices]\n    if len(y_sampled) == 0: return None\n    _, labels_map = cv2.connectedComponents(mask_uint8, connectivity=8)\n    labels_map = labels_map.astype(np.int32)\n    pca_ny, pca_nx = compute_pca_field_numba(labels_map, pca_radius, step=pca_step)\n    found, r1, c1, r2, c2 = find_first_collision_pair(\n        y_sampled, x_sampled, pca_ny, pca_nx, mask_uint8, labels_map, ray_len\n    )\n    if found: return (r1, c1), (r2, c2)\n    return None\n\n\ndef analyze_shortest_path(mask_3d, start_pt, end_pt):\n    dist = distance_transform_edt(mask_3d)\n    cost_map = np.max(dist) - dist + 1\n    cost_map[mask_3d == 0] = np.inf\n    try:\n        indices, _ = route_through_array(cost_map, start_pt, end_pt)\n        path_indices = np.array(indices)\n        if len(path_indices) == 0: return None, None\n        midpoint = tuple(map(int, path_indices[len(path_indices) // 2]))\n        return path_indices, midpoint\n    except:\n        return None, None\n\n\n# =========================================================================\n#  主函數: split_paper\n# =========================================================================\n\ndef get_split_paper(mask_3d, ray_len=64, max_iter=5, cleanup_iter=2):\n    \"\"\"\n    對 3D Binary Mask 執行迭代式幾何分割與邊界清理。\n\n    Args:\n        mask_3d (np.ndarray): 輸入的 3D Mask (boolean 或 uint8)。\n        ray_len (int): 射線偵測長度，用於尋找凹陷處。\n        max_iter (int): 最大迭代次數 (重複分割的輪數)。\n        cleanup_iter (int): 邊界清理強度 (侵蝕層數)，數值越大縫隙越寬。\n\n    Returns:\n        np.ndarray: 分割後的 3D Binary Mask (uint8, 0=背景/縫隙, 1=前景)。\n    \"\"\"\n\n    # 確保輸入為 Instance Labels (方便內部運算)\n    # 如果輸入是 Binary，先做一次連通分量標記\n    if mask_3d.dtype == bool or np.max(mask_3d) == 1:\n        refined_labels = measure_label(mask_3d)\n    else:\n        # 假設已經是 label map，但為了安全起見，重新標記以確保 ID 連續性\n        refined_labels = measure_label(mask_3d > 0)\n\n    refined_labels = refined_labels.astype(np.int32)\n    next_new_label = refined_labels.max() + 1\n\n    print(f\">> 啟動 split_paper: Ray={ray_len}, Iter={max_iter}, Cleanup={cleanup_iter}\")\n\n    # 迭代分割迴圈\n    for iteration in range(max_iter):\n        unique_ids = np.unique(refined_labels)\n        unique_ids = unique_ids[unique_ids != 0]\n        splits_this_round = 0\n\n        # 針對每個物件檢查\n        for obj_id in unique_ids:\n            obj_mask = (refined_labels == obj_id)\n            slices = np.where(obj_mask)\n            if len(slices[0]) == 0: continue\n\n            # ROI 提取\n            z_min, z_max = np.min(slices[0]), np.max(slices[0])\n            y_min, y_max = np.min(slices[1]), np.max(slices[1])\n            x_min, x_max = np.min(slices[2]), np.max(slices[2])\n\n            roi_mask = obj_mask[z_min:z_max + 1, y_min:y_max + 1, x_min:x_max + 1]\n            d, h, w = roi_mask.shape\n\n            # 3點取樣 (上中下)\n            target_z_indices = sorted(list(set([0, d // 2, d - 1])))\n            valid_midpoints = []\n            seed_pairs = {}\n\n            # A. 尋找種子點與路徑中點\n            for z_local in target_z_indices:\n                slice_2d = roi_mask[z_local, :, :]\n                if np.sum(slice_2d) < 50: continue\n\n                seeds = get_split_seeds(slice_2d, ray_len=ray_len)\n                if seeds:\n                    (y1, x1), (y2, x2) = seeds\n                    p1 = (z_local, y1, x1)\n                    p2 = (z_local, y2, x2)\n                    _, midpoint = analyze_shortest_path(roi_mask, p1, p2)\n                    if midpoint is not None:\n                        valid_midpoints.append(midpoint)\n                        seed_pairs[z_local] = (p1, p2)\n\n            # B. 若有足夠證據，進行 Geodesic Distance Fusion\n            if len(valid_midpoints) > 0:\n                fused_dist_map = np.full(roi_mask.shape, np.inf, dtype=np.float32)\n                roi_mask_int = roi_mask.astype(np.uint8)\n\n                for midpoint in valid_midpoints:\n                    start_pt = (int(midpoint[0]), int(midpoint[1]), int(midpoint[2]))\n                    d_map = compute_geodesic_distance_numba(roi_mask_int, start_pt)\n                    fused_dist_map = np.minimum(fused_dist_map, d_map)\n\n                fused_dist_map[fused_dist_map == np.inf] = 0\n\n                # C. 分水嶺分割\n                chosen_seeds = None\n                if (d // 2) in seed_pairs:\n                    chosen_seeds = seed_pairs[d // 2]\n                elif len(seed_pairs) > 0:\n                    chosen_seeds = seed_pairs[list(seed_pairs.keys())[0]]\n\n                if chosen_seeds:\n                    p1, p2 = chosen_seeds\n                    basin_map = -1.0 * fused_dist_map\n                    markers = np.zeros_like(roi_mask, dtype=np.int32)\n                    markers[p1] = 1\n                    markers[p2] = 2\n\n                    labels_ws = watershed(basin_map, markers, mask=roi_mask)\n\n                    if np.max(labels_ws) == 2:\n                        # 標記新的 ID\n                        split_mask_2 = (labels_ws == 2)\n                        roi_view = refined_labels[z_min:z_max + 1, y_min:y_max + 1, x_min:x_max + 1]\n                        roi_view[split_mask_2] = next_new_label\n\n                        # D. 邊界清理 (Boundary Cleanup / Erosion)\n                        mask_old = (roi_view == obj_id)\n                        mask_new = (roi_view == next_new_label)\n                        struct_26 = np.ones((3, 3, 3), dtype=bool)\n\n                        # 雙向擴張找交界\n                        neighbors_have_old = binary_dilation(mask_old, structure=struct_26, iterations=cleanup_iter)\n                        boundary_pixels_in_new = mask_new & neighbors_have_old\n\n                        neighbors_have_new = binary_dilation(mask_new, structure=struct_26, iterations=cleanup_iter)\n                        boundary_pixels_in_old = mask_old & neighbors_have_new\n\n                        # 將交界處設為背景 (0)\n                        roi_view[boundary_pixels_in_new] = 0\n                        roi_view[boundary_pixels_in_old] = 0\n\n                        next_new_label += 1\n                        splits_this_round += 1\n\n        print(f\"   Iter {iteration + 1}/{max_iter}: 分割了 {splits_this_round} 處。\")\n        if splits_this_round == 0:\n            break\n\n    # 轉回 Binary Mask 輸出 (uint8)\n    final_binary = (refined_labels > 0).astype(np.uint8)\n    return final_binary\n\n\ndef get_outward_vector(skel, ep_xy, pca_vec):\n    \"\"\"計算骨架端點的朝外向量\"\"\"\n    start_x, start_y = int(ep_xy[0]), int(ep_xy[1])\n    h, w = skel.shape\n    STEPS = 5\n    visited = set()\n    visited.add((start_y, start_x))\n    curr_x, curr_y = start_x, start_y\n    for _ in range(STEPS):\n        found_next = False\n        for dy in [-1, 0, 1]:\n            for dx in [-1, 0, 1]:\n                if dx == 0 and dy == 0: continue\n                ny, nx = curr_y + dy, curr_x + dx\n                if ny < 0 or ny >= h or nx < 0 or nx >= w: continue\n                if skel[ny, nx] > 0 and (ny, nx) not in visited:\n                    curr_x, curr_y = nx, ny\n                    visited.add((ny, nx))\n                    found_next = True\n                    break\n            if found_next: break\n        if not found_next: break\n    ref_vec_x, ref_vec_y = start_x - curr_x, start_y - curr_y\n    if ref_vec_x == 0 and ref_vec_y == 0: return pca_vec\n    if ref_vec_x * pca_vec[0] + ref_vec_y * pca_vec[1] < 0: return -pca_vec\n    return pca_vec\n\n\ndef check_path_blocked(p1, p2, mask_img):\n    \"\"\"檢查兩點之間的路徑是否被障礙物阻擋\"\"\"\n    x1, y1 = int(p1[0]), int(p1[1])\n    x2, y2 = int(p2[0]), int(p2[1])\n    rr, cc = skimage_line(y1, x1, y2, x2)\n    valid_mask = (rr >= 0) & (rr < mask_img.shape[0]) & (cc >= 0) & (cc < mask_img.shape[1])\n    rr, cc = rr[valid_mask], cc[valid_mask]\n    if len(rr) == 0: return False\n    safe_radius_sq = 4 ** 2\n\n    def is_obstacle(px, py):\n        if px < 0 or px >= mask_img.shape[1] or py < 0 or py >= mask_img.shape[0]: return False\n        if mask_img[py, px] == 0: return False\n        d1 = (px - x1) ** 2 + (py - y1) ** 2\n        d2 = (px - x2) ** 2 + (py - y2) ** 2\n        return (d1 > safe_radius_sq) and (d2 > safe_radius_sq)\n\n    for i in range(len(rr)):\n        curr_y, curr_x = rr[i], cc[i]\n        if is_obstacle(curr_x, curr_y): return True\n        if i > 0:\n            prev_y, prev_x = rr[i - 1], cc[i - 1]\n            if abs(curr_x - prev_x) == 1 and abs(curr_y - prev_y) == 1:\n                if is_obstacle(prev_x, curr_y) and is_obstacle(curr_x, prev_y): return True\n    return False\n\n\ndef get_pca_tangent(skeleton, center_yx, radius=16, line_len=16):\n    \"\"\"使用 PCA 計算端點的切線方向\"\"\"\n    #\n    cy, cx = center_yx\n    h, w = skeleton.shape\n    y_min, y_max = max(0, cy - radius), min(h, cy + radius + 1)\n    x_min, x_max = max(0, cx - radius), min(w, cx + radius + 1)\n    roi = skeleton[y_min:y_max, x_min:x_max].copy()\n    if not np.any(roi): return None\n    roi_labels = label(roi > 0, connectivity=2)\n    local_cy, local_cx = cy - y_min, cx - x_min\n    target_label = roi_labels[local_cy, local_cx]\n    if target_label == 0: return None\n    pts_y, pts_x = np.where(roi_labels == target_label)\n    pts_global_y, pts_global_x = pts_y + y_min, pts_x + x_min\n    dist_sq = (pts_global_y - cy) ** 2 + (pts_global_x - cx) ** 2\n    mask = dist_sq <= (radius ** 2)\n    valid_y, valid_x = pts_global_y[mask], pts_global_x[mask]\n    if len(valid_x) < 2: return None\n    data = np.vstack([valid_x, valid_y]).T\n    mean = np.mean(data, axis=0)\n    cov = np.cov((data - mean).T)\n    if np.isnan(cov).any() or np.isinf(cov).any(): return None\n    eig_vals, eig_vecs = np.linalg.eigh(cov)\n    principal_vec = eig_vecs[:, -1]\n    vec_len = np.linalg.norm(principal_vec)\n    if vec_len == 0: return None\n    vx, vy = principal_vec / vec_len\n    pt1 = (int(cx - vx * line_len), int(cy - vy * line_len))\n    pt2 = (int(cx + vx * line_len), int(cy + vy * line_len))\n    return pt1, pt2\n\n\ndef get_endpoints(mask_2d):\n    \"\"\"取得 2D mask 的骨架端點\"\"\"\n    #\n    if not np.any(mask_2d): return [], None\n    skeleton = skeletonize(mask_2d.astype(bool))\n    skeleton_uint8 = skeleton.astype(np.uint8) * 255\n    kernel = np.array([[1, 1, 1], [1, 0, 1], [1, 1, 1]], dtype=np.uint8)\n    neighbors = convolve(skeleton.astype(np.uint8), kernel, mode='constant', cval=0)\n    coords = np.argwhere(skeleton & (neighbors <= 1))\n    if len(coords) == 0: return [], skeleton_uint8\n    raw_eps = [[pt[1], pt[0]] for pt in coords]\n    raw_eps.sort(key=lambda p: (p[1], p[0]))\n    return raw_eps, skeleton_uint8\n\n\ndef get_bezier_points(p0, p3, v0, v3, scale=0.5, num_points=30):\n    \"\"\"生成貝茲曲線點集\"\"\"\n    #\n    dist = np.linalg.norm(p3 - p0)\n    control_dist = dist * scale\n    p1 = p0 + v0 * control_dist\n    p2 = p3 + v3 * control_dist\n    t = np.linspace(0, 1, num_points).reshape(-1, 1)\n    points = (1 - t) ** 3 * p0 + 3 * (1 - t) ** 2 * t * p1 + 3 * (1 - t) * t ** 2 * p2 + t ** 3 * p3\n    return points.astype(np.int32)\n\n\ndef compute_candidate_info(all_eps, skel_img, radius=16):\n    \"\"\"計算所有候選端點的詳細資訊 (位置、向量)\"\"\"\n    ep_info_list = []\n    for i, (x, y) in enumerate(all_eps):\n        t_pts = get_pca_tangent(skel_img, (y, x), radius=radius, line_len=radius)\n        if not t_pts: continue\n        pt1, pt2 = t_pts\n        vx, vy = pt2[0] - pt1[0], pt2[1] - pt1[1]\n        norm = np.sqrt(vx ** 2 + vy ** 2)\n        if norm == 0: continue\n        pca_vec = np.array([vx / norm, vy / norm])\n        out_vec = get_outward_vector(skel_img, (x, y), pca_vec)\n        ep_info_list.append({'id': i, 'xy': np.array([x, y]), 'vec': out_vec, 'tangent_pts': t_pts})\n    return ep_info_list\n\n\ndef dynamic_filter(ep_info_list, skel_img, use_angle=True, use_block=True, use_pair=True):\n    \"\"\"篩選並配對需要連接的端點\"\"\"\n    num_eps = len(ep_info_list)\n    if num_eps < 2: return []\n\n    candidate_links = []\n    ANGLE_WEIGHT = 10.0\n    COS_THRES = np.cos(np.deg2rad(60))\n\n    for i in range(num_eps):\n        for j in range(i + 1, num_eps):\n            p1 = ep_info_list[i]\n            p2 = ep_info_list[j]\n            vec_p1_to_p2 = p2['xy'] - p1['xy']\n            dist = np.linalg.norm(vec_p1_to_p2)\n            if dist == 0: continue\n            dir_1_to_2 = vec_p1_to_p2 / dist\n            dir_2_to_1 = -dir_1_to_2\n            cos_p1 = np.dot(p1['vec'], dir_1_to_2)\n            cos_p2 = np.dot(p2['vec'], dir_2_to_1)\n\n            if use_angle and (cos_p1 < COS_THRES or cos_p2 < COS_THRES): continue\n            if use_block and check_path_blocked(tuple(p1['xy']), tuple(p2['xy']), skel_img): continue\n\n            avg_cos = (cos_p1 + cos_p2) / 2.0\n            score = dist * (1.0 + ANGLE_WEIGHT * (1.0 - min(avg_cos, 1.0)))\n            candidate_links.append((score, i, j))\n\n    candidate_links.sort(key=lambda x: x[0])\n    matched_indices = set()\n    final_connections = []\n\n    if use_pair:\n        for score, idx1, idx2 in candidate_links:\n            if idx1 in matched_indices or idx2 in matched_indices: continue\n            p1, p2 = ep_info_list[idx1], ep_info_list[idx2]\n            # 再次檢查遮擋 (使用原始骨架)\n            if use_block and check_path_blocked(tuple(p1['xy']), tuple(p2['xy']), skel_img): continue\n            matched_indices.add(idx1);\n            matched_indices.add(idx2)\n            final_connections.append((p1, p2))\n    else:\n        # 如果不強制配對，這段邏輯在原代碼中行為不同，這裡保留配對邏輯為主\n        pass\n\n    return final_connections\n\n\ndef execute_line_repair(vol_3d, settings=None):\n    \"\"\"\n    執行 3D 體積的逐層補線操作\n\n    Args:\n        vol_3d: 3D numpy array (Boolean or Label)\n        settings: 字典，包含 'curve_str', 'use_angle' 等參數\n\n    Returns:\n        patched_vol: 補線後的 3D numpy array (uint8, 0 or 255)\n    \"\"\"\n    if settings is None:\n        settings = {\n            'use_angle': True,\n            'use_block': True,\n            'use_pair': True,\n            'trim_ends': True,\n            'curve_str': 50\n        }\n\n    print(\"正在執行斷線修復 (Line Endpoint Repair)...\")\n\n    # 確保輸入是 Label 格式以區分 3D 物體，避免錯誤連接不同物體\n    if vol_3d.dtype == bool:\n        labeled_vol = label(vol_3d)\n    else:\n        # 假設已經是 label 或者 uint8\n        labeled_vol = label(vol_3d > 0)\n\n    # 準備輸出容器 (uint8)\n    patched_vol = (vol_3d > 0).astype(np.uint8) * 255\n    curve_scale = settings['curve_str'] / 100.0\n\n    d, h, w = labeled_vol.shape\n\n    # 逐層處理\n    for z in range(d):\n        slice_labels = labeled_vol[z, :, :]\n        if not np.any(slice_labels): continue\n\n        slice_img = patched_vol[z, :, :]  # 這是要畫圖的畫布\n\n        present_ids = np.unique(slice_labels)\n        present_ids = present_ids[present_ids > 0]  # 排除背景 0\n\n        for pid in present_ids:\n            # 取出單一物件的 Mask\n            obj_mask = (slice_labels == pid)\n\n            # 1. 取得端點與骨架\n            raw_eps, skel_img = get_endpoints(obj_mask)\n            if len(raw_eps) < 2: continue\n\n            # 2. 計算向量資訊\n            ep_info_all = compute_candidate_info(raw_eps, skel_img, radius=16)\n\n            # 3. 篩選與配對\n            candidates = list(ep_info_all)\n            if settings['trim_ends'] and len(candidates) > 2:\n                candidates.sort(key=lambda p: (p['xy'][1], p['xy'][0]))\n                candidates = candidates[1:-1]\n\n            connections = dynamic_filter(\n                candidates, skel_img,\n                use_angle=settings['use_angle'],\n                use_block=settings['use_block'],\n                use_pair=settings['use_pair']\n            )\n\n            # 4. 畫貝茲曲線\n            # for (p1_obj, p2_obj) in connections:\n            #     curve_pts = get_bezier_points(\n            #         p1_obj['xy'], p2_obj['xy'],\n            #         p1_obj['vec'], p2_obj['vec'],\n            #         scale=curve_scale\n            #     )\n            #     # 直接畫在該層的結果圖上\n            #     cv2.polylines(slice_img, [curve_pts], False, 255, 2, cv2.LINE_AA)\n\n            # 4. 畫直線\n            for (p1_obj, p2_obj) in connections:\n                # 確保座標格式為整數 Tuple (x, y)\n                # 如果 p1_obj['xy'] 已經是整數 Tuple，可以直接用 p1_obj['xy']\n                pt1 = tuple(map(int, p1_obj['xy']))\n                pt2 = tuple(map(int, p2_obj['xy']))\n\n                # 使用 cv2.line 直接畫直線\n                # 參數：圖片, 起點, 終點, 顏色(255), 粗細(2), 反鋸齒(LINE_AA)\n                cv2.line(slice_img, pt1, pt2, 255, 2, cv2.LINE_AA)\n\n        patched_vol[z, :, :] = slice_img\n\n    return patched_vol\n\n\ndef line_repair_msk(vol_bool: np.ndarray, pass_iters: int = 1, radius: float = 1.5, axis: int = 0, max_dist: int = 30,\n                    max_angle_deg: int = 30, safe_distance: int = 4, border: int = 8):\n    \"\"\"\n    核心修補邏輯：檢查 3D 連通物件是否在 2D 切片上斷裂並修復。\n    加入 Dijkstra 成本圖演算法以避開原 Mask，降低沾黏機率。\n    \"\"\"\n\n    def get_2d_endpoints(skeleton):\n        # 建立一個中間為 0，周圍為 1 的 kernel\n        kernel = np.array([[1, 1, 1],\n                           [1, 0, 1],\n                           [1, 1, 1]])\n        neighbor_count = convolve(skeleton.astype(int), kernel, mode='constant', cval=0)\n\n        # 骨架點且鄰居只有 1 個\n        endpoints = np.argwhere((skeleton == 1) & (neighbor_count == 1))\n        return endpoints\n\n    def trace_back_and_get_vector_2d(skeleton, endpoint, steps=3):\n        \"\"\"\n        從端點沿著骨架往回走指定的步數，計算生長方向向量。\n        \"\"\"\n        current_pt = endpoint\n        visited = {tuple(endpoint)}\n\n        for _ in range(steps):\n            y, x = current_pt\n            y_min, y_max = max(0, y - 1), min(skeleton.shape[0], y + 2)\n            x_min, x_max = max(0, x - 1), min(skeleton.shape[1], x + 2)\n\n            neighborhood = skeleton[y_min:y_max, x_min:x_max]\n            neighbors = np.argwhere(neighborhood == 1)\n            neighbors += np.array([y_min, x_min])\n\n            next_pt = None\n            for n in neighbors:\n                if tuple(n) not in visited:\n                    next_pt = n\n                    break\n\n            if next_pt is None:\n                return None\n\n            visited.add(tuple(next_pt))\n            current_pt = next_pt\n\n        # 向量方向：從回溯點指向端點 (即未來的生長趨勢)\n        vec = endpoint - current_pt\n        norm = np.linalg.norm(vec)\n        if norm == 0:\n            return np.array([0.0, 0.0])\n        return vec / norm\n\n    def match_endpoints_2d(endpoints, vectors, label_slice, max_dist=30, max_angle_deg=45):\n        \"\"\"\n        根據距離和方向夾角配對端點\n        \"\"\"\n        matched_pairs = []\n        used_indices = set()\n        max_angle_rad = math.radians(max_angle_deg)\n\n        for i in range(len(endpoints)):\n            if i in used_indices: continue\n\n            best_match = -1\n            min_dist = max_dist\n\n            for j in range(len(endpoints)):\n                if i == j or j in used_indices: continue\n\n                p_A, v_A = endpoints[i], vectors[i]\n                p_B, v_B = endpoints[j], vectors[j]\n\n                label_A = label_slice[p_A[0], p_A[1]]\n                label_B = label_slice[p_B[0], p_B[1]]\n\n                if label_A != label_B:\n                    continue\n\n                # 1. 檢查距離\n                dist = np.linalg.norm(p_A - p_B)\n                if dist > max_dist: continue\n\n                # 2. 檢查方向性\n                vec_AB = (p_B - p_A) / dist\n                vec_BA = -vec_AB\n\n                # 確保向量長度不為 0 才計算夾角\n                if np.linalg.norm(v_A) == 0 or np.linalg.norm(v_B) == 0: continue\n\n                # v_A 應該要指向 B，v_B 應該要指向 A\n                # cos_A = np.clip(np.dot(v_A, vec_AB), -1.0, 1.0)\n                # cos_B = np.clip(np.dot(v_B, vec_BA), -1.0, 1.0)\n                #\n                # if math.acos(cos_A) < max_angle_rad and math.acos(cos_B) < max_angle_rad:\n                #     if dist < min_dist:  # 尋找距離最近的合格者\n                #         min_dist = dist\n                #         best_match = j\n\n                dot_A_forward = np.dot(v_A, vec_AB)\n                dot_B_forward = np.dot(v_B, vec_BA)\n\n                # 2. 檢查兩個端點本身的生長向量是否「相向」(例如夾角大於 135 度)\n                # v_B 應該要和 v_A 大致反向，所以 v_A 和 -v_B 應該要大致同向\n                cos_dirs = np.clip(np.dot(v_A, -v_B), -1.0, 1.0)\n                dirs_aligned = math.acos(cos_dirs) < max_angle_rad  # 這裡的 max_angle_rad 可以設為 45 度 (即允許 45 度的方向誤差)\n\n                if dot_A_forward > 0 and dot_B_forward > 0 and dirs_aligned:\n                    if dist < min_dist:\n                        min_dist = dist\n                        best_match = j\n\n            if best_match != -1:\n                matched_pairs.append({\n                    'pA': endpoints[i], 'pB': endpoints[best_match],\n                    'vA': vectors[i], 'vB': vectors[best_match]\n                })\n                used_indices.add(i)\n                used_indices.add(best_match)\n\n        return matched_pairs\n\n    def connect_skeleton_holes_2d(skeleton, label_slice, max_dist=30, max_angle_deg=45, trace_steps=5):\n        \"\"\"\n        主函數：輸入 2D 骨架矩陣，輸出修補好的骨架矩陣\n        \"\"\"\n        # 複製一份準備畫線用\n        repaired_skeleton = skeleton.copy()\n\n        # 1. 抓取端點\n        raw_endpoints = get_2d_endpoints(skeleton)\n        if len(raw_endpoints) < 2:\n            return repaired_skeleton\n\n        # 2. 【修改處】：直接先拿掉太短的骨架，不參與配對\n        valid_endpoints = []\n        valid_vectors = []\n        for ep in raw_endpoints:\n            # 這裡的 trace_steps 可以依需求設為 3 或 5\n            vec = trace_back_and_get_vector_2d(skeleton, ep, steps=trace_steps)\n\n            # 如果 vec 不是 None，代表這條骨架夠長，才允許加入配對池\n            if vec is not None:\n                valid_endpoints.append(ep)\n                valid_vectors.append(vec)\n\n        if len(valid_endpoints) < 2:\n            return repaired_skeleton\n\n        # 3. 執行配對\n        matched_pairs = match_endpoints_2d(valid_endpoints, valid_vectors, label_slice, max_dist, max_angle_deg)\n\n        # 4. Bresenham 畫線補洞\n        for match in matched_pairs:\n            pA, pB = match['pA'], match['pB']\n            vA, vB = match['vA'], match['vB']\n            # skimage.draw.line 會回傳線上所有點的 (y, x) 座標\n            rr, cc = line(pA[0], pA[1], pB[0], pB[1])\n            repaired_skeleton[rr, cc] = 1\n\n        # 4. 貝茲曲線畫線補洞\n        # for match in matched_pairs:\n        #     pA, pB = match['pA'], match['pB']\n        #     vA, vB = match['vA'], match['vB']\n        #\n        #     # 取得曲線的 y, x 座標\n        #     rr, cc = generate_bezier_curve_2d(pA, pB, vA, vB, weight=0.4)\n        #\n        #     # 安全機制：確保座標沒有超出圖片邊界\n        #     valid = (rr >= 0) & (rr < skeleton.shape[0]) & \\\n        #             (cc >= 0) & (cc < skeleton.shape[1])\n        #\n        #     repaired_skeleton[rr[valid], cc[valid]] = 1\n\n        return repaired_skeleton\n\n    def connect_skeleton_holes_dijkstra(skeleton, slice_2d, label_slice, max_dist=30, max_angle_deg=45, trace_steps=5):\n        \"\"\"\n        主函數：利用 Dijkstra 與成本圖連接骨架\n        \"\"\"\n        repaired_skeleton = skeleton.copy()\n\n        # 1. 抓取端點\n        raw_endpoints = get_2d_endpoints(skeleton)\n        if len(raw_endpoints) < 2:\n            return repaired_skeleton\n\n        # 2. 【修改處】：直接先拿掉太短的骨架，不參與配對\n        valid_endpoints = []\n        valid_vectors = []\n        for ep in raw_endpoints:\n            # 這裡的 trace_steps 可以依需求設為 3 或 5\n            vec = trace_back_and_get_vector_2d(skeleton, ep, steps=trace_steps)\n\n            # 如果 vec 不是 None，代表這條骨架夠長，才允許加入配對池\n            if vec is not None:\n                valid_endpoints.append(ep)\n                valid_vectors.append(vec)\n\n        if len(valid_endpoints) < 2:\n            return repaired_skeleton\n\n        # 3. 執行配對\n        matched_pairs = match_endpoints_2d(valid_endpoints, valid_vectors, label_slice, max_dist, max_angle_deg)\n\n        if not matched_pairs:\n            return repaired_skeleton\n\n        # 3. 建立成本圖 (Cost Map) 的全新設計\n        # 初始化：預設所有地方的基礎成本為 1.0\n        cost_map = np.ones_like(slice_2d, dtype=np.float32)\n\n        # 定義區域\n        is_background = (slice_2d == 0)\n        is_current_mask = (slice_2d == label_slice)\n        is_other_mask = (~is_background) & (~is_current_mask)\n\n        # --- A. 建立避開「其他 Mask」的護城河 ---\n        # 若畫面上真的有其他的 Mask，才需要建立護城河\n        if np.any(is_other_mask):\n            # 計算每個像素到「其他 Mask」的距離\n            # edt 找 0，所以要把 is_other_mask 反轉 (其他 mask 處為 0，其餘為 1)\n            dist_to_others = distance_transform_edt(~is_other_mask)\n            max_warning_penalty = 20.0\n\n            # 距離越近，懲罰越高；大於 safe_distance 則懲罰為 0\n            warning_penalty = np.clip(max_warning_penalty * (1 - dist_to_others / safe_distance), 0,\n                                      max_warning_penalty)\n            cost_map += warning_penalty\n\n        # --- B. 設定各地形基礎成本 ---\n\n        # 1. 自己的 Mask 內部：保持最低成本，不加任何懲罰\n        # (cost_map 在這裡維持為 1.0 + 可能的護城河懲罰)\n\n        # 2. 黑色背景：給予基礎懲罰，強迫它在背景中只走最短路徑\n        bg_penalty = 1\n        cost_map[is_background] += bg_penalty\n\n        # 3. 其他 Mask：設定為絕對高牆，嚴禁穿越\n        cost_map[is_other_mask] = 1e6\n\n        # 4. 使用 Dijkstra 畫線補洞\n        for match in matched_pairs:\n            pA, pB = tuple(match['pA']), tuple(match['pB'])\n\n            try:\n                # 尋找兩點之間最小成本的路徑\n                path, cost = route_through_array(cost_map, pA, pB, fully_connected=True)\n\n                if cost >= max_dist * 2:\n                    # 代表成本過高，放棄修補\n                    continue\n\n                path_y = [p[0] for p in path]\n                path_x = [p[1] for p in path]\n\n                repaired_skeleton[path_y, path_x] = 1\n            except ValueError:\n                continue\n\n        return repaired_skeleton\n\n    if border > 0:\n        working_vol = vol_bool[border:-border, border:-border, border:-border].copy()\n    else:\n        working_vol = vol_bool.copy()\n\n    working_skel = np.zeros_like(working_vol)\n    labeled_vol, num_features = ndi.label(working_vol, structure=generate_binary_structure(rank=3, connectivity=1))\n\n    for _ in range(pass_iters):\n        for a in range(working_vol.shape[0]):\n            if axis == 0:\n                slice_2d = working_vol[a, :, :]\n                label_slice_2d = labeled_vol[a, :, :]\n            elif axis == 1:\n                slice_2d = working_vol[:, a, :]\n                label_slice_2d = labeled_vol[:, a, :]\n            else:\n                slice_2d = working_vol[:, :, a]\n                label_slice_2d = labeled_vol[:, :, a]\n\n            if not np.any(slice_2d): continue\n\n            # skel = skeletonize(slice_2d, method='lee')\n            skel = medial_axis(slice_2d)\n            # 傳入原來的 slice_2d 來計算成本圖\n            skel = connect_skeleton_holes_dijkstra(skeleton=skel, slice_2d=slice_2d, label_slice=label_slice_2d,\n                                                   max_dist=max_dist, max_angle_deg=max_angle_deg)\n            if np.any(skel):\n                # 距離變換擴張\n                reconstructed = distance_transform_edt(~skel) <= radius\n                if axis == 0:\n                    working_vol[a, :, :] = reconstructed\n                    working_skel[a, :, :] = skel\n                elif axis == 1:\n                    working_vol[:, a, :] = reconstructed\n                    working_skel[:, a, :] = skel\n                else:\n                    working_vol[:, :, a] = reconstructed\n                    working_skel[:, :, a] = skel\n\n    if border > 0:\n        # 建立與輸入相同尺寸的全 0 (False) 陣列\n        final_vol = np.zeros_like(vol_bool)\n        final_skel = np.zeros_like(vol_bool)\n\n        # 把處理好的區域填入中心\n        final_vol[border:-border, border:-border, border:-border] = working_vol\n        final_skel[border:-border, border:-border, border:-border] = working_skel\n\n        # 如果你希望邊緣區域保留「原始的 Mask」，可以使用 OR 運算：\n        final_vol = final_vol | vol_bool\n        # 若你只想保留處理過的中心區域，邊緣維持空白，就直接 return 下面這行\n        return final_vol, final_skel\n\n    return working_vol, working_skel\n\n\ndef normalize_segments_3d(vol_bool: np.ndarray, radius: float = 1.5, border: int = 8, axis: int = 0, sigma: float = 1.0,\n                          repair: bool = False, max_dist: int = 30, max_angle_deg: int = 30,\n                          safe_distance: int = 4, repair_dijkstra: bool = False, iterations: int = 1) -> np.ndarray:\n    \"\"\"\n    將 3D 體積切成 8 塊後，對每塊分別進行線段正規化。\n    axis: 0 為 Z 軸, 1 為 Y 軸, 2 為 X 軸\n    \"\"\"\n\n    normalized_vol = vol_bool.copy()\n    skel_vol = np.zeros_like(vol_bool)\n    labeled_vol, num_features = ndi.label(normalized_vol, structure=generate_binary_structure(rank=3, connectivity=1))\n    slices_list = _octant_slices(normalized_vol.shape, splits=(2, 2, 2))\n\n    # slices_list = [(slice(None), slice(None), slice(None))]\n\n    def generate_bezier_curve_2d(pA, pB, vA, vB, weight=0.4):\n        \"\"\"\n        生成 2D 三次方貝茲曲線的像素座標。\n        pA, pB: 起點與終點 (y, x)\n        vA, vB: 兩端的生長方向向量 (需為朝向對方的向外向量)\n        weight: 控制點延伸的力度，建議設為起終點距離的 0.3 ~ 0.5\n        \"\"\"\n        p0 = np.array(pA, dtype=float)\n        p3 = np.array(pB, dtype=float)\n\n        # 計算兩點距離，決定控制點要推多遠 (力度)\n        dist = np.linalg.norm(p3 - p0)\n        k = dist * weight\n\n        # 計算控制點 P1, P2\n        p1 = p0 + vA * k\n        p2 = p3 + vB * k\n\n        # 根據距離決定取樣點數量，確保畫出來的線是連續的 (多取一點保證不斷線)\n        num_points = int(dist * 2.5)\n        if num_points < 5:\n            num_points = 5\n\n        t = np.linspace(0, 1, num_points)[:, np.newaxis]\n\n        # 三次方貝茲曲線公式\n        curve = ((1 - t) ** 3 * p0 +\n                 3 * (1 - t) ** 2 * t * p1 +\n                 3 * (1 - t) * t ** 2 * p2 +\n                 t ** 3 * p3)\n\n        # 四捨五入轉為整數像素座標\n        curve_pixels = np.round(curve).astype(int)\n\n        # 去除重複的像素點，減少後續寫入陣列的負擔\n        curve_pixels = np.unique(curve_pixels, axis=0)\n\n        return curve_pixels[:, 0], curve_pixels[:, 1]\n\n    def get_2d_endpoints(skeleton):\n        # 建立一個中間為 0，周圍為 1 的 kernel\n        kernel = np.array([[1, 1, 1],\n                           [1, 0, 1],\n                           [1, 1, 1]])\n        neighbor_count = convolve(skeleton.astype(int), kernel, mode='constant', cval=0)\n\n        # 骨架點且鄰居只有 1 個\n        endpoints = np.argwhere((skeleton == 1) & (neighbor_count == 1))\n        return endpoints\n\n    def trace_back_and_get_vector_2d(skeleton, endpoint, steps=3):\n        \"\"\"\n        從端點沿著骨架往回走指定的步數，計算生長方向向量。\n        \"\"\"\n        current_pt = endpoint\n        visited = {tuple(endpoint)}\n\n        for _ in range(steps):\n            y, x = current_pt\n            y_min, y_max = max(0, y - 1), min(skeleton.shape[0], y + 2)\n            x_min, x_max = max(0, x - 1), min(skeleton.shape[1], x + 2)\n\n            neighborhood = skeleton[y_min:y_max, x_min:x_max]\n            neighbors = np.argwhere(neighborhood == 1)\n\n            # 轉換為全局座標\n            neighbors += np.array([y_min, x_min])\n\n            next_pt = None\n            for n in neighbors:\n                if tuple(n) not in visited:\n                    next_pt = n\n                    break\n\n            if next_pt is None:\n                break  # 骨架太短，提早走到盡頭\n\n            visited.add(tuple(next_pt))\n            current_pt = next_pt\n\n        # 向量方向：從回溯點指向端點 (即未來的生長趨勢)\n        vec = endpoint - current_pt\n        norm = np.linalg.norm(vec)\n        if norm == 0:\n            return np.array([0.0, 0.0])\n        return vec / norm\n\n    def match_endpoints_2d(endpoints, vectors, label_slice, max_dist=30, max_angle_deg=45):\n        \"\"\"\n        根據距離和方向夾角配對端點\n        \"\"\"\n        matched_pairs = []\n        used_indices = set()\n        max_angle_rad = math.radians(max_angle_deg)\n\n        for i in range(len(endpoints)):\n            if i in used_indices: continue\n\n            best_match = -1\n            min_dist = max_dist\n\n            for j in range(len(endpoints)):\n                if i == j or j in used_indices: continue\n\n                p_A, v_A = endpoints[i], vectors[i]\n                p_B, v_B = endpoints[j], vectors[j]\n\n                label_A = label_slice[p_A[0], p_A[1]]\n                label_B = label_slice[p_B[0], p_B[1]]\n\n                if label_A != label_B:\n                    continue\n\n                # 1. 檢查距離\n                dist = np.linalg.norm(p_A - p_B)\n                if dist > max_dist: continue\n\n                # 2. 檢查方向性\n                vec_AB = (p_B - p_A) / dist\n                vec_BA = -vec_AB\n\n                # 確保向量長度不為 0 才計算夾角\n                if np.linalg.norm(v_A) == 0 or np.linalg.norm(v_B) == 0: continue\n\n                # v_A 應該要指向 B，v_B 應該要指向 A\n                # cos_A = np.clip(np.dot(v_A, vec_AB), -1.0, 1.0)\n                # cos_B = np.clip(np.dot(v_B, vec_BA), -1.0, 1.0)\n                #\n                # if math.acos(cos_A) < max_angle_rad and math.acos(cos_B) < max_angle_rad:\n                #     if dist < min_dist:  # 尋找距離最近的合格者\n                #         min_dist = dist\n                #         best_match = j\n\n                dot_A_forward = np.dot(v_A, vec_AB)\n                dot_B_forward = np.dot(v_B, vec_BA)\n\n                # 2. 檢查兩個端點本身的生長向量是否「相向」(例如夾角大於 135 度)\n                # v_B 應該要和 v_A 大致反向，所以 v_A 和 -v_B 應該要大致同向\n                cos_dirs = np.clip(np.dot(v_A, -v_B), -1.0, 1.0)\n                dirs_aligned = math.acos(cos_dirs) < max_angle_rad  # 這裡的 max_angle_rad 可以設為 45 度 (即允許 45 度的方向誤差)\n\n                if dot_A_forward > 0 and dot_B_forward > 0 and dirs_aligned:\n                    if dist < min_dist:\n                        min_dist = dist\n                        best_match = j\n\n            if best_match != -1:\n                matched_pairs.append({\n                    'pA': endpoints[i], 'pB': endpoints[best_match],\n                    'vA': vectors[i], 'vB': vectors[best_match]\n                })\n                used_indices.add(i)\n                used_indices.add(best_match)\n\n        return matched_pairs\n\n    def connect_skeleton_holes_2d(skeleton, label_slice, max_dist=30, max_angle_deg=45, trace_steps=5):\n        \"\"\"\n        主函數：輸入 2D 骨架矩陣，輸出修補好的骨架矩陣\n        \"\"\"\n        # 複製一份準備畫線用\n        repaired_skeleton = skeleton.copy()\n\n        # 1. 抓取端點\n        raw_endpoints = get_2d_endpoints(skeleton)\n        if len(raw_endpoints) < 2:\n            return repaired_skeleton\n\n        # 2. 【修改處】：直接先拿掉太短的骨架，不參與配對\n        valid_endpoints = []\n        valid_vectors = []\n        for ep in raw_endpoints:\n            # 這裡的 trace_steps 可以依需求設為 3 或 5\n            vec = trace_back_and_get_vector_2d(skeleton, ep, steps=trace_steps)\n\n            # 如果 vec 不是 None，代表這條骨架夠長，才允許加入配對池\n            if vec is not None:\n                valid_endpoints.append(ep)\n                valid_vectors.append(vec)\n\n        if len(valid_endpoints) < 2:\n            return repaired_skeleton\n\n        # 3. 執行配對\n        matched_pairs = match_endpoints_2d(valid_endpoints, valid_vectors, label_slice, max_dist, max_angle_deg)\n\n        # 4. Bresenham 畫線補洞\n        for match in matched_pairs:\n            pA, pB = match['pA'], match['pB']\n            vA, vB = match['vA'], match['vB']\n            # skimage.draw.line 會回傳線上所有點的 (y, x) 座標\n            rr, cc = line(pA[0], pA[1], pB[0], pB[1])\n            repaired_skeleton[rr, cc] = 1\n\n        # 4. 貝茲曲線畫線補洞\n        # for match in matched_pairs:\n        #     pA, pB = match['pA'], match['pB']\n        #     vA, vB = match['vA'], match['vB']\n        #\n        #     # 取得曲線的 y, x 座標\n        #     rr, cc = generate_bezier_curve_2d(pA, pB, vA, vB, weight=0.4)\n        #\n        #     # 安全機制：確保座標沒有超出圖片邊界\n        #     valid = (rr >= 0) & (rr < skeleton.shape[0]) & \\\n        #             (cc >= 0) & (cc < skeleton.shape[1])\n        #\n        #     repaired_skeleton[rr[valid], cc[valid]] = 1\n\n        return repaired_skeleton\n\n    def connect_skeleton_holes_dijkstra(skeleton, slice_2d, label_slice, max_dist=30, max_angle_deg=45, trace_steps=5):\n        \"\"\"\n        主函數：利用 Dijkstra 與成本圖連接骨架\n        \"\"\"\n        repaired_skeleton = skeleton.copy()\n\n        # 1. 抓取端點\n        raw_endpoints = get_2d_endpoints(skeleton)\n        if len(raw_endpoints) < 2:\n            return repaired_skeleton\n\n        # 2. 直接先拿掉太短的骨架，不參與配對\n        valid_endpoints = []\n        valid_vectors = []\n        for ep in raw_endpoints:\n            vec = trace_back_and_get_vector_2d(skeleton, ep, steps=trace_steps)\n            if vec is not None:\n                valid_endpoints.append(ep)\n                valid_vectors.append(vec)\n\n        if len(valid_endpoints) < 2:\n            return repaired_skeleton\n\n        # 3. 執行配對\n        matched_pairs = match_endpoints_2d(valid_endpoints, valid_vectors, label_slice, max_dist, max_angle_deg)\n        if not matched_pairs:\n            return repaired_skeleton\n\n        # ==========================================\n        # 【新增處】：利用連通集區分真正的目標區域與無關孤島\n        # ==========================================\n        # 取出當前 label 的所有區域\n        current_mask_full = (slice_2d == label_slice)\n\n        # 標記所有連通的孤島\n        labeled_mask, num_features = ndi.label(current_mask_full)\n\n        # 找出包含 valid_endpoints 的所有孤島 ID\n        target_island_ids = set()\n        for ep in valid_endpoints:\n            island_id = labeled_mask[ep[0], ep[1]]\n            if island_id != 0:\n                target_island_ids.add(island_id)\n\n        # 定義區域\n        # is_target_mask: 只有包含端點的孤島才享有最低成本\n        is_target_mask = np.isin(labeled_mask, list(target_island_ids))\n\n        # is_disconnected_same_label: 同 label 但不含端點的孤島 (視同背景處理，避免被借道)\n        is_disconnected_same_label = current_mask_full & (~is_target_mask)\n\n        is_background = (slice_2d == 0) | is_disconnected_same_label\n        is_other_mask = (~is_background) & (~is_target_mask)\n\n        # ==========================================\n\n        # 4. 建立成本圖 (Cost Map)\n        cost_map = np.ones_like(slice_2d, dtype=np.float32)\n\n        # --- A. 建立避開「其他 Mask」的護城河 ---\n        if np.any(is_other_mask):\n            dist_to_others = distance_transform_edt(~is_other_mask)\n            max_warning_penalty = 15.0\n\n            warning_penalty = np.clip(max_warning_penalty * (1 - dist_to_others / safe_distance), 0,\n                                      max_warning_penalty)\n            cost_map += warning_penalty\n\n        # --- B. 設定各地形基礎成本 ---\n\n        # 1. 自己的目標 Mask 內部：保持最低成本 1.0\n        # (cost_map[is_target_mask] 維持原樣)\n\n        # 2. 黑色背景與「無關的同 label 孤島」：給予基礎懲罰\n        bg_penalty = 2.0  # 建議稍微提高背景懲罰，強迫走最短直線\n        cost_map[is_background] = bg_penalty\n\n        # 3. 其他 Mask：絕對高牆\n        cost_map[is_other_mask] = 1e6\n\n        # 5. 使用 Dijkstra 畫線補洞\n        for match in matched_pairs:\n            pA, pB = tuple(match['pA']), tuple(match['pB'])\n\n            try:\n                path, cost = route_through_array(cost_map, pA, pB, fully_connected=True)\n\n                if cost >= max_dist * bg_penalty * 2:  # 注意這裡的閾值可能需要配合 bg_penalty 調整\n                    continue\n\n                path_y = [p[0] for p in path]\n                path_x = [p[1] for p in path]\n\n                repaired_skeleton[path_y, path_x] = 1\n            except ValueError:\n                continue\n\n        return repaired_skeleton\n\n    def _process_block(block: np.ndarray, labeled_block: np.ndarray, radius: float, border: int,\n                       axis: int, sigma: float = 0.0) -> tuple[np.ndarray, np.ndarray]:\n\n        block_out = block.copy()\n        skel_out = block.copy()\n        shape = block.shape\n\n        # 1. 檢查非處理軸的其餘兩個維度是否夠大\n        other_dims = [shape[i] for i in range(3) if i != axis]\n        if border > 0 and any(d <= 2 * border for d in other_dims):\n            return block_out, skel_out  # 修正：確保提早結束時也回傳兩個陣列\n\n        # 2. 定義統一的內部區域切片邊界\n        inner_slice = slice(border, -border if border > 0 else None)\n\n        # 3. 沿著指定的 axis 進行迭代\n        for i in range(shape[axis]):\n            # 動態生成 3D 讀取切片 (例如 axis=1, i=5 時，等同於 [:, 5, :])\n            read_idx = [slice(None)] * 3\n            read_idx[axis] = i\n            read_idx = tuple(read_idx)\n\n            slice_2d = block[read_idx]\n            if not np.any(slice_2d):\n                continue\n\n            label_slice_2d = labeled_block[read_idx]\n            skel = skeletonize(slice_2d, method='lee')\n            if repair:\n                for it in range(iterations):\n                    if repair_dijkstra:\n                        skel = connect_skeleton_holes_dijkstra(skeleton=skel, slice_2d=slice_2d,\n                                                               label_slice=label_slice_2d,\n                                                               max_dist=max_dist, max_angle_deg=max_angle_deg)\n                    else:\n                        skel = connect_skeleton_holes_2d(skeleton=skel, label_slice=label_slice_2d, max_dist=max_dist,\n                                                         max_angle_deg=max_angle_deg)\n\n            if np.any(skel):\n                # 距離變換擴張\n                reconstructed = distance_transform_edt(~skel) <= radius\n\n                # 加入 2D 高斯模糊\n                if sigma > 0:\n                    reconstructed = gaussian_filter(reconstructed.astype(float), sigma=sigma)\n                    reconstructed = reconstructed > 0.3\n\n                # 動態生成 3D 寫回切片 (保留 border，並指定當前層 i)\n                write_idx = [inner_slice] * 3\n                write_idx[axis] = i\n                write_idx = tuple(write_idx)\n\n                # reconstructed 是 2D，我們只需要擷取它內部的部分來填入 3D 結構\n                recon_idx = tuple([inner_slice] * 2)\n\n                block_out[write_idx] = reconstructed[recon_idx]\n                skel_out[write_idx] = skel[recon_idx]\n\n        return block_out, skel_out\n\n    # 疊代處理 8 個子塊\n    for sl in slices_list:\n        sub_vol = vol_bool[sl]\n        labeled_sub_vol = labeled_vol[sl]\n        processed_sub_vol, skel_sub_vol = _process_block(sub_vol, labeled_sub_vol, radius, border, axis)\n        normalized_vol[sl] = processed_sub_vol\n        skel_vol[sl] = skel_sub_vol\n\n    return normalized_vol, skel_vol\n\n\ndef apply_closing(vol, iterations: int = 1, conn: int = 3):\n    \"\"\"3D方向閉運算\"\"\"\n    m = vol.astype(bool)\n    struct = ndi.generate_binary_structure(3, conn)\n    closed = ndi.binary_closing(m, structure=struct, iterations=iterations)\n    return closed.astype(np.uint8)\n\n\ndef apply_z_axis_closing(vol, iterations=1):\n    \"\"\"Z 軸方向閉運算\"\"\"\n    processed_vol = np.copy(vol)\n    for z in range(vol.shape[1]):\n        slice_2d = vol[z, :, :]\n        if np.any(slice_2d):\n            closed_slice = binary_closing(slice_2d, iterations=iterations)\n            processed_vol[z, :, :] = closed_slice\n    return processed_vol\n\n\ndef apply_y_axis_closing(vol, iterations=1):\n    \"\"\"Y 軸方向閉運算\"\"\"\n    processed_vol = np.copy(vol)\n    for y in range(vol.shape[1]):\n        slice_2d = vol[:, y, :]\n        if np.any(slice_2d):\n            closed_slice = binary_closing(slice_2d, iterations=iterations)\n            processed_vol[:, y, :] = closed_slice\n    return processed_vol\n\n\ndef apply_z_axis_fill_holes(vol):\n    \"\"\"Z 軸方向孔洞填充\"\"\"\n    processed_vol = np.copy(vol)\n    for z in range(vol.shape[0]):\n        slice_2d = vol[z, :, :]\n        if np.any(slice_2d):\n            filled_slice = binary_fill_holes(slice_2d)\n            processed_vol[z, :, :] = filled_slice\n    return processed_vol\n\n\ndef apply_y_axis_fill_holes(vol):\n    \"\"\"Y 軸方向孔洞填充\"\"\"\n    processed_vol = np.copy(vol)\n    for y in range(vol.shape[1]):\n        slice_2d = vol[:, y, :]\n        if np.any(slice_2d):\n            filled_slice = binary_fill_holes(slice_2d)\n            processed_vol[:, y, :] = filled_slice\n    return processed_vol\n\n\ndef compute_separation_constraint(mask):\n    \"\"\"\n    計算連通域間的不可侵犯領域與邊界。\n\n    規則：\n    1. 定義物件：使用 6 連通 (Face-connected) 區分物件。\n    2. 劃分領土：使用 Voronoi (EDT) 擴張，填滿背景。\n    3. 定義邊界：若體素的 26 連通鄰域 (3x3x3) 內包含 \">=2 個不同的連通域 ID\"，則視為邊界。\n    \"\"\"\n\n    # -----------------------------------------------------------\n    # 1. 計算 3D 連通域 (Labeling) - 嚴格區分對角線接觸\n    # -----------------------------------------------------------\n    mask = mask.copy()\n    mask = remove_small_objects(mask > 0, 5000)\n    # 使用 6 連通 (connectivity=1)，確保僅在對角線接觸的物件被視為不同個體\n    structure_6conn = ndi.generate_binary_structure(rank=3, connectivity=1)\n    labeled_array, num_features = ndi.label(mask, structure=structure_6conn)\n\n    # 如果全場只有 0 或 1 個物件，不存在\"交界\"，直接回傳全 True\n    if num_features <= 1:\n        return np.ones_like(mask, dtype=bool)\n\n    # 2. 擴張勢力範圍 (Voronoi / EDT)\n    # 我們需要一個\"填滿\"的空間圖，每個點都知道自己最近的物件是誰 (Territory ID)\n    # 對\"背景 (0)\"做距離變換，找到最近的前景索引\n    # indices shape: (3, D, H, W)\n    _, indices = ndi.distance_transform_edt(labeled_array == 0, return_indices=True, return_distances=True)\n\n    # 映射回 Label ID，得到全空間的領土圖\n    # 這裡 territory 不再有 0 (背景)，全都是 1~N 的 ID\n    territory = labeled_array[indices[0], indices[1], indices[2]]\n\n    # 3. 標記交界處 (26-Connectivity Multi-Label Detection)\n    # 規則：如果一個體素的 26 連通鄰域內，包含 2 個以上的不同 ID，它就是邊界。\n    # 數學實作：在 territory map 上，若 Max(鄰域) != Min(鄰域)，則必有至少兩個不同 ID。\n\n    footprint_26conn = np.ones((3, 3, 3), dtype=int)  # 26 連通核心\n\n    # 找出鄰域內的最大 ID 與最小 ID\n    max_labels = ndi.maximum_filter(territory, footprint=footprint_26conn)\n    min_labels = ndi.minimum_filter(territory, footprint=footprint_26conn)\n\n    # 邊界判定：只要最大值不等於最小值，代表鄰域內混雜了不同的領土\n    # 這會偵測到：\n    # 1. 兩個物件擴張後的接觸面 (Voronoi 邊界)\n    # 2. 兩個原始物件物理上靠得很近或接觸的地方 (Original Contact)\n    boundary_mask = (max_labels != min_labels)\n\n    # 4. 取得最終 Valid Mask\n    # 只要不是邊界，就是合法填充區\n    valid_mask = ~boundary_mask\n\n    return valid_mask\n\n\ndef recover_weak_signals(\n        vol: np.ndarray,\n        prob: np.ndarray,\n        fg_threshold: float = 0.3,\n        fg_bg_diff_threshold: float = 0.3\n) -> np.ndarray:\n    \"\"\"\n    基於機率圖的弱訊號恢復 (Weak Signal Recovery)。\n    將「模型預測為前景(>fg_thresh) 但信心不足((fg-bg) <= diff_thresh)」的區域撈回。\n\n    Args:\n        vol: 原始預測的 3D Mask (D, H, W)，通常是 argmax 後的結果 (0=BG, >0=FG)。\n        prob: 機率圖 (Probability Map)。\n              支援 (C, D, H, W) 或 (D, H, W, C) 格式。\n              假設 C=2 (Index 0=BG, Index 1=FG)。\n        fg_threshold: 前景機率門檻值。\n        fg_bg_diff_threshold: (前景 - 背景) 的信心差門檻值。\n\n    Returns:\n        np.ndarray: 疊加恢復訊號後的 3D Mask (uint8, 0 or 1)。\n    \"\"\"\n    # 檢查是否傳入機率圖\n    if prob is None:\n        return vol.astype(np.uint8)\n\n    # 檢查維度是否匹配 (prob 應該比 vol 多一個 channel 維度)\n    if prob.ndim != vol.ndim + 1:\n        print(f\"[Warning] Prob shape {prob.shape} does not match Vol shape {vol.shape} (+1 dim). Skipping recovery.\")\n        return vol.astype(np.uint8)\n\n    print(f\"正在執行弱訊號恢復 (FG > {fg_threshold}, Diff <= {fg_bg_diff_threshold})...\")\n\n    # 自動判斷 Channel 維度 (假設 2 類別: 0=BG, 1=FG)\n    # 判斷 Channel 是在第一個維度 (PyTorch style) 還是最後一個維度 (TensorFlow style)\n    # 這裡假設 D, H, W 通常都大於 Channel 數 (2)\n    if prob.shape[0] < prob.shape[-1]:\n        # Shape: (C, D, H, W)\n        bg_prob = prob[0, ...]\n        fg_prob = prob[1, ...]\n    else:\n        # Shape: (D, H, W, C)\n        bg_prob = prob[..., 0]\n        fg_prob = prob[..., 1]\n\n    # 計算恢復條件\n    # 1. 前景機率需大於基本門檻\n    cond_fg_strong = (fg_prob > fg_threshold)\n    # 2. 前景與背景的差距小於設定值 (代表模型猶豫)\n    cond_diff_small = ((fg_prob - bg_prob) > fg_bg_diff_threshold)\n\n    recover_mask = cond_fg_strong & cond_diff_small\n\n    # 將恢復的像素疊加到原有的預測上 (聯集 OR 操作)\n    # vol > 0 確保將輸入視為二值\n    final_vol = np.logical_or(vol > 0, recover_mask)\n\n    return final_vol.astype(np.uint8)\n\n\ndef _axis_cuts(n: int, parts: int) -> List[Tuple[int, int]]:\n    \"\"\"將維度 n 切分為 parts 份，處理餘數確保覆蓋全域。\"\"\"\n    parts = max(1, min(parts, n if n > 0 else 1))\n    base = n // parts\n    extra = n % parts\n    cuts = []\n    start = 0\n    for i in range(parts):\n        size = base + (1 if i < extra else 0)\n        end = start + size\n        cuts.append((start, end))\n        start = end\n    return cuts\n\n\ndef _octant_slices(shape: Tuple[int, int, int], splits: Tuple[int, int, int]) -> List[Tuple[slice, slice, slice]]:\n    \"\"\"生成 3D 切片物件清單。\"\"\"\n    zcuts = _axis_cuts(shape[0], splits[0])\n    ycuts = _axis_cuts(shape[1], splits[1])\n    xcuts = _axis_cuts(shape[2], splits[2])\n    out = []\n    for z0, z1 in zcuts:\n        for y0, y1 in ycuts:\n            for x0, x1 in xcuts:\n                if (z1 - z0) > 0 and (y1 - y0) > 0 and (x1 - x0) > 0:\n                    out.append((slice(z0, z1), slice(y0, y1), slice(x0, x1)))\n    return out\n\n\ndef remove_small_8slice(vol, min_size=100, connectivity=1):\n    \"\"\"使用 8-slice 邏輯平行化/分段移除小物件。\"\"\"\n    vol_bool = vol > 0\n    output = np.zeros_like(vol_bool)\n\n    # 取得 8 個區塊的 slice (2x2x2)\n    slices = _octant_slices(vol_bool.shape, (2, 2, 2))\n\n    for slc in slices:\n        block = vol_bool[slc]\n        if np.any(block):\n            # 處理區塊並塞回對應位置\n            output[slc] = remove_small_objects(\n                block,\n                min_size=min_size,\n                connectivity=connectivity\n            )\n\n    return output.astype(np.uint8)\n\n\ndef robust_mask_sandwich_v3(mask, max_gap=2, axis_weights=[1, 1, 1], min_neighbors=2, kernel_size=5):\n    \"\"\"\n    V3 優化強健版（支援擴大 Kernel）：\n    1. 預計算支撐圖 (Support Map)，大幅提升效能。\n    2. 全軸向支持：所有軸向均可進行 2D 支撐檢查。\n    3. 記憶體優化：最小化陣列複製次數。\n    4. 動態 Kernel：可透過 kernel_size 參數（需為奇數）擴大鄰居檢查範圍。\n    \"\"\"\n    if kernel_size % 2 == 0:\n        raise ValueError(\"kernel_size 必須為奇數（例如 3, 5, 7）\")\n\n    # 統一轉換為 bool 進行運算，減少記憶體壓力\n    refined = mask.astype(bool, copy=True)\n    ndim = 3\n    half_k = kernel_size // 2\n\n    # --- 預計算：每一層的支撐圖 (只算一次) ---\n    # 這樣在後續遍歷不同的 gap 時，不需要重複捲積\n    support_masks = []\n    for axis in range(ndim):\n        if axis_weights[axis] == 0:\n            support_masks.append(None)\n            continue\n\n        # 建立該軸向的 2D 鄰居統計 Kernel\n        # 若 kernel_size=5，則 k_shape 預設為 [5, 5, 5]\n        k_shape = [kernel_size] * ndim\n        k_shape[axis] = 1  # 讓 Kernel 在該軸向上是扁平的 (例如 [1, 5, 5])\n        kernel = np.ones(k_shape, dtype=np.uint8)\n\n        # 移除中心點 (不計算自己)\n        # 動態定位中心，例如 kernel_size=5，half_k=2，center=[2, 2, 2]\n        center = [half_k] * ndim\n        center[axis] = 0\n        kernel[tuple(center)] = 0\n\n        # 一次性計算整體的鄰居數量\n        neighbor_count = convolve(refined.astype(np.uint8), kernel, mode='constant', cval=0)\n\n        # 只有「自己是1」且「鄰居夠多」的才算有效錨點\n        support_masks.append((refined) & (neighbor_count >= min_neighbors))\n\n    # --- 執行填補 ---\n    for axis in range(ndim):\n        if axis_weights[axis] == 0:\n            continue\n\n        valid_map = support_masks[axis]\n        shape_at_axis = refined.shape[axis]\n\n        for gap in range(1, max_gap + 1):\n            stride = gap + 1\n\n            # 取得兩端的「強健錨點」\n            s_prev = [slice(None)] * ndim\n            s_next = [slice(None)] * ndim\n            s_prev[axis] = slice(0, -stride)\n            s_next[axis] = slice(stride, None)\n\n            # 只有兩端都是「結實像素」時，才建立橋樑\n            bridge = valid_map[tuple(s_prev)] & valid_map[tuple(s_next)]\n\n            if not np.any(bridge):\n                continue\n\n            # 填補中間空隙\n            for i in range(1, stride):\n                s_fill = [slice(None)] * ndim\n                s_fill[axis] = slice(i, i + (shape_at_axis - stride))\n                refined[tuple(s_fill)] |= bridge\n\n    return refined\n\n\ndef robust_sandwich_fill(mask, max_gap=2, fix_planes=('XZ', 'YZ'), kernel_size=3, min_neighbors=2):\n    \"\"\"\n    支援大空隙 (max_gap) 的對角線強健版填補（動態更新錨點模式）。\n    每次填補後會即時更新可用的「強健錨點」，允許產生連鎖填補效應。\n    \"\"\"\n    if kernel_size % 2 == 0:\n        raise ValueError(\"kernel_size 必須為奇數（例如 3, 5, 7）\")\n\n    # 統一轉換為 bool 進行運算，減少記憶體壓力\n    img = mask.astype(bool, copy=True)\n    refined = img.copy()\n    ndim = 3\n    shape = refined.shape\n\n    # 準備 3D 鄰居統計 Kernel (若需要判斷鄰居才會用到)\n    if min_neighbors > 0:\n        support_kernel = np.ones((kernel_size, kernel_size, kernel_size), dtype=np.uint8)\n        half_k = kernel_size // 2\n        support_kernel[half_k, half_k, half_k] = 0\n\n    # 定義要修復的平面對角線方向 (dz, dy, dx)\n    plane_map = {\n        'XZ': [(1, 0, 1), (1, 0, -1)],\n        'YZ': [(1, 1, 0), (1, -1, 0)],\n        'XY': [(0, 1, 1), (0, 1, -1)]\n    }\n\n    selected_dirs = []\n    for p in fix_planes:\n        if p in plane_map:\n            selected_dirs.extend(plane_map[p])\n\n    # 執行切片位移對角線填補\n    for dz, dy, dx in selected_dirs:\n        dir_vector = (dz, dy, dx)\n\n        for gap in range(1, max_gap + 1):\n\n            # --- 【核心修改】：動態更新錨點 ---\n            # 每次處理新方向或新間距前，根據「最新的 refined 狀態」重新計算錨點\n            if min_neighbors > 0:\n                # 重新執行卷積以獲取最新鄰居數量\n                neighbor_count = convolve(refined.astype(np.uint8), support_kernel, mode='constant', cval=0)\n                valid_anchors = refined & (neighbor_count >= min_neighbors)\n            else:\n                # 效能優化：如果 min_neighbors=0，最新的錨點就等於當前的 refined 狀態，免算卷積\n                valid_anchors = refined.copy()\n            # ----------------------------------\n\n            stride = gap + 1\n            stride_offset = tuple(d * stride for d in dir_vector)\n\n            s_base = []\n            s_target = []\n            lengths = []\n            starts = []\n\n            for a in range(ndim):\n                offset_a = stride_offset[a]\n                if offset_a > 0:\n                    s_base.append(slice(0, -offset_a))\n                    s_target.append(slice(offset_a, None))\n                    lengths.append(shape[a] - offset_a)\n                    starts.append(0)\n                elif offset_a < 0:\n                    s_base.append(slice(-offset_a, None))\n                    s_target.append(slice(0, offset_a))\n                    lengths.append(shape[a] - abs(offset_a))\n                    starts.append(-offset_a)\n                else:\n                    s_base.append(slice(None))\n                    s_target.append(slice(None))\n                    lengths.append(shape[a])\n                    starts.append(0)\n\n            # 找出兩端都是「最新強健錨點」的橋樑\n            bridge = valid_anchors[tuple(s_base)] & valid_anchors[tuple(s_target)]\n\n            if not np.any(bridge):\n                continue\n\n            # 將橋樑中間的「所有」空隙點補齊\n            for i in range(1, stride):\n                fill_offset = tuple(d * i for d in dir_vector)\n                s_fill = []\n\n                for a in range(ndim):\n                    start = starts[a] + fill_offset[a]\n                    end = start + lengths[a]\n                    s_fill.append(slice(start, end))\n\n                # 更新 refined，供下一個迴圈使用\n                refined[tuple(s_fill)] |= bridge\n\n    return refined\n\n\ndef robust_mask_refine(mask, max_gap=2, iterations=1, min_neighbors=2, kernel_size=5, keep_connectivity=True):\n    \"\"\"\n    包含連通性保護的綜合修復流程。\n\n    Args:\n        keep_connectivity (bool): 若為 True，當填補導致連通數量改變時，撤銷該次操作。\n    \"\"\"\n\n    def get_conn_count(mask, structure=None):\n        \"\"\"\n        計算 3D 遮罩中的連通域數量。\n        structure=None 預設使用 6-連通 (面相鄰)。\n        若要考慮對角線連通，可自定義 3x3x3 結構。\n        \"\"\"\n        _, count = ndi.label(mask, structure=structure)\n        return count\n\n    res = mask.copy()\n\n    # 預先定義連通性結構 (3D 26-connectivity 較為嚴謹)\n    conn_structure = np.ones((3, 3, 3), dtype=int)\n    structure_6_connected = np.array([\n        [[0, 0, 0], [0, 1, 0], [0, 0, 0]],\n        [[0, 1, 0], [1, 1, 1], [0, 1, 0]],\n        [[0, 0, 0], [0, 1, 0], [0, 0, 0]]\n    ])\n\n    for i in range(iterations):\n        # 紀錄本次迭代前的狀態與連通數\n        prev_res = res.copy()\n        if keep_connectivity:\n            prev_count = get_conn_count(prev_res, structure=structure_6_connected)\n\n        res = robust_mask_sandwich_v3(res, max_gap=max_gap, min_neighbors=min_neighbors, kernel_size=kernel_size)\n        res = robust_sandwich_fill(res, max_gap=2, min_neighbors=0, kernel_size=5,\n                                   fix_planes=('XZ', 'YZ', 'XY'))\n        res = robust_sandwich_fill(res, max_gap=2, min_neighbors=0, kernel_size=5,\n                                   fix_planes=('XZ', 'YZ', 'XY'))\n\n        # --- 連通性檢查 ---\n        if keep_connectivity:\n            current_count = get_conn_count(res, structure=structure_6_connected)\n\n            if current_count != prev_count:\n                # 發生合併(數量減少)或產生孤島(數量增加，雖然在此情境較少見)\n                # 撤銷此輪迭代\n                # print(f\"Iteration {i+1}: Connectivity changed ({prev_count} -> {current_count}). Rollback.\")\n                res = prev_res\n            else:\n                pass\n                # print(f\"Iteration {i+1}: Success. Connectivity maintained at {current_count}.\")\n\n    return res\n\n\ndef fill_2x2diag(\n        mask_3d: np.ndarray,\n        axis: int = 0,\n        gap: int = 0,\n) -> np.ndarray:\n    \"\"\"\n    掃描相鄰切片，修補 2x2 對角跳空 (diagonal gap) 問題。\n\n    Args:\n        mask_3d: 3D binary mask，shape = (z, y, x)\n        axis:    沿哪個軸做切片比較。0 = z 軸（預設），1 = y 軸，2 = x 軸\n        gap:     prev 與 curr 之間允許的中間層數。\n                 0 = 原始行為（相鄰兩層）。\n                 N = prev 與 curr 相距 N+1 層，中間 N 層需滿足\n                     連續相鄰對在該 2x2 的 AND 總和 > 1，\n                     確認有連續性後再判斷 prev/curr 的 3/3/4 對角條件，\n                     全部通過才將整段全設為 True。\n\n    Returns:\n        修補後的 3D binary mask（不修改原始輸入）\n    \"\"\"\n    if axis not in (0, 1, 2):\n        raise ValueError(f\"axis 必須是 0、1 或 2，收到 {axis}\")\n    if gap < 0:\n        raise ValueError(f\"gap 必須 >= 0，收到 {gap}\")\n\n    work = np.moveaxis(mask_3d, axis, 0).copy().astype(bool)\n    d, h, w = work.shape\n\n    step = gap + 1  # prev 與 curr 的 index 距離\n    corner_offsets = [(0, 0), (0, 1), (1, 0), (1, 1)]\n\n    for z in range(step, d):\n        prev_idx = z - step\n        curr_idx = z  # inclusive，共 step+1 層需要填補\n\n        prev = work[prev_idx].astype(np.uint8)\n        curr = work[curr_idx].astype(np.uint8)\n        orr = prev | curr\n\n        for dr, dc in corner_offsets:\n            r0, r1 = dr, h - 1 - dr  # 2x2 左上角的 row 範圍 [r0, r1)\n            c0, c1 = dc, w - 1 - dc\n\n            if r1 <= r0 or c1 <= c0:\n                continue\n\n            def box_sum(arr: np.ndarray) -> np.ndarray:\n                return (arr[r0: r1, c0: c1]\n                        + arr[r0: r1, c0 + 1: c1 + 1]\n                        + arr[r0 + 1: r1 + 1, c0: c1]\n                        + arr[r0 + 1: r1 + 1, c0 + 1: c1 + 1])\n\n            # ── 條件一：prev / curr 滿足 3/3/4 對角規則 ──────────────────\n            cond = (box_sum(prev) == 3) & (box_sum(curr) == 3) & (box_sum(orr) == 4)\n\n            if not np.any(cond):\n                continue\n\n            # ── 條件二：中間所有相鄰對的 AND 在該 2x2 內 > 1 ─────────────\n            # 遍歷 (prev_idx, prev_idx+1), ..., (curr_idx-1, curr_idx)\n            # gap=0 時此迴圈不執行\n            if gap > 0:\n                for i in range(prev_idx, curr_idx):\n                    a = work[i].astype(np.uint8)\n                    b = work[i + 1].astype(np.uint8)\n                    cond = cond & (box_sum(a & b) > 1)\n\n                    if not np.any(cond):\n                        break  # 提早結束，此 corner 已無命中點\n\n            # ── 填補：整段 [prev_idx, curr_idx] 的該 2x2 全設為 True ─────\n            hit_rows, hit_cols = np.where(cond)\n            if len(hit_rows) == 0:\n                continue\n\n            abs_rows = hit_rows + r0\n            abs_cols = hit_cols + c0\n\n            for rr, cc in zip(abs_rows, abs_cols):\n                work[prev_idx: curr_idx + 1, rr: rr + 2, cc: cc + 2] = True\n\n    return np.moveaxis(work, 0, axis)\n\n\ndef process_gaussian_octants_dual_threshold(\n        vol: np.ndarray,\n        sigma: tuple = (1.0, 1.0, 1.0),\n        inner_threshold: float = 0.5,\n        border_threshold: float = 0.1,\n        border: int = 8,\n        cval: float = 0.0,\n        mode: str = \"reflect\"\n) -> np.ndarray:\n    \"\"\"\n    將 3D 體積切塊處理：\n    1. 內部區域 (Inner): 使用 inner_threshold (預設 0.5)\n    2. 邊界區域 (Border): 使用 border_threshold (預設 0.1)\n    \"\"\"\n    output_vol = vol.copy()\n\n    # 取得 8 塊區域的切片範圍\n    slices_list = _octant_slices(vol.shape, splits=(2, 2, 2))\n\n    for sl in slices_list:\n        # 取出該區塊\n        sub_vol = vol[sl]\n        depth, height, width = sub_vol.shape\n\n        # 防呆：確保區塊大於兩倍 border\n        if height <= 2 * border or width <= 2 * border or depth <= 2 * border:\n            # 若區塊太小無法區分 border，則統一使用 border_threshold 處理或跳過\n            blurred_small = ndi.gaussian_filter(sub_vol.astype(float), sigma=sigma, mode=mode,\n                                                truncate=4.0)\n            output_vol[sl] = (blurred_small > border_threshold).astype(vol.dtype)\n            continue\n\n        # 核心：高斯模糊運算\n        blurred_sub = ndi.gaussian_filter(sub_vol.astype(float), sigma=sigma, mode=mode,\n                                          truncate=4.0, cval=cval)\n        # 建立該區塊的結果容器，先統一用 border_threshold 處理整塊\n        processed_block = (blurred_sub > border_threshold).astype(vol.dtype)\n\n        # 定義中心區域 (Inner area) 的切片\n        inner_z = slice(border, -border if border > 0 else None)\n        inner_h = slice(border, -border if border > 0 else None)\n        inner_w = slice(border, -border if border > 0 else None)\n        inner_slice = (inner_z, inner_h, inner_w)\n\n        # 針對「中心區域」覆蓋使用 inner_threshold 的結果\n        # 這會讓中心區域與邊界區域有不同的二值化敏感度\n        processed_block[inner_slice] = (blurred_sub[inner_slice] > inner_threshold).astype(vol.dtype)\n\n        # 將處理完的塊填回大圖\n        output_vol[sl] = processed_block\n\n    return output_vol\n\n\ndef apply_post_processing(\n        pred: np.ndarray,\n        skip: bool,\n        min_size: int,\n        small_holes_size: int,\n        split_parer: bool,\n        split_paper_iter: int = 1,\n        prob: np.ndarray = None,  # 機率圖\n        fg_threshold: float = 0.3,  # 參數\n        fg_bg_diff_threshold: float = 0.3,  # 參數\n        line_norm: bool = False,\n        y_axis_closing: bool = False,\n        z_axis_closing: bool = False,\n        closing: bool = False,\n        y_axis_fill_holes: bool = False,\n        z_axis_fill_holes: bool = False,\n        fill_hole: bool = False,\n        line_endpoints_repair: bool = False,\n        line_mask_repair: bool = False,\n        repair_settings: dict = None,\n        gaussian: bool = False,\n        sandwich: bool = False,\n        diag: bool = False,\n) -> np.ndarray:\n    \"\"\"\n    整合後處理流程：\n    1. 機率弱訊號恢復 (Weak Signal Recovery)\n    2. 移除小物件 (Remove Small Objects)\n    3. 斷線修復 (Line Repair)\n    4. 形態學優化 (Morphological Ops)\n    \"\"\"\n    # 複製並初始化\n    vol = pred.copy()\n    vol[vol == 2] = 0  # 清除 Class 2\n\n    if skip:\n        return vol\n\n    if prob is not None:\n        vol = recover_weak_signals(\n            vol,\n            prob,\n            fg_threshold=fg_threshold,\n            fg_bg_diff_threshold=fg_bg_diff_threshold\n        )\n    else:\n        vol = (vol > 0).astype(np.uint8)\n    if min_size > 0:\n        vol_bool = vol > 0\n        vol_bool = remove_small_objects(vol_bool, min_size=min_size, connectivity=1)\n        vol = vol_bool.astype(np.uint8)\n    if line_endpoints_repair:\n        vol = execute_line_repair(vol, settings=repair_settings)\n        vol = (vol > 0).astype(np.uint8)\n    if line_mask_repair:\n        vol, _ = vol | line_repair_msk(vol_bool=vol, pass_iters=1, radius=1.5, axis=0, max_dist=30, max_angle_deg=30,\n                                       safe_distance=3, border=8)\n        # vol = line_repair_msk(vol_bool=vol, pass_iters=1, radius=1.5, axis=1, max_dist=30, max_angle_deg=30)\n        vol = (vol > 0).astype(np.uint8)\n    if line_norm:\n        Z_15, _ = normalize_segments_3d(vol, radius=1.5, border=8, axis=0, sigma=0, repair=True, repair_dijkstra=False,\n                                        max_dist=30, max_angle_deg=45, safe_distance=2, iterations=1)\n        # Z_20, _ = normalize_segments_3d(vol, radius=2, border=8, axis=0, sigma=0, repair=False, repair_dijkstra=False,\n        #                                 max_dist=30, max_angle_deg=45, safe_distance=2, iterations=1)\n        Y_15, _ = normalize_segments_3d(vol, radius=1.5, border=8, axis=1, sigma=0, repair=True, repair_dijkstra=False,\n                                        max_dist=30, max_angle_deg=45, safe_distance=2, iterations=1)\n        # Y_20, _ = normalize_segments_3d(vol, radius=2, border=8, axis=1, sigma=0, repair=False, repair_dijkstra=False,\n        #                                 max_dist=30, max_angle_deg=45, safe_distance=2, iterations=1)\n        # X_15, _ = normalize_segments_3d(vol, radius=1.5, border=8, axis=2, sigma=0, repair=True, repair_dijkstra=False,\n        #                                 max_dist=30, max_angle_deg=45, safe_distance=2, iterations=1)\n        # vol = vol | Z_15\n        vol = vol | Z_15 | Y_15\n        # vol = (vol | Z_15) & (Z_20 | Y_20)\n        vol = vol.astype(np.uint8)\n    if sandwich:\n        vol = robust_mask_refine(vol, max_gap=3, iterations=5, min_neighbors=9, kernel_size=5, keep_connectivity=False)\n        vol = vol.astype(np.uint8)\n    if split_parer:\n        vol_bool = vol > 0\n        vol_bool = remove_small_objects(vol_bool, min_size=min_size, connectivity=1)\n        vol = vol_bool.astype(np.uint8)\n        vol = get_split_paper(vol, ray_len=64, max_iter=split_paper_iter, cleanup_iter=2)\n        vol = vol.astype(np.uint8)\n    if y_axis_closing:\n        vol = apply_y_axis_closing(vol, iterations=1)\n        vol = vol.astype(np.uint8)\n    if z_axis_closing:\n        vol = apply_z_axis_closing(vol, iterations=1)\n        vol = vol.astype(np.uint8)\n    if closing:\n        vol = apply_closing(vol, iterations=1, conn=3)\n        vol = vol.astype(np.uint8)\n\n    if z_axis_fill_holes:\n        vol = apply_z_axis_fill_holes(vol)\n        vol = vol.astype(np.uint8)\n    if y_axis_fill_holes:\n        vol = apply_y_axis_fill_holes(vol)\n        vol = vol.astype(np.uint8)\n    if fill_hole:\n        structure_6_connected = np.array([\n            [[0, 0, 0], [0, 1, 0], [0, 0, 0]],\n            [[0, 1, 0], [1, 1, 1], [0, 1, 0]],\n            [[0, 0, 0], [0, 1, 0], [0, 0, 0]]\n        ])\n        structure_18_connected = np.array([\n            [[0, 1, 0], [1, 1, 1], [0, 1, 0]],\n            [[1, 1, 1], [1, 1, 1], [1, 1, 1]],\n            [[0, 1, 0], [1, 1, 1], [0, 1, 0]]\n        ])\n        structure_26_connected = np.array([\n            [[1, 1, 1], [1, 1, 1], [1, 1, 1]],\n            [[1, 1, 1], [1, 1, 1], [1, 1, 1]],\n            [[1, 1, 1], [1, 1, 1], [1, 1, 1]]\n        ])\n        vol = binary_fill_holes(vol, structure=structure_6_connected)\n        vol = vol.astype(np.uint8)\n\n    if gaussian:\n        vol = process_gaussian_octants_dual_threshold(vol=vol, inner_threshold=0.5, border_threshold=0.5,\n                                                      sigma=(1.0, 1.0, 1.0), border=8, mode=\"reflect\", cval=0.0)\n        # mode = \"mirror\"  # reflect constant nearest mirror wrap\n    if min_size > 0:\n        vol_bool = vol > 0\n        vol_bool = remove_small_objects(vol_bool, min_size=min_size, connectivity=1)\n        vol_bool = remove_small_8slice(vol_bool, min_size=10)\n        vol = vol_bool.astype(np.uint8)\n    if small_holes_size > 0:\n        vol_bool = vol > 0\n        vol_bool = remove_small_holes(vol_bool, area_threshold=small_holes_size)\n        vol = vol_bool.astype(np.uint8)\n    if diag:\n        vol = fill_2x2diag(vol, axis=0, gap=0)\n        vol = fill_2x2diag(vol, axis=0, gap=1)\n        vol = fill_2x2diag(vol, axis=1, gap=0)\n        vol = fill_2x2diag(vol, axis=1, gap=1)\n        vol = vol.astype(np.uint8)\n    return vol\n\n\nif __name__ == \"__main__\":\n    a = np.array(([1, 1, 1, 1, 1, 1, 1],\n                  [1, 1, 1, 1, 1, 1, 1],\n                  [1, 1, 1, 1, 1, 1, 1],\n                  [1, 1, 1, 0, 1, 1, 1],\n                  [1, 1, 1, 1, 1, 1, 1],\n                  [1, 1, 1, 1, 1, 1, 1],\n                  [1, 1, 1, 1, 1, 1, 1]))\n    np.set_printoptions(precision=3)\n    print(distance_transform_edt(a))\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Post-Processing","metadata":{}},{"cell_type":"code","source":"from pathlib import Path\nimport tifffile as tiff\nfrom skimage import measure\nimport zipfile\nimport csv\n\n# ==========================================\n# 2. 主程式設定與執行 (Main Execution)\n# ==========================================\n\n# --- 設定路徑 ---\nsource_dir = '/kaggle/working/nnUNet_upload/3d_fullres_TTA'\nprocessed_dir = '/kaggle/working/processed_labels'\noutput_zip = '/kaggle/working/submission.zip'\nstats_file = '/kaggle/working/processing_stats.csv'\n\nos.makedirs(processed_dir, exist_ok=True)\n\n# --- 設定後處理參數 ---\nMIN_SIZE = 5000         # 移除小於此體積 (voxel數) 的雜訊\n\n\n# 準備統計數據列表\nstats_data = []\n\nprint(f\"{'Filename':<30} | {'Before':<10} | {'After':<10} | {'Diff':<10}\")\nprint(\"-\" * 70)\n\ntif_files = list(Path(source_dir).glob('*.tif'))\n\nfor tif_path in tif_files:\n    # A. 讀取數據\n    data = tiff.imread(str(tif_path))\n    \n    # 確保是二值化 (0 或 1)\n    binary_mask = (data > 0).astype(np.uint8)\n    \n    # --- 統計處理前物件數量 ---\n    # label() 用於計算連通域數量\n    _, count_before = measure.label(binary_mask, return_num=True)\n    \n    # B. 執行新的後處理函數\n    SKIP = False\n    min_size = 5000\n    small_holes_size = 0\n    split_parer = False\n    split_paper_iter = 1\n    line_norm = True\n    line_mask_repair = False\n    line_endpoints_repair = False\n    y_axis_closing = False\n    z_axis_closing = False\n    closing = False\n    z_axis_fill_holes = False\n    y_axis_fill_holes = False\n    fill_hole = False\n    gaussian = True\n    sandwich = True\n    diag = True\n\n    processed_data = apply_post_processing(\n        binary_mask, \n        skip=SKIP,\n        min_size=min_size,\n        small_holes_size=small_holes_size,\n        split_parer=split_parer,\n        split_paper_iter=split_paper_iter,\n        line_norm=line_norm,\n        line_mask_repair=line_mask_repair,\n        line_endpoints_repair=line_endpoints_repair,\n        y_axis_closing=y_axis_closing,\n        closing=closing,\n        z_axis_fill_holes=z_axis_fill_holes,\n        y_axis_fill_holes=y_axis_fill_holes,\n        fill_hole=fill_hole,\n        gaussian=gaussian,\n        sandwich=sandwich,\n        diag=diag,\n    )\n    \n    # --- 統計處理後物件數量 ---\n    _, count_after = measure.label(processed_data, return_num=True)\n        \n    # 記錄統計數據\n    diff = count_after - count_before\n    stats_data.append([tif_path.name, count_before, count_after, diff])\n    \n    # C. 儲存處理後的圖檔 (只有 0 與 1)\n    save_path = os.path.join(processed_dir, tif_path.name)\n    tiff.imwrite(save_path, processed_data, compression='zlib')\n    \n    print(f\"{tif_path.name:<30} | {count_before:<10} | {count_after:<10} | {diff:<+10}\")\n\n# ==========================================\n# 3. 匯出結果 (ZIP & CSV)\n# ==========================================\n\n# 建立 ZIP 檔案\nwith zipfile.ZipFile(output_zip, 'w', zipfile.ZIP_DEFLATED) as zf:\n    for tif_file in Path(processed_dir).glob('*.tif'):\n        zf.write(tif_file, arcname=tif_file.name)\n\n# 輸出 CSV 統計報表\nwith open(stats_file, 'w', newline='') as csvfile:\n    writer = csv.writer(csvfile)\n    writer.writerow(['Filename', 'Count_Before', 'Count_After', 'Difference'])\n    writer.writerows(stats_data)\n\nprint(\"\\n\" + \"=\"*70)\nprint(f\"處理完成！使用新演算法\")\nprint(f\"總處理檔案數: {len(tif_files)}\")\nprint(f\"ZIP 檔案路徑: {output_zip}\")\nprint(f\"CSV 統計報表: {stats_file}\")\nprint(\"=\"*70)","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}