{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":117682,"databundleVersionId":15062069,"sourceType":"competition"},{"sourceId":290917305,"sourceType":"kernelVersion"},{"sourceId":655294,"sourceType":"modelInstanceVersion","modelInstanceId":495238,"modelId":510647},{"sourceId":660383,"sourceType":"modelInstanceVersion","modelInstanceId":499479,"modelId":510647},{"sourceId":665589,"sourceType":"modelInstanceVersion","modelInstanceId":503784,"modelId":510647},{"sourceId":665924,"sourceType":"modelInstanceVersion","modelInstanceId":504051,"modelId":510647},{"sourceId":672178,"sourceType":"modelInstanceVersion","modelInstanceId":495238,"modelId":510647},{"sourceId":673516,"sourceType":"modelInstanceVersion","modelInstanceId":499479,"modelId":510647},{"sourceId":674747,"sourceType":"modelInstanceVersion","modelInstanceId":503784,"modelId":510647},{"sourceId":681152,"sourceType":"modelInstanceVersion","modelInstanceId":516822,"modelId":510647},{"sourceId":726582,"sourceType":"modelInstanceVersion","modelInstanceId":553092,"modelId":565646}],"dockerImageVersionId":31193,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Training Notebooks\n\n- [Vesuvius Surface 3D Detection in Keras-JAX](https://www.kaggle.com/code/ipythonx/vesuvius-surface-3d-detection-in-jax)\n- [Vesuvius Surface 3D Detection in PyTorch](https://www.kaggle.com/code/ipythonx/vesuvius-surface-3d-detection-in-pytorch)\n- [Vesuvius Surface 3D Detection in PyTorch Lightning](https://www.kaggle.com/code/ipythonx/train-vesuvius-surface-3d-detection-in-lightning)\n- [[WIP] Vesuvius Surface 2.5D Detection](https://www.kaggle.com/code/ipythonx/wip-vesuvius-surface-2-5d-detection)\n\n**Note**\n1. The inference code below is adapted from the **Keras-JAX** version. The PyTorch and Lightning implementations follow the same workflow. Training was performed on a single Tesla T4 (16 GB VRAM) with extended epochs.\n2. Both the training and inference pipelines are implemented using [`medicai`](https://github.com/innat/medic-ai), a **Keras 3** based multi-backend medical ML library designed for 2D and 3D classification and segmentation tasks. However, please note, `medicai` project is still new and actively evolving.","metadata":{}},{"cell_type":"markdown","source":"# Inference","metadata":{}},{"cell_type":"code","source":"var=\"/kaggle/input/vsdetection-packages-offline-installer-only/whls\"\n!pip install \\\n    \"$var\"/keras_nightly-3.12.0.dev2025100703-py3-none-any.whl \\\n    \"$var\"/tifffile-2025.12.20-py3-none-any.whl \\\n    \"$var\"/imagecodecs-2026.1.1-cp311-abi3-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl \\\n    \"$var\"/medicai-0.0.3-py3-none-any.whl \\\n    --no-index \\\n    --find-links \"$var\"","metadata":{"trusted":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2026-01-27T04:08:34.220862Z","iopub.execute_input":"2026-01-27T04:08:34.221566Z","iopub.status.idle":"2026-01-27T04:08:42.081772Z","shell.execute_reply.started":"2026-01-27T04:08:34.221542Z","shell.execute_reply":"2026-01-27T04:08:42.081070Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# --- Protobuf compatibility patch (for old code using MessageFactory.GetPrototype) ---\n\ntry:\n    from google.protobuf import message_factory as _message_factory\n\n    # Only patch if the method is missing (protobuf >= 5)\n    if not hasattr(_message_factory.MessageFactory, \"GetPrototype\"):\n        from google.protobuf.message_factory import GetMessageClass\n\n        def _GetPrototype(self, descriptor):\n            # Old API used MessageFactory().GetPrototype(descriptor)\n            # New API is GetMessageClass(descriptor). We just bridge them.\n            return GetMessageClass(descriptor)\n\n        _message_factory.MessageFactory.GetPrototype = _GetPrototype\n        print(\"Patched protobuf: added MessageFactory.GetPrototype\")\n    else:\n        print(\"protobuf already has MessageFactory.GetPrototype; no patch needed.\")\nexcept Exception as e:\n    print(\"Could not patch protobuf MessageFactory:\", e)\n\nimport os\nos.environ[\"KERAS_BACKEND\"] = \"jax\"\n\nimport keras\nfrom medicai.transforms import (\n    Compose,\n    ScaleIntensityRange,\n)\nfrom medicai.models import SegFormer, TransUNet\nfrom medicai.utils.inference import SlidingWindowInference\n\nimport numpy as np\nimport pandas as pd\nimport zipfile\nimport tifffile\nfrom matplotlib import pyplot as plt\n\nkeras.config.backend(), keras.version()","metadata":{"trusted":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2026-01-27T04:08:42.083404Z","iopub.execute_input":"2026-01-27T04:08:42.083831Z","iopub.status.idle":"2026-01-27T04:08:57.953075Z","shell.execute_reply.started":"2026-01-27T04:08:42.083806Z","shell.execute_reply":"2026-01-27T04:08:57.952416Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"**Dataset**","metadata":{}},{"cell_type":"code","source":"root_dir = \"/kaggle/input/vesuvius-challenge-surface-detection\"\ntest_dir = f\"{root_dir}/test_images\"\noutput_dir = \"/kaggle/working/submission_masks\"\nzip_path = \"/kaggle/working/submission.zip\"\nos.makedirs(output_dir, exist_ok=True)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-27T04:08:57.953808Z","iopub.execute_input":"2026-01-27T04:08:57.954303Z","iopub.status.idle":"2026-01-27T04:08:57.958739Z","shell.execute_reply.started":"2026-01-27T04:08:57.954275Z","shell.execute_reply":"2026-01-27T04:08:57.957839Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"test_df = pd.read_csv(f\"{root_dir}/test.csv\")\ntest_df.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-27T04:08:57.960692Z","iopub.execute_input":"2026-01-27T04:08:57.961052Z","iopub.status.idle":"2026-01-27T04:08:58.081620Z","shell.execute_reply.started":"2026-01-27T04:08:57.961027Z","shell.execute_reply":"2026-01-27T04:08:58.080768Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"**Transformation**","metadata":{}},{"cell_type":"code","source":"def val_transformation(image):\n    data = {\"image\": image}\n    pipeline = Compose([\n        ScaleIntensityRange(\n            keys=[\"image\"],\n            a_min = 0,\n            a_max = 255,\n            b_min = 0,\n            b_max = 1,\n            clip = True,\n        ),\n    ])\n    result = pipeline(data)\n    return result[\"image\"]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-27T04:08:58.082321Z","iopub.execute_input":"2026-01-27T04:08:58.082569Z","iopub.status.idle":"2026-01-27T04:08:58.088502Z","shell.execute_reply.started":"2026-01-27T04:08:58.082551Z","shell.execute_reply":"2026-01-27T04:08:58.087535Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"**Model**","metadata":{}},{"cell_type":"code","source":"num_classes=3\n\ndef get_model():\n    ## LB: 0.486\n    # model = SegFormer(\n    #     input_shape=(128, 128, 128, 1),\n    #     encoder_name='mit_b2',\n    #     classifier_activation='softmax',\n    #     num_classes=num_classes,\n    # )\n    # model.load_weights(\n    #     \"/kaggle/input/vsd-model/keras/segformer.mit.b2/2/segformer.mit.b2.weights.h5\"\n    # )\n\n    ## LB: 0.5 \n    model = TransUNet(\n        input_shape=(160, 160, 160, 1),\n        encoder_name='seresnext50',\n        classifier_activation='softmax',\n        num_classes=num_classes,\n    )\n    model.load_weights(\n        \"/kaggle/input/opentransunet/pytorch/default/1/model.weights.h5\"\n    )\n    return model","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-27T04:08:58.089319Z","iopub.execute_input":"2026-01-27T04:08:58.089646Z","iopub.status.idle":"2026-01-27T04:08:58.104643Z","shell.execute_reply.started":"2026-01-27T04:08:58.089627Z","shell.execute_reply":"2026-01-27T04:08:58.103963Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"model = get_model()\nmodel.count_params() / 1e6\n# predictor = tf.function(model, jit_compile=True)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-27T04:08:58.105822Z","iopub.execute_input":"2026-01-27T04:08:58.106094Z","iopub.status.idle":"2026-01-27T04:09:18.546833Z","shell.execute_reply.started":"2026-01-27T04:08:58.106075Z","shell.execute_reply":"2026-01-27T04:09:18.546170Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"model.instance_describe()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-27T04:09:18.547589Z","iopub.execute_input":"2026-01-27T04:09:18.547850Z","iopub.status.idle":"2026-01-27T04:09:18.614256Z","shell.execute_reply.started":"2026-01-27T04:09:18.547830Z","shell.execute_reply":"2026-01-27T04:09:18.613654Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"model","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-27T04:09:18.615036Z","iopub.execute_input":"2026-01-27T04:09:18.615325Z","iopub.status.idle":"2026-01-27T04:09:18.620662Z","shell.execute_reply.started":"2026-01-27T04:09:18.615303Z","shell.execute_reply":"2026-01-27T04:09:18.619858Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"**Sliding Window Inference**","metadata":{}},{"cell_type":"code","source":"pred = SlidingWindowInference(\n    model,\n    roi_size=(160,160,160),\n    num_classes = 3,\n    mode=\"gaussian\",\n    overlap=0.6,\n    sw_batch_size = 1\n)\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-27T04:09:18.622827Z","iopub.execute_input":"2026-01-27T04:09:18.623077Z","iopub.status.idle":"2026-01-27T04:09:18.632136Z","shell.execute_reply.started":"2026-01-27T04:09:18.623060Z","shell.execute_reply":"2026-01-27T04:09:18.631376Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nimport numpy as np\nimport tensorflow as tf\nfrom skimage.morphology import skeletonize\nfrom scipy.ndimage import binary_dilation\nfrom scipy import ndimage\n\ndef frangi_filter_3d(image, sigmas=(1, 3, 5), beta1=0.5, beta2=15, gamma=None, black_ridges=True):\n    \"\"\"\n    Frangi滤波器 - 用于增强3D图像中的管状结构（如血管）。\n    \n    Frangi滤波器基于Hessian矩阵的特征值来检测管状结构。它通过分析局部二阶导数\n    来识别具有\"管状\"几何特征的区域。\n    \n    原理：\n    1. 计算图像在不同尺度下的Hessian矩阵\n    2. 计算Hessian矩阵的特征值（λ1, λ2, λ3）\n    3. 根据特征值的关系判断是否为管状结构\n    4. 使用Frangi响应函数计算增强值\n    \n    参数:\n        image: numpy数组，形状为 (D, H, W) 或 (D, H, W, 1) 的3D图像\n        sigmas: tuple，高斯核的标准差范围，用于多尺度检测\n                例如 (1, 3, 5) 表示检测1到5像素宽度的血管\n        beta1: float，控制对偏离管状结构的敏感度（默认0.5）\n        beta2: float，控制对背景的敏感度（默认15）\n        gamma: float，归一化因子，如果为None则自动计算\n        black_ridges: bool，True表示检测暗色管状结构（血管），False表示检测亮色结构\n    \n    返回:\n        enhanced: numpy数组，增强后的图像，形状与输入相同\n    \"\"\"\n    # 确保输入是3D数组\n    if len(image.shape) == 4:\n        image = image[..., 0]\n    \n    # 确保是float类型\n    image = image.astype(np.float64)\n    \n    # 如果检测暗色结构，反转图像\n    if black_ridges:\n        image = -image\n    \n    # 初始化输出\n    enhanced = np.zeros_like(image)\n    \n    # 对每个尺度计算Frangi响应\n    for sigma in sigmas:\n        # 计算Hessian矩阵的各个分量\n        # 使用高斯滤波的导数\n        hxx = ndimage.gaussian_filter1d(\n            ndimage.gaussian_filter1d(\n                ndimage.gaussian_filter1d(image, sigma, axis=2, order=2), \n                sigma, axis=1, order=0), \n            sigma, axis=0, order=0)\n        \n        hyy = ndimage.gaussian_filter1d(\n            ndimage.gaussian_filter1d(\n                ndimage.gaussian_filter1d(image, sigma, axis=2, order=0), \n                sigma, axis=1, order=2), \n            sigma, axis=0, order=0)\n        \n        hzz = ndimage.gaussian_filter1d(\n            ndimage.gaussian_filter1d(\n                ndimage.gaussian_filter1d(image, sigma, axis=2, order=0), \n                sigma, axis=1, order=0), \n            sigma, axis=0, order=2)\n        \n        hxy = ndimage.gaussian_filter1d(\n            ndimage.gaussian_filter1d(\n                ndimage.gaussian_filter1d(image, sigma, axis=2, order=1), \n                sigma, axis=1, order=1), \n            sigma, axis=0, order=0)\n        \n        hxz = ndimage.gaussian_filter1d(\n            ndimage.gaussian_filter1d(\n                ndimage.gaussian_filter1d(image, sigma, axis=2, order=1), \n                sigma, axis=1, order=0), \n            sigma, axis=0, order=1)\n        \n        hyz = ndimage.gaussian_filter1d(\n            ndimage.gaussian_filter1d(\n                ndimage.gaussian_filter1d(image, sigma, axis=2, order=0), \n                sigma, axis=1, order=1), \n            sigma, axis=0, order=1)\n        \n        # 构建Hessian矩阵并计算特征值\n        # 对于3D，我们需要在每个体素处计算3x3矩阵的特征值\n        # 这里使用近似方法：计算特征值的平方和\n        # 更精确的方法需要逐像素计算特征值，但计算量大\n        \n        # 简化的Frangi响应计算（基于Hessian矩阵的迹和行列式）\n        # 对于管状结构：一个特征值接近0，另外两个较大且符号相同\n        trace = hxx + hyy + hzz\n        det = (hxx * hyy * hzz + \n               2 * hxy * hxz * hyz - \n               hxx * hyz * hyz - \n               hyy * hxz * hxz - \n               hzz * hxy * hxy)\n        \n        # Frangi响应函数\n        # 对于管状结构，我们希望：\n        # - 一个特征值接近0（沿管轴方向）\n        # - 另外两个特征值较大且同号（垂直于管轴）\n        \n        # 使用迹和行列式的组合来近似\n        # 归一化\n        if gamma is None:\n            gamma = np.max(np.abs(trace))\n        \n        # 计算响应\n        # 简化版本：基于Hessian矩阵的Frobenius范数\n        hessian_norm = np.sqrt(hxx**2 + hyy**2 + hzz**2 + \n                              2 * (hxy**2 + hxz**2 + hyz**2))\n        \n        # Frangi响应\n        response = np.exp(-beta1 * (trace**2) / (hessian_norm + 1e-10)) * \\\n                   (1 - np.exp(-beta2 * hessian_norm**2 / (gamma**2 + 1e-10)))\n        \n        # 只保留正值（管状结构）\n        response = np.maximum(response, 0)\n        \n        # 取所有尺度中的最大值\n        enhanced = np.maximum(enhanced, response)\n    \n    # 归一化到[0, 1]\n    if enhanced.max() > 0:\n        enhanced = enhanced / enhanced.max()\n    \n    return enhanced.astype(np.float32)\n\n\ndef apply_frangi_postprocessing(prediction, threshold=0.5, sigmas=(1, 3, 5), \n                                beta1=0.5, beta2=15, use_frangi=True):\n    \"\"\"\n    使用Frangi滤波器对分割预测结果进行后处理。\n    \n    这个方法可以：\n    1. 使用Frangi滤波器增强预测结果中的管状结构\n    2. 结合原始预测和Frangi增强结果\n    3. 应用阈值得到最终的分割掩码\n    \n    参数:\n        prediction: numpy数组，模型预测的概率图，形状为 (D, H, W) 或 (D, H, W, num_classes)\n        threshold: float，二值化阈值（默认0.5）\n        sigmas: tuple，Frangi滤波器的尺度参数\n        beta1: float，Frangi参数\n        beta2: float，Frangi参数\n        use_frangi: bool，是否使用Frangi滤波器（如果False，只做阈值化）\n    \n    返回:\n        processed_mask: numpy数组，处理后的二值掩码，形状为 (D, H, W)\n    \"\"\"\n    # 处理多类别预测\n    if len(prediction.shape) == 4:\n        # 如果是多类别，取前景类别（假设类别1是血管/墨水）\n        if prediction.shape[-1] > 1:\n            pred_prob = prediction[..., 1]  # 取类别1的概率\n        else:\n            pred_prob = prediction[..., 0]\n    else:\n        pred_prob = prediction\n    \n    if use_frangi:\n        # 应用Frangi滤波器增强管状结构\n        frangi_enhanced = frangi_filter_3d(\n            pred_prob, \n            sigmas=sigmas, \n            beta1=beta1, \n            beta2=beta2,\n            black_ridges=True  # 假设检测暗色结构\n        )\n        \n        # 结合原始预测和Frangi增强结果\n        # 方法1：加权平均\n        combined = 0.7 * pred_prob + 0.3 * frangi_enhanced\n        \n        # 方法2：取最大值（保留更强的响应）\n        # combined = np.maximum(pred_prob, frangi_enhanced)\n        \n        # 方法3：只在Frangi检测到的区域增强\n        # frangi_mask = frangi_enhanced > 0.3\n        # combined = pred_prob.copy()\n        # combined[frangi_mask] = np.maximum(pred_prob[frangi_mask], frangi_enhanced[frangi_mask])\n        \n        processed_prob = combined\n    else:\n        processed_prob = pred_prob\n    \n    # 二值化\n    processed_mask = (processed_prob > threshold).astype(np.uint8)\n    \n    return processed_mask\n\ndef build_anisotropic_struct(z_radius, xy_radius):\n    \"\"\"\n    构建各向异性结构元素用于3D闭运算\n    \n    参数:\n        z_radius: z方向的半径\n        xy_radius: xy平面的半径\n    \n    返回:\n        structure: 3D结构元素，如果半径都为0则返回None\n    \"\"\"\n    if z_radius == 0 and xy_radius == 0:\n        return None\n    \n    size_z = 2 * z_radius + 1\n    size_xy = 2 * xy_radius + 1\n    \n    struct = np.zeros((size_z, size_xy, size_xy), dtype=bool)\n    center_z = z_radius\n    center_xy = xy_radius\n    \n    # 创建各向异性结构\n    for z in range(size_z):\n        for y in range(size_xy):\n            for x in range(size_xy):\n                dz = abs(z - center_z)\n                dxy = np.sqrt((x - center_xy)**2 + (y - center_xy)**2)\n                if dz <= z_radius and dxy <= xy_radius:\n                    struct[z, y, x] = True\n    \n    return struct","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-27T04:09:18.633030Z","iopub.execute_input":"2026-01-27T04:09:18.633369Z","iopub.status.idle":"2026-01-27T04:09:18.776010Z","shell.execute_reply.started":"2026-01-27T04:09:18.633329Z","shell.execute_reply":"2026-01-27T04:09:18.775386Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nimport tifffile\nimport scipy.ndimage as ndi\nfrom skimage.morphology import remove_small_objects\nfrom scipy.ndimage import generate_binary_structure\nfrom skimage.morphology import remove_small_objects\n\ndef load_volume(path):\n    vol = tifffile.imread(path)          # (D, H, W)\n    vol = vol.astype(np.float32)\n    vol = vol[None, ..., None]           # (1, D, H, W, 1)\n    return vol\n\n\n# ==========================================\n# ROTATION TTA HELPERS (CLOCKWISE)\n# ==========================================\ndef rot90_volume(vol, k):\n    \"\"\"\n    Rotate volume k times 90° clockwise in HW plane.\n    vol:\n      (1, D, H, W, 1) OR (D, H, W)\n    \"\"\"\n    if vol.ndim == 5:\n        return np.rot90(vol, k=-k, axes=(2, 3))\n    else:\n        return np.rot90(vol, k=-k, axes=(1, 2))\n\n\ndef unrot90_volume(vol, k):\n    return rot90_volume(vol, (4 - k) % 4)\n\n\ndef predict_probs_tta_rot(sample):\n    \"\"\"\n    4x rotation TTA: 0°, 90°, 180°, 270°\n    sample: (1, D, H, W, 1)\n    returns: averaged probs (D, H, W)\n    \"\"\"\n    probs_accum = []\n\n    for k in range(4):\n        s_rot = rot90_volume(sample, k)\n\n        out = pred(s_rot)              # (1, D, H, W, 2)\n        out = np.asarray(out)\n        probs = out[0, ..., 1]         # (D, H, W)\n\n        probs = unrot90_volume(probs, k)\n        probs_accum.append(probs)\n\n    return np.mean(probs_accum, axis=0)\n\n\n# ==========================================\n# HELPER: Anisotropic Structure Builder\n# ==========================================\ndef build_anisotropic_struct(z_radius: int, xy_radius: int):\n    z, r = z_radius, xy_radius\n\n    if z == 0 and r == 0:\n        return None\n\n    if z == 0 and r > 0:\n        size = 2 * r + 1\n        struct = np.zeros((1, size, size), dtype=bool)\n        cy, cx = r, r\n        for dy in range(-r, r + 1):\n            for dx in range(-r, r + 1):\n                if dy * dy + dx * dx <= r * r:\n                    struct[0, cy + dy, cx + dx] = True\n        return struct\n\n    if z > 0 and r == 0:\n        struct = np.zeros((2 * z + 1, 1, 1), dtype=bool)\n        struct[:, 0, 0] = True\n        return struct\n\n    depth = 2 * z + 1\n    size = 2 * r + 1\n    struct = np.zeros((depth, size, size), dtype=bool)\n    cz, cy, cx = z, r, r\n    for dz in range(-z, z + 1):\n        for dy in range(-r, r + 1):\n            for dx in range(-r, r + 1):\n                if dy * dy + dx * dx <= r * r:\n                    struct[cz + dz, cy + dy, cx + dx] = True\n    return struct\n\n\n# ==========================================\n# MAIN POST-PROCESSING LOGIC\n# ==========================================\ndef topo_postprocess(\n    probs,          # (D, H, W)\n    T_low=0.90,\n    T_high=0.90,\n    z_radius=1,\n    xy_radius=0,\n    dust_min_size=100,\n):\n    # --- Step 1: 3D Hysteresis ---\n    strong = probs >= T_high\n    weak   = probs >= T_low\n\n    if not strong.any():\n        return np.zeros_like(probs, dtype=np.uint8)\n\n    struct_hyst = ndi.generate_binary_structure(3, 3)\n    mask = ndi.binary_propagation(strong, mask=weak, structure=struct_hyst)\n\n    if not mask.any():\n        return np.zeros_like(probs, dtype=np.uint8)\n\n    # --- Step 2: 3D Anisotropic Closing ---\n    if z_radius > 0 or xy_radius > 0:\n        struct_close = build_anisotropic_struct(z_radius, xy_radius)\n        if struct_close is not None:\n            mask = ndi.binary_closing(mask, structure=struct_close)\n\n    # --- Step 3: Dust Removal ---\n    if dust_min_size > 0:\n        mask = remove_small_objects(mask.astype(bool), min_size=dust_min_size)\n\n    return mask.astype(np.uint8)\n\n\ndef topo_postprocess_with_frangi(\n    probs,          # (D, H, W)\n    T_low=0.90,\n    T_high=0.90,\n    z_radius=1,\n    xy_radius=0,\n    dust_min_size=100,\n    # Frangi参数\n    use_frangi=True,\n    frangi_weight=0.3,      # Frangi增强的权重 (0-1)\n    frangi_sigmas=(1, 3, 5),\n    frangi_beta1=0.5,\n    frangi_beta2=15,\n    frangi_mode='enhance_before',  # 'enhance_before', 'enhance_after', 'replace_hyst'\n):\n    \"\"\"\n    集成Frangi滤波器的拓扑后处理\n    \n    参数:\n        probs: 概率图 (D, H, W)\n        T_low, T_high: Hysteresis阈值\n        z_radius, xy_radius: 各向异性闭运算参数\n        dust_min_size: 小对象去除的最小尺寸\n        use_frangi: 是否使用Frangi滤波器\n        frangi_weight: Frangi增强的权重（0-1之间）\n        frangi_sigmas: Frangi滤波器的尺度参数\n        frangi_beta1, frangi_beta2: Frangi参数\n        frangi_mode: Frangi集成模式\n            - 'enhance_before': 在Hysteresis之前增强概率图（推荐）\n            - 'enhance_after': 在Hysteresis之后增强掩码\n            - 'replace_hyst': 用Frangi增强的概率图替代Hysteresis\n    \n    返回:\n        mask: 处理后的二值掩码 (D, H, W) uint8\n    \"\"\"\n    if not use_frangi:\n        # 如果不使用Frangi，直接使用原始后处理\n        return topo_postprocess(\n            probs, T_low, T_high, z_radius, xy_radius, dust_min_size\n        )\n    \n    # 应用Frangi滤波器增强\n    frangi_enhanced = frangi_filter_3d(\n        probs,\n        sigmas=frangi_sigmas,\n        beta1=frangi_beta1,\n        beta2=frangi_beta2,\n        black_ridges=True  # 假设检测暗色结构（血管）\n    )\n    \n    if frangi_mode == 'enhance_before':\n        # 模式1: 在Hysteresis之前增强概率图（推荐）\n        # 结合原始概率和Frangi增强结果\n        enhanced_probs = (1 - frangi_weight) * probs + frangi_weight * frangi_enhanced\n        \n        # 然后进行标准的拓扑后处理\n        return topo_postprocess(\n            enhanced_probs, T_low, T_high, z_radius, xy_radius, dust_min_size\n        )\n    \n    elif frangi_mode == 'enhance_after':\n        # 模式2: 在Hysteresis之后增强掩码\n        # 先进行标准拓扑后处理\n        mask = topo_postprocess(\n            probs, T_low, T_high, z_radius, xy_radius, dust_min_size\n        )\n        \n        # 将掩码转换为概率图（用于Frangi）\n        mask_probs = mask.astype(np.float32)\n        \n        # 应用Frangi增强\n        frangi_mask = frangi_filter_3d(\n            mask_probs,\n            sigmas=frangi_sigmas,\n            beta1=frangi_beta1,\n            beta2=frangi_beta2,\n            black_ridges=True\n        )\n        \n        # 结合原始掩码和Frangi增强结果\n        enhanced_mask = np.maximum(\n            mask_probs,\n            frangi_mask * frangi_weight\n        )\n        \n        # 二值化\n        return (enhanced_mask > 0.5).astype(np.uint8)\n    \n    elif frangi_mode == 'replace_hyst':\n        # 模式3: 用Frangi增强的概率图替代Hysteresis\n        # 结合原始概率和Frangi增强结果\n        enhanced_probs = (1 - frangi_weight) * probs + frangi_weight * frangi_enhanced\n        \n        # 使用简单的阈值化替代Hysteresis\n        mask = enhanced_probs >= T_high\n        \n        if not mask.any():\n            return np.zeros_like(probs, dtype=np.uint8)\n        \n        # 然后进行闭运算和小对象去除\n        if z_radius > 0 or xy_radius > 0:\n            struct_close = build_anisotropic_struct(z_radius, xy_radius)\n            if struct_close is not None:\n                mask = ndi.binary_closing(mask, structure=struct_close)\n        \n        if dust_min_size > 0:\n            mask = remove_small_objects(mask.astype(bool), min_size=dust_min_size)\n        \n        return mask.astype(np.uint8)\n    \n    else:\n        raise ValueError(f\"Unknown frangi_mode: {frangi_mode}\")\n\n# ==========================================\n# PREDICT (WITH ROTATION TTA)\n# ==========================================\n# def predict(\n#     sample,\n#     iid=None,\n#     T_low=0.50,\n#     T_high=0.90,\n#     z_radius=1,\n#     xy_radius=0,\n#     dust_min_size=100,\n# ):\n#     \"\"\"\n#     sample: (1, D, H, W, 1)\n#     \"\"\"\n\n#     # --------- ROTATION TTA PROBS ---------\n#     probs_fg = predict_probs_tta_rot(sample)   # (D, H, W)\n\n#     if iid is not None:\n#         np.save(iid, probs_fg)\n\n#     # --------- POSTPROCESS (UNCHANGED) ----\n#     final = topo_postprocess(\n#         probs_fg,\n#         T_low=T_low,\n#         T_high=T_high,\n#         z_radius=z_radius,\n#         xy_radius=xy_radius,\n#         dust_min_size=dust_min_size,\n#     )\n\n#     return final  # (D, H, W) uint8 {0,1}\n\ndef predict(\n    sample,  # TTA预测函数\n    iid=None,\n    T_low=0.50,\n    T_high=0.90,\n    z_radius=1,\n    xy_radius=0,\n    dust_min_size=100,\n    # Frangi参数\n    use_frangi=True,\n    frangi_weight=0.3,\n    frangi_sigmas=(1, 3, 5),\n    frangi_beta1=0.5,\n    frangi_beta2=15,\n    frangi_mode='enhance_after',\n):\n    \"\"\"\n    完整的预测函数，集成Frangi滤波器\n    \n    参数:\n        sample: 输入样本 (1, D, H, W, 1)\n        predict_probs_tta_rot_func: TTA预测函数\n        iid: 保存概率图的路径（可选）\n        其他参数: 后处理参数\n    \n    返回:\n        final: 最终掩码 (D, H, W) uint8 {0,1}\n    \"\"\"\n    # --------- ROTATION TTA PROBS ---------\n    probs_fg = predict_probs_tta_rot(sample)   # (D, H, W)\n\n    if iid is not None:\n        np.save(iid, probs_fg)\n\n    # --------- POSTPROCESS WITH FRANGI ----\n    final = topo_postprocess_with_frangi(\n        probs_fg,\n        T_low=T_low,\n        T_high=T_high,\n        z_radius=z_radius,\n        xy_radius=xy_radius,\n        dust_min_size=dust_min_size,\n        use_frangi=use_frangi,\n        frangi_weight=frangi_weight,\n        frangi_sigmas=frangi_sigmas,\n        frangi_beta1=frangi_beta1,\n        frangi_beta2=frangi_beta2,\n        frangi_mode=frangi_mode,\n    )\n\n    return final  # (D, H, W) uint8 {0,1}\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-27T04:09:18.776862Z","iopub.execute_input":"2026-01-27T04:09:18.777474Z","iopub.status.idle":"2026-01-27T04:09:18.798984Z","shell.execute_reply.started":"2026-01-27T04:09:18.777450Z","shell.execute_reply":"2026-01-27T04:09:18.798218Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"**Prediction and Zip Submission**","metadata":{}},{"cell_type":"code","source":"testing = False","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-27T04:09:18.799651Z","iopub.execute_input":"2026-01-27T04:09:18.799900Z","iopub.status.idle":"2026-01-27T04:09:18.817183Z","shell.execute_reply.started":"2026-01-27T04:09:18.799882Z","shell.execute_reply":"2026-01-27T04:09:18.816377Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"if testing:\n    \n    test_dir = \"/kaggle/input/vesuvius-challenge-surface-detection/train_images\"\n    test_df = pd.read_csv(f\"{root_dir}/train.csv\")\n    test_ids = {956073442, 961304774,969293709,975031774,985841575,992852942}\n    test_df = (\n        test_df\n        .loc[test_df[\"id\"].isin(test_ids)]\n        .reset_index(drop=True)\n    )","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-27T04:09:18.817838Z","iopub.execute_input":"2026-01-27T04:09:18.818025Z","iopub.status.idle":"2026-01-27T04:09:18.830940Z","shell.execute_reply.started":"2026-01-27T04:09:18.818009Z","shell.execute_reply":"2026-01-27T04:09:18.830314Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"with zipfile.ZipFile(\n    zip_path, \"w\", compression=zipfile.ZIP_DEFLATED\n) as z:\n    for image_id in test_df[\"id\"]:\n        tif_path = f\"{test_dir}/{image_id}.tif\"\n            \n        volume = load_volume(tif_path)\n        volume = val_transformation(volume)\n        if testing :\n            output = predict(volume,f\"{image_id}\") \n        else :\n            output = predict(volume)\n        \n        out_path = f\"{output_dir}/{image_id}.tif\"\n        tifffile.imwrite(out_path, output.astype(np.uint8))\n\n        z.write(out_path, arcname=f\"{image_id}.tif\")\n        os.remove(out_path)\n\nprint(\"Submission ZIP:\", zip_path)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-27T04:09:18.831809Z","iopub.execute_input":"2026-01-27T04:09:18.832651Z","iopub.status.idle":"2026-01-27T04:12:30.467563Z","shell.execute_reply.started":"2026-01-27T04:09:18.832632Z","shell.execute_reply":"2026-01-27T04:12:30.466756Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}