{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.11.0"},"kaggle":{"accelerator":"nvidiaTeslaT4","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":732880,"sourceType":"modelInstanceVersion","modelInstanceId":516822,"modelId":510647}],"dockerImageVersionId":31193,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"
\n

\n Vesuvius Surface Detection\n

\n
","metadata":{}},{"cell_type":"markdown","source":"
\n3D Segmentation of Ancient Papyrus Scroll Surfaces from CT Scans
\nTransUNet + TTA + Post-Processing Pipeline\n
","metadata":{}},{"cell_type":"markdown","source":"## Task Overview\n\nThe **Vesuvius Challenge - Surface Detection** competition asks participants to identify the **recto (front) surface** of ancient papyrus scrolls from 3D CT scan volumes. These scrolls, carbonized by the eruption of Mount Vesuvius in 79 AD, are too fragile to physically unroll. By accurately detecting the scroll surface in 3D, we enable downstream ink detection algorithms to read the ancient text.\n\n### Pipeline Architecture\n\n```\nInput (3D TIFF) → Normalize → TransUNet → Sliding Window → TTA (x8) → Hysteresis → Closing → Output\n```\n\n| Component | Description | Configuration |\n|-----------|-------------|---------------|\n| **Model** | TransUNet with SEResNeXt50 encoder + Vision Transformer bottleneck | 70M parameters |\n| **Input** | 3D CT volume (160×160×160 voxels) | Float32, z-score normalized |\n| **Output** | 3 classes: Background (0), Surface (1), Unlabeled (2) | Logits → Argmax |\n| **Sliding Window** | Gaussian-weighted inference with overlap | overlap=0.46 |\n| **TTA** | 3 flips + 3 rotations = 8 views averaged | Logit-space averaging |\n| **Hysteresis** | Two-threshold segmentation | T_low=0.30, T_high=0.80 |\n| **Morphological Closing** | Anisotropic structuring element | z_radius=3, xy_radius=2 |\n| **Dust Removal** | Remove small connected components | min_size=100 voxels |\n\n### Key Design Choices\n\n| Choice | Rationale |\n|--------|----------|\n| **Single model** (no ensemble) | Comboloss model alone outperforms multi-model ensembles |\n| **classifier_activation=None** | Logit averaging preserves sharper class boundaries than softmax |\n| **overlap=0.46** | Tuned value; balances boundary consistency vs computation |\n| **Aggressive closing (z=3)** | Fills gaps in thin 3-voxel sheet predictions |","metadata":{}},{"cell_type":"markdown","source":"#
Import Libraries
","metadata":{}},{"cell_type":"code","source":"from IPython.display import clear_output\n\nvar=\"/kaggle/input/vsdetection-packages-offline-installer-only/whls\"\n!pip install \\\n \"$var\"/keras_nightly-*.whl \\\n \"$var\"/tifffile-*.whl \\\n \"$var\"/imagecodecs-*.whl \\\n \"$var\"/medicai-*.whl \\\n --no-index \\\n --find-links \"$var\"\n\nclear_output()","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nos.environ[\"KERAS_BACKEND\"] = \"jax\"\n\nimport keras\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()","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"#
Load 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":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"test_df = pd.read_csv(f\"{root_dir}/test.csv\")\nprint(f\"Test samples: {len(test_df)}\")\ntest_df.head()","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"#
Data Transformation
","metadata":{}},{"cell_type":"code","source":"def val_transformation(image):\n \"\"\"Z-score normalization on nonzero voxels.\"\"\"\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\"]","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"#
Model Configuration
","metadata":{}},{"cell_type":"code","source":"num_classes = 3\ninput_shape = (160, 160, 160)\nkaggle_model_path = \"/kaggle/input/vsd-model/keras/\"","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def get_model():\n \"\"\"Load TransUNet with SEResNeXt50 encoder (comboloss weights).\"\"\"\n model = TransUNet(\n input_shape=(160, 160, 160, 1),\n encoder_name='seresnext50',\n classifier_activation=None,\n num_classes=3,\n )\n model.load_weights(\n f\"{kaggle_model_path}/transunet/3/transunet.seresnext50.160px.comboloss.weights.h5\"\n )\n return model","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"model = get_model()\nprint(f\"Model parameters: {model.count_params() / 1e6:.1f}M\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"model.instance_describe()","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"#
Sliding Window Inference
","metadata":{}},{"cell_type":"code","source":"swi = SlidingWindowInference(\n model,\n num_classes=3,\n roi_size=input_shape,\n sw_batch_size=1,\n mode='gaussian',\n overlap=0.46,\n)","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def load_volume(path):\n \"\"\"Load 3D TIFF and add batch/channel dimensions.\"\"\"\n vol = tifffile.imread(path)\n vol = vol.astype(np.float32)\n vol = vol[None, ..., None]\n return vol","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"#
Test Time Augmentation
","metadata":{}},{"cell_type":"code","source":"def predict_with_tta(inputs, swi):\n \"\"\"TTA with 3 flips + 3 rotations = 8 views, averaged in logit space.\"\"\"\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.argmax(-1).astype(np.uint8).squeeze()","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"#
Post Processing
","metadata":{}},{"cell_type":"code","source":"def build_anisotropic_struct(z_radius: int, xy_radius: int):\n \"\"\"Build 3D structuring element with different radii for z and xy axes.\"\"\"\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 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):\n \"\"\"3-stage post-processing: Hysteresis → Closing → Dust removal.\"\"\"\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 return mask.astype(np.uint8)","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"#
Inference Pipeline
","metadata":{}},{"cell_type":"code","source":"def inference_pipelines(\n volume,\n T_low=0.30,\n T_high=0.80,\n z_radius=3,\n xy_radius=2,\n dust_min_size=100,\n):\n \"\"\"Complete inference: TTA prediction → Post-processing.\"\"\"\n probs = predict_with_tta(volume, swi)\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 )\n return final","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"#
Generate Submission
","metadata":{}},{"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 output = inference_pipelines(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":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"#
Visualization
","metadata":{}},{"cell_type":"code","source":"def plot_sample(x, y, sample_idx=0, max_slices=16):\n \"\"\"Visualize input volume and predicted mask slices.\"\"\"\n img = np.squeeze(x[sample_idx])\n mask = np.squeeze(y[sample_idx])\n D = img.shape[0]\n\n step = max(1, D // max_slices)\n slices = range(0, D, step)\n\n n_slices = len(slices)\n fig, axes = plt.subplots(2, n_slices, figsize=(3*n_slices, 6))\n\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\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\n plt.suptitle(f\"Sample {sample_idx}\")\n plt.tight_layout()\n plt.show()","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plot_sample(\n volume.numpy(), output[None], sample_idx=0, max_slices=5\n)","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"
\nThank you for reading! If you find this notebook useful, please consider giving it an upvote.\n
","metadata":{}}]}