{"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":[{"sourceType":"competition","sourceId":117682,"databundleVersionId":15062069},{"sourceType":"modelInstanceVersion","sourceId":732880,"databundleVersionId":15477237,"modelInstanceId":516822,"modelId":510647},{"sourceType":"kernelVersion","sourceId":290917305}],"dockerImageVersionId":31193,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## Inference From Tony Li · PaulG · Yiheng Wang\nhttps://www.kaggle.com/code/tonylica/vesuvius-0-552?scriptVersionId=295481156\n\n## Training Notebooks From Innat\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","metadata":{}},{"cell_type":"code","source":"from IPython.display import clear_output\nimport os\n\n# 环境配置\nos.environ.setdefault(\"PROTOCOL_BUFFERS_PYTHON_IMPLEMENTATION\", \"python\")\nos.environ[\"KERAS_BACKEND\"] = \"jax\"\nos.environ.setdefault(\"OMP_NUM_THREADS\", \"4\")\nos.environ.setdefault(\"MKL_NUM_THREADS\", \"4\")\nos.environ.setdefault(\"OPENBLAS_NUM_THREADS\", \"4\")\nos.environ.setdefault(\"NUMEXPR_NUM_THREADS\", \"4\")\n\n# 安装依赖\nvar=\"/kaggle/input/vsdetection-packages-offline-installer-only/whls\"\n!pip install \"$var\"/keras_nightly-*.whl \"$var\"/tifffile-*.whl \"$var\"/imagecodecs-*.whl \"$var\"/medicai-*.whl --no-index --find-links \"$var\"\nclear_output()\n\nimport time, zipfile\nimport numpy as np\nimport pandas as pd\nimport tifffile\nimport scipy.ndimage as ndi\nfrom skimage.morphology import remove_small_objects\nimport keras\nfrom medicai.transforms import Compose, NormalizeIntensity\nfrom medicai.models import TransUNet\nfrom medicai.utils.inference import SlidingWindowInference\nimport keras\nfrom keras import ops\nfrom medicai.transforms import (\n    Compose,\n    ScaleIntensityRange,\n    NormalizeIntensity\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\nimport scipy.ndimage as ndi\nfrom skimage.morphology import remove_small_objects\nfrom matplotlib import pyplot as plt\n\nkeras.config.backend(), keras.version()\n\nroot_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)\ntest_df = pd.read_csv(f\"{root_dir}/test.csv\")\ntest_df.head()\n\ndef val_transformation(image):\n    data = {\"image\": image}\n    pipeline = Compose([\n        NormalizeIntensity(\n            keys=[\"image\"], \n            nonzero=True,\n            channel_wise=False\n        ),\n    ])\n    result = pipeline(data)\n    return result[\"image\"]\n\ntta=1\nnum_classes=3\ninput_shape=(160, 160, 160)\ndef get_model_a():\n    # 0.551 (tta+pp)\n    model = TransUNet(\n        input_shape=(160, 160, 160, 1),\n        encoder_name='seresnext50',\n        classifier_activation=None, # For tta, do softmax later.\n        num_classes=3,\n    )\n    model.load_weights(\n        \"/kaggle/input/models/ipythonx/vsd-model/keras/transunet/4/transunet.seresnext50.160px.duel.weights.h5\"\n    )\n    return model\n\ndef get_model_b():\n    # 0.522 (tta+pp)\n    model = SegFormer(\n        input_shape=input_shape + (1,),\n        encoder_name='mit_b4',\n        classifier_activation=None, # For tta, do softmax later.\n        num_classes=num_classes,\n    )\n    model.load_weights(\n        \"/kaggle/input/models/ipythonx/vsd-model/keras/segformer/1/segformer.mit.b4.weights.h5\"\n    )\n    return model\nmodels = [get_model_a(), get_model_b()]\n[model.count_params() / 1e6 for model in models]\n\nfor model in models:\n    print(model.instance_describe())\n\nswis = [\n    SlidingWindowInference(\n        model=m,\n        num_classes=num_classes,\n        roi_size=input_shape,\n        sw_batch_size=1,\n        overlap=0.5,\n        mode=\"gaussian\",\n    )\n    for m in models\n]\ndef load_volume(path):\n    vol = tifffile.imread(path)\n    vol = vol.astype(np.float32)\n    vol = vol[None, ..., None]\n    return vol\n\ndef predict_with_tta(inputs, swi):\n    logits = []\n\n    # Original\n    logits.append(swi(inputs))\n\n    # Flips (spatial only)\n    for axis in [1, 2, 3]:\n        img_f = np.flip(inputs, axis=axis)\n        p = swi(img_f)\n        p = np.flip(p, axis=axis)\n        logits.append(p)\n\n    # Axial rotations (H, W)\n    for k in [1, 2, 3]:\n        img_r = np.rot90(inputs, k=k, axes=(2, 3))\n        p = swi(img_r)\n        p = np.rot90(p, k=-k, axes=(2, 3))\n        logits.append(p)\n\n    mean_logits = np.mean(logits, axis=0)\n    return mean_logits\n# https://www.kaggle.com/code/choudharymanas/inference-baseline-transunet-lb-0-537\ndef build_anisotropic_struct(z_radius: int, xy_radius: int):\n    z, r = z_radius, xy_radius\n    if z == 0 and r == 0:\n        return None\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    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    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# 边缘掩码函数\ndef mask_edge_artifacts(mask, margin=10):\n    if margin <= 0:\n        return mask\n    d, h, w = mask.shape\n    mask[:margin, :, :] = 0\n    mask[-margin:, :, :] = 0\n    mask[:, :margin, :] = 0\n    mask[:, -margin:, :] = 0\n    mask[:, :, :margin] = 0\n    mask[:, :, -margin:] = 0\n    return mask\n\ndef topo_postprocess(\n    probs,\n    T_low=0.90,\n    T_high=0.90,\n    z_radius=1,\n    xy_radius=0,\n    dust_min_size=100,\n    edge_margin=10,\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(\n        strong, mask=weak, structure=struct_hyst\n    )\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(\n            mask.astype(bool), min_size=dust_min_size\n        )\n\n    # Step 4: Edge Masking\n    mask = mask_edge_artifacts(mask, margin=edge_margin)\n\n    return mask.astype(np.uint8)\n\nweights = np.array([0.8, 0.2], dtype=np.float32)\ndef predict_weighted_ensemble(volume):\n    logits = []\n    for swi in swis:\n        logits.append(\n            predict_with_tta(volume, swi) if tta else swi(volume) \n        )\n    logits = np.stack(logits, axis=0)\n    weighted_logits = logits * weights[:, None, None, None, None, None]\n    mean_logits = weighted_logits.sum(axis=0) / weights.sum()\n    probs = ops.softmax(mean_logits, axis=-1)\n    return probs.argmax(axis=-1).astype(np.uint8).squeeze()\n\ndef inference_pipelines(\n    volume,\n    T_low=0.50,\n    T_high=0.90,\n    z_radius=1,\n    xy_radius=0,\n    dust_min_size=100,\n    edge_margin=10,\n):\n    probs = predict_weighted_ensemble(volume)\n    final = topo_postprocess(\n        probs,\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        edge_margin=edge_margin,\n    )\n    return final\n\n# 主执行流程\nprint(f\"开始预测，共{len(test_df['id'])}个测试样本\")\nt_global0 = time.perf_counter()\n\nwith zipfile.ZipFile(\n    zip_path, \"w\", compression=zipfile.ZIP_DEFLATED\n) as z:\n    for i, image_id in enumerate(test_df[\"id\"]):\n        t0 = time.perf_counter()\n        tif_path = f\"{test_dir}/{image_id}.tif\"\n\n        volume = load_volume(tif_path)\n        volume = val_transformation(volume)\n        output = inference_pipelines(\n            volume, T_low=0.35, T_high=0.85, z_radius=3, xy_radius=2, dust_min_size=200, edge_margin=5\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\n        dt = time.perf_counter() - t0\n        elapsed = time.perf_counter() - t_global0\n        print(f\"[{i+1}/{len(test_df['id'])}] {image_id} | 耗时{dt/60:.2f}分钟 | 累计{elapsed/3600:.2f}小时 | 阳性像素数={int(output.sum())}\")\n\nprint(\"Submission ZIP:\", zip_path)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-22T05:37:06.413574Z","iopub.execute_input":"2026-02-22T05:37:06.41412Z","iopub.status.idle":"2026-02-22T05:42:36.506816Z","shell.execute_reply.started":"2026-02-22T05:37:06.414096Z","shell.execute_reply":"2026-02-22T05:42:36.506075Z"},"_kg_hide-output":true},"outputs":[],"execution_count":null}]}