{"cells":[{"cell_type":"markdown","metadata":{},"source":"# Vesuvius Surface Detection - Winner Strategy v1\n\n**Objetivo: Subir del puesto 108 al TOP**\n\n### Técnicas implementadas (copiadas del TOP):\n1. **TransUNet + SEResNeXt101** - backbone potente\n2. **Rotation TTA 4x** (0°, 90°, 180°, 270°) - +0.02-0.03 score\n3. **Topology-aware Post-processing**:\n   - 3D Hysteresis thresholding\n   - Anisotropic morphological closing\n   - Connected components filtering\n4. **Overlap 0.52 Gaussian**\n\n### Score Formula:\n```\nScore = 0.30×TopoScore + 0.35×SurfaceDice@τ + 0.35×VOI\n```"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# Install packages\nimport os\nvar = \"/kaggle/input/vesuvius25-packages-offline-installer-v20251226/whls\"\nif os.path.exists(var):\n    !pip install -q \\\n        \"$var\"/keras_nightly-3.12.0.dev2025100703-py3-none-any.whl \\\n        \"$var\"/tifffile-2025.10.16-py3-none-any.whl \\\n        \"$var\"/imagecodecs-2025.11.11-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 --find-links \"$var\" 2>/dev/null\n    print(\"Packages installed\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"import os\nos.environ[\"KERAS_BACKEND\"] = \"jax\"\n\nimport numpy as np\nimport pandas as pd\nimport zipfile\nimport tifffile\nimport time\nimport scipy.ndimage as ndi\nfrom skimage.morphology import remove_small_objects\n\nimport keras\nfrom medicai.transforms import Compose, ScaleIntensityRange\nfrom medicai.models import TransUNet\nfrom medicai.utils.inference import SlidingWindowInference\n\nprint(f\"Keras backend: {keras.config.backend()}\")\nprint(f\"Keras version: {keras.version()}\")"},{"cell_type":"markdown","metadata":{},"source":"## Configuration"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"class CFG:\n    # Model\n    ENCODER = 'seresnext101'\n    INPUT_SHAPE = (160, 160, 160, 1)\n    NUM_CLASSES = 3\n    WEIGHTS = '/kaggle/input/colab-a-162v4-gpu-transunet-seresnext101-x160/model.weights.h5'\n\n    # Inference\n    OVERLAP = 0.52  # Gaussian overlap\n    SW_BATCH_SIZE = 1\n\n    # TTA - CRITICAL FOR WINNING\n    USE_ROTATION_TTA = True  # 4x rotation average (+0.02-0.03)\n\n    # Post-processing - TOPOLOGY AWARE (critical for TopoScore)\n    USE_TOPO_POSTPROC = True\n    T_LOW = 0.45          # Hysteresis low threshold\n    T_HIGH = 0.85         # Hysteresis high threshold\n    Z_RADIUS = 1          # Anisotropic closing Z\n    XY_RADIUS = 0         # Anisotropic closing XY\n    DUST_MIN_SIZE = 100   # Remove small components\n\n    # Paths\n    ROOT_DIR = \"/kaggle/input/vesuvius-challenge-surface-detection\"\n    TEST_DIR = f\"{ROOT_DIR}/test_images\"\n    OUTPUT_DIR = \"/kaggle/working/submission_masks\"\n    ZIP_PATH = \"/kaggle/working/submission.zip\"\n\nos.makedirs(CFG.OUTPUT_DIR, exist_ok=True)\nprint(\"Config loaded\")"},{"cell_type":"markdown","metadata":{},"source":"## Rotation TTA (Test-Time Augmentation)"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"def rot90_volume(vol, k):\n    \"\"\"\n    Rotate volume k times 90 degrees clockwise in HW plane.\n    vol: (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    \"\"\"Undo rotation\"\"\"\n    return rot90_volume(vol, (4 - k) % 4)\n\n\ndef predict_with_tta(predictor, sample):\n    \"\"\"\n    4x rotation TTA: 0, 90, 180, 270 degrees\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        # Rotate input\n        s_rot = rot90_volume(sample, k)\n\n        # Predict\n        out = predictor(s_rot)\n        out = np.asarray(out)\n        probs = out[0, ..., 1]  # (D, H, W) - class 1 (surface)\n\n        # Unrotate output\n        probs = unrot90_volume(probs, k)\n        probs_accum.append(probs)\n\n    # Average all rotations\n    return np.mean(probs_accum, axis=0)\n\nprint(\"TTA functions ready\")"},{"cell_type":"markdown","metadata":{},"source":"## Topology-Aware Post-Processing"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"def build_anisotropic_struct(z_radius: int, xy_radius: int):\n    \"\"\"Build anisotropic structuring element for 3D morphology\"\"\"\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\ndef topo_postprocess(probs, cfg=CFG):\n    \"\"\"\n    Topology-aware post-processing optimizado para TopoScore.\n    \n    Score formula: 0.30 x TopoScore + 0.35 x SurfaceDice + 0.35 x VOI\n    \n    Steps:\n    1. 3D Hysteresis thresholding (strong seeds propagated through weak regions)\n    2. 3D Anisotropic morphological closing (fill small gaps - improves topology)\n    3. Connected components filtering (remove dust - improves VOI)\n    \"\"\"\n    T_low = cfg.T_LOW\n    T_high = cfg.T_HIGH\n    z_radius = cfg.Z_RADIUS\n    xy_radius = cfg.XY_RADIUS\n    dust_min_size = cfg.DUST_MIN_SIZE\n\n    # Step 1: 3D Hysteresis Thresholding\n    # Strong regions seed the segmentation, weak regions can connect them\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)  # 26-connectivity\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 (fills small holes/gaps)\n    # Critical for TopoScore - fewer tunnels (beta_1)\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: Remove small connected components (dust)\n    # Critical for VOI score - fewer spurious components\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\nprint(\"Post-processing functions ready\")"},{"cell_type":"markdown","metadata":{},"source":"## Data Transform"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"def transform(image):\n    \"\"\"Normalize image to [0, 1]\"\"\"\n    data = {\"image\": image}\n    pipeline = Compose([\n        ScaleIntensityRange(\n            keys=[\"image\"],\n            a_min=0, a_max=255,\n            b_min=0, b_max=1,\n            clip=True\n        )\n    ])\n    return pipeline(data)[\"image\"]"},{"cell_type":"markdown","metadata":{},"source":"## Load Model"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# Load test data\ntest_df = pd.read_csv(f\"{CFG.ROOT_DIR}/test.csv\")\nprint(f\"Test samples: {len(test_df)}\")\n\n# Load model\nprint(f\"\\nLoading model: TransUNet + {CFG.ENCODER}\")\nmodel = TransUNet(\n    input_shape=CFG.INPUT_SHAPE,\n    encoder_name=CFG.ENCODER,\n    classifier_activation='softmax',\n    num_classes=CFG.NUM_CLASSES,\n)\nmodel.load_weights(CFG.WEIGHTS)\nprint(f\"Parameters: {model.count_params()/1e6:.1f}M\")\n\n# Sliding Window Inference\nswi = SlidingWindowInference(\n    model,\n    num_classes=CFG.NUM_CLASSES,\n    roi_size=CFG.INPUT_SHAPE[:3],\n    sw_batch_size=CFG.SW_BATCH_SIZE,\n    mode='gaussian',\n    overlap=CFG.OVERLAP\n)\nprint(\"Model loaded\")"},{"cell_type":"markdown","metadata":{},"source":"## Inference"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# Config summary\nprint(\"=\" * 60)\nprint(\"WINNER STRATEGY v1 - Subiendo del puesto 108 al TOP\")\nprint(\"=\" * 60)\nprint(f\"Model: TransUNet + {CFG.ENCODER}\")\nprint(f\"Overlap: {CFG.OVERLAP}\")\nprint(f\"Rotation TTA: {CFG.USE_ROTATION_TTA} (4x boost)\")\nprint(f\"Topo Post-proc: {CFG.USE_TOPO_POSTPROC}\")\nif CFG.USE_TOPO_POSTPROC:\n    print(f\"  - Hysteresis: T_low={CFG.T_LOW}, T_high={CFG.T_HIGH}\")\n    print(f\"  - Closing: z={CFG.Z_RADIUS}, xy={CFG.XY_RADIUS}\")\n    print(f\"  - Dust removal: min_size={CFG.DUST_MIN_SIZE}\")\nprint(\"=\" * 60)"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# Main inference loop\nwith zipfile.ZipFile(CFG.ZIP_PATH, \"w\", zipfile.ZIP_DEFLATED) as z:\n    for idx, row in enumerate(test_df.itertuples()):\n        image_id = row.id\n        print(f\"[{idx+1}/{len(test_df)}] {image_id}\")\n        t0 = time.time()\n\n        # Load volume\n        vol = tifffile.imread(f\"{CFG.TEST_DIR}/{image_id}.tif\").astype(np.float32)\n        vol = vol[None, ..., None]  # (1, D, H, W, 1)\n        vol = transform(vol)\n\n        # Inference with/without TTA\n        if CFG.USE_ROTATION_TTA:\n            probs = predict_with_tta(swi, vol)\n        else:\n            out = swi(vol)\n            probs = np.asarray(out)[0, ..., 1]  # (D, H, W)\n\n        # Post-processing\n        if CFG.USE_TOPO_POSTPROC:\n            output = topo_postprocess(probs)\n        else:\n            output = (probs > 0.5).astype(np.uint8)\n\n        # Stats\n        surface_voxels = output.sum()\n        elapsed = time.time() - t0\n        print(f\"  Surface voxels: {surface_voxels:,} | Time: {elapsed:.1f}s\")\n\n        # Save\n        out_path = f\"{CFG.OUTPUT_DIR}/{image_id}.tif\"\n        tifffile.imwrite(out_path, output)\n        z.write(out_path, arcname=f\"{image_id}.tif\")\n        os.remove(out_path)\n\nprint(f\"\\nSubmission ready: {CFG.ZIP_PATH}\")"},{"cell_type":"markdown","metadata":{},"source":"## Summary\n\n**Mejoras implementadas vs tu version anterior:**\n\n| Tecnica | Antes | Ahora | Impacto estimado |\n|---------|-------|-------|------------------|\n| Rotation TTA | No | 4x | +0.02-0.03 |\n| Hysteresis threshold | No | Si | +0.01-0.02 |\n| Morphological closing | No | Si | +0.005-0.01 |\n| Dust removal | No | Si | +0.005-0.01 |\n\n**Esperado: Subir de ~0.46 a ~0.52-0.55**"}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.10.0"}},"nbformat":4,"nbformat_minor":4}