{"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,"isSourceIdPinned":false},{"sourceType":"modelInstanceVersion","sourceId":732880,"databundleVersionId":15477237,"modelInstanceId":516822,"modelId":510647,"isSourceIdPinned":false},{"sourceType":"kernelVersion","sourceId":290917305,"isSourceIdPinned":false},{"sourceType":"kernelVersion","sourceId":299209845,"isSourceIdPinned":false}],"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\n\n# 配置参数\nCFG = dict(\n    kaggle_model_path=\"/kaggle/input/notebooks/tonyai007/\",\n    weights_relpath=\"train-vesuvius-seresnext50-comboloss-5/final_model.weights.h5\", # 830轮\n    overlap_public=0.42,\n    overlap_base=0.48,\n    overlap_hi=0.6,\n    OV06_MAIN_ONLY=True,\n    USE_TTA=True,\n    INK_MODE=\"fg12\",\n    T_low=0.35,\n    T_high=0.85,\n    z_radius=0,\n    xy_radius=2,\n    dust_min_size=250,\n    DO_WARMUP=True,\n)\n\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)\n\n# 加载测试数据\ntest_df = pd.read_csv(f\"{root_dir}/test.csv\")\nids = test_df[\"id\"].tolist()\nROI = (160, 160, 160)\n\n# 数据预处理\n_val_pipeline = Compose([NormalizeIntensity(keys=[\"image\"], nonzero=True, channel_wise=False)])\ndef val_transformation(image):\n    return _val_pipeline({\"image\": image})[\"image\"]\n\ndef load_volume(path):\n    vol = tifffile.imread(path).astype(np.float32)\n    return vol[None, ..., None]\n\n# 数值计算函数\ndef sigmoid_stable(x):\n    x = np.asarray(x, dtype=np.float32)\n    out = np.empty_like(x, dtype=np.float32)\n    pos = x >= 0\n    out[pos] = 1.0 / (1.0 + np.exp(-x[pos]))\n    ex = np.exp(x[~pos])\n    out[~pos] = ex / (1.0 + ex)\n    return out\n\ndef logsumexp2(a, b):\n    a = np.asarray(a, dtype=np.float32)\n    b = np.asarray(b, dtype=np.float32)\n    m = np.maximum(a, b)\n    return m + np.log(np.exp(a - m) + np.exp(b - m) + 1e-12)\n\ndef binary_logit_from_multiclass_logits(logits_5d, mode=\"fg12\"):\n    x = np.asarray(logits_5d, dtype=np.float32)[0]\n    L0, L1, L2 = x[...,0], x[...,1], x[...,2]\n    if mode == \"fg12\":\n        return (logsumexp2(L1, L2) - L0).astype(np.float32, copy=False)\n    elif mode == \"class1\":\n        return (L1 - logsumexp2(L0, L2)).astype(np.float32, copy=False)\n    else:\n        raise ValueError(f\"Unknown INK_MODE={mode}\")\n\ndef build_anisotropic_struct(z_radius: int, xy_radius: int):\n    z, r = int(z_radius), int(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\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 = z; cy = cx = 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 seeded_hysteresis_with_topology(prob, pub_fg_bool, T_low=0.35, T_high=0.85, z_radius=0, xy_radius=2, dust_min_size=250):\n    prob = np.asarray(prob, dtype=np.float32)\n    strong = prob >= float(T_high)\n    if not strong.any():\n        return np.zeros_like(prob, dtype=np.uint8)\n    \n    weak = (prob >= float(T_low)) | pub_fg_bool\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(prob, dtype=np.uint8)\n    \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    mask = remove_small_objects(mask.astype(bool), min_size=int(dust_min_size))\n    return mask.astype(np.uint8)\n\n# 加载模型和滑动窗口推理\nweights_path = f\"{CFG['kaggle_model_path']}/{CFG['weights_relpath']}\"\nmodel = TransUNet(input_shape=(160, 160, 160, 1), encoder_name=\"seresnext50\", classifier_activation=None, num_classes=3)\nmodel.load_weights(weights_path)\n\ndef build_swi(overlap):\n    return SlidingWindowInference(model, num_classes=3, roi_size=ROI, sw_batch_size=1, mode=\"gaussian\", overlap=float(overlap))\n\nswi_public = build_swi(CFG[\"overlap_public\"])\nswi_base   = build_swi(CFG[\"overlap_base\"])\nswi_hi     = build_swi(CFG[\"overlap_hi\"])\n\n# TTA数据增强\ndef iter_tta(volume):\n    yield volume, (lambda y: y)\n    for axis in [1, 2, 3]:\n        v = np.flip(volume, axis=axis)\n        inv = (lambda y, axis=axis: np.flip(y, axis=axis))\n        yield v, inv\n    for k in [1, 2, 3]:\n        v = np.rot90(volume, k=k, axes=(2, 3))\n        inv = (lambda y, k=k: np.rot90(y, k=-k, axes=(2, 3)))\n        yield v, inv\n\n# 预测函数\ndef predict_pub_labels_and_private_prob(volume):\n    mode = CFG[\"INK_MODE\"]\n    if not CFG[\"USE_TTA\"]:\n        l_pub = np.asarray(swi_public(volume))\n        pub_labels = l_pub.argmax(-1).astype(np.uint8).squeeze()\n        l_prv = np.asarray(swi_hi(volume))\n        s = binary_logit_from_multiclass_logits(l_prv, mode=mode)\n        prob = sigmoid_stable(s)\n        return pub_labels, prob\n\n    logits_sum = None\n    s_sum = None\n    n = 0\n    for t, (v, inv) in enumerate(iter_tta(volume)):\n        l_pub = np.asarray(swi_public(v))\n        l_pub = inv(l_pub)\n        logits_sum = l_pub.astype(np.float32) if logits_sum is None else (logits_sum + l_pub.astype(np.float32))\n        \n        swi_use = swi_hi if (t == 0) else swi_base if CFG[\"OV06_MAIN_ONLY\"] else swi_hi\n        l_prv = np.asarray(swi_use(v))\n        l_prv = inv(l_prv)\n        s = binary_logit_from_multiclass_logits(l_prv, mode=mode)\n        s_sum = s.astype(np.float32) if s_sum is None else (s_sum + s.astype(np.float32))\n        n += 1\n\n    mean_logits = logits_sum / float(n)\n    pub_labels = mean_logits.argmax(-1).astype(np.uint8).squeeze()\n    s_mean = (s_sum / float(n)).astype(np.float32, copy=False)\n    prob = sigmoid_stable(s_mean)\n    return pub_labels, prob\n\n# 预热函数\ndef warmup(volume):\n    _ = np.asarray(swi_public(volume))\n    _ = np.asarray(swi_base(volume))\n    _ = np.asarray(swi_hi(volume))\n\n# 主执行流程\nprint(f\"开始预测，共{len(ids)}个测试样本\")\nt_global0 = time.perf_counter()\n\nwith zipfile.ZipFile(zip_path, \"w\", compression=zipfile.ZIP_DEFLATED) as zf:\n    for i, image_id in enumerate(ids):\n        t0 = time.perf_counter()\n        tif_path = f\"{test_dir}/{image_id}.tif\"\n        volume = load_volume(tif_path)\n        volume = val_transformation(volume)\n\n        if i == 0 and CFG[\"DO_WARMUP\"]:\n            print(\"JAX预热中...\")\n            warmup(volume)\n\n        pub_labels, prob = predict_pub_labels_and_private_prob(volume)\n        pub_fg = (pub_labels != 0)\n        output = seeded_hysteresis_with_topology(prob, pub_fg_bool=pub_fg, T_low=CFG[\"T_low\"], T_high=CFG[\"T_high\"], z_radius=CFG[\"z_radius\"], xy_radius=CFG[\"xy_radius\"], dust_min_size=CFG[\"dust_min_size\"])\n\n        out_path = f\"{output_dir}/{image_id}.tif\"\n        tifffile.imwrite(out_path, output.astype(np.uint8))\n        zf.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(ids)}] {image_id} | 耗时{dt/60:.2f}分钟 | 累计{elapsed/3600:.2f}小时 | 阳性像素数={int(output.sum())}\")\n\nprint(f\"生成提交文件: {zip_path}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-27T03:48:49.299236Z","iopub.execute_input":"2026-02-27T03:48:49.299914Z","iopub.status.idle":"2026-02-27T03:54:23.461622Z","shell.execute_reply.started":"2026-02-27T03:48:49.299888Z","shell.execute_reply":"2026-02-27T03:54:23.460652Z"},"_kg_hide-output":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import matplotlib.pyplot as plt\ndef plot_sample(x, y, sample_idx=0, max_slices=16):\n    img = np.squeeze(x[sample_idx])\n    mask = np.squeeze(y[sample_idx])\n    D = img.shape[0]\n    step = max(1, D // max_slices)\n    slices = range(0, D, step)\n    n_slices = len(slices)\n    fig, axes = plt.subplots(2, n_slices, figsize=(3*n_slices, 6))\n    for i, s in enumerate(slices):\n        axes[0, i].imshow(img[s], cmap='gray')\n        axes[0, i].set_title(f\"Slice {s}\")\n        axes[0, i].axis('off')\n        axes[1, i].imshow(mask[s], cmap='gray')\n        axes[1, i].set_title(f\"Mask {s}\")\n        axes[1, i].axis('off')\n    plt.suptitle(f\"Sample {sample_idx}\")\n    plt.tight_layout()\n    plt.show()\nplot_sample(volume.numpy(), output[None], sample_idx=0, max_slices=5)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-27T03:54:23.462923Z","iopub.execute_input":"2026-02-27T03:54:23.463159Z","iopub.status.idle":"2026-02-27T03:54:24.36917Z","shell.execute_reply.started":"2026-02-27T03:54:23.463138Z","shell.execute_reply":"2026-02-27T03:54:24.368318Z"}},"outputs":[],"execution_count":null}]}