{"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":"nvidiaTeslaT4","dataSources":[{"sourceId":117682,"databundleVersionId":14443416,"sourceType":"competition"}],"dockerImageVersionId":31192,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Let's Read Ancient Papers with Cold Coffee: The Simplest Way","metadata":{}},{"cell_type":"markdown","source":"## The Ancient Coffee Mystery ☕🏺\n\n**Picture this:** You're a barista-archaeologist staring at a coffee scroll (papyrus) that's been carbonized into volcanic ash toast for 2,000 years. \n☠️ It's so fragile that touching it = instant coffee dust. \n\nBut inside? Priceless Roman coffee recipes!\n\n**Your Mission:** Use CT scans (X-ray vision for coffee!) to find the papyrus layers without physically unrolling. You're mapping the coffee surface through folds, gaps, and ash damage, like finding a single coffee bean in a crumpled napkin! 🎯","metadata":{}},{"cell_type":"markdown","source":"**Why This is HARD :**\n\nSuper fragile - One wrong move = ☠️\n\nBurnt to a crisp - Looks like ash, hard to distinguish\n\nCrazy folded - Like origami from hell\n\nTopology matters - Can't glue layers together or split them apart!\n\nNoise everywhere - Ash creates fake coffee artifacts","metadata":{}},{"cell_type":"markdown","source":"## 1: ☕  ☕ Let's Gather Our Coffee Ingredients (Imports & Setup)\nLoad all necessary libraries and define global parameters used throughout the pipeline","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport torch\nimport torch.nn as nn\nfrom pathlib import Path\nfrom tqdm.auto import tqdm\nfrom PIL import Image, ImageSequence\nfrom io import BytesIO\nimport zipfile\nfrom scipy import ndimage\n\n\n# Our coffee recipe parameters\n\nTOP_PERCENT = 25.0          # How much coffee foam to keep (top 25% brightest)\nMIN_COMPONENT_VOXELS = 400  # Remove tiny coffee grounds (noise)\nCLOSING_ITERS = 1           # Gentle stirring (gap filling)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-20T10:03:45.030243Z","iopub.execute_input":"2025-11-20T10:03:45.030525Z","iopub.status.idle":"2025-11-20T10:03:45.035370Z","shell.execute_reply.started":"2025-11-20T10:03:45.030505Z","shell.execute_reply":"2025-11-20T10:03:45.034588Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 2: 🏗️ Building Our Espresso Machine (Neural Network)\nEvery great cold coffee needs a powerful espresso machine! This is our 3D U-Net - it's like a super smart coffee filter that knows exactly which beans (voxels) are the good stuff :","metadata":{}},{"cell_type":"code","source":"class DoubleConv3D(nn.Module):\n    def __init__(self, in_ch, out_ch):\n        super().__init__()\n        self.conv = nn.Sequential(\n            nn.Conv3d(in_ch, out_ch, 3, padding=1),\n            nn.BatchNorm3d(out_ch),\n            nn.ReLU(inplace=True),\n            nn.Conv3d(out_ch, out_ch, 3, padding=1),\n            nn.BatchNorm3d(out_ch),\n            nn.ReLU(inplace=True)\n        )\n    def forward(self, x):\n        return self.conv(x)\n\nclass UNet3D(nn.Module):\n    def __init__(self, in_channels=1, out_channels=1, init_features=24):\n        super().__init__()\n        f = init_features\n        self.enc1 = DoubleConv3D(in_channels, f)\n        self.pool1 = nn.MaxPool3d(2)\n        self.enc2 = DoubleConv3D(f, f*2)\n        self.pool2 = nn.MaxPool3d(2)\n        self.enc3 = DoubleConv3D(f*2, f*4)\n        self.pool3 = nn.MaxPool3d(2)\n        self.bottleneck = DoubleConv3D(f*4, f*8)\n        self.upconv3 = nn.ConvTranspose3d(f*8, f*4, 2, stride=2)\n        self.dec3 = DoubleConv3D(f*8, f*4)\n        self.upconv2 = nn.ConvTranspose3d(f*4, f*2, 2, stride=2)\n        self.dec2 = DoubleConv3D(f*4, f*2)\n        self.upconv1 = nn.ConvTranspose3d(f*2, f, 2, stride=2)\n        self.dec1 = DoubleConv3D(f*2, f)\n        self.out = nn.Conv3d(f, out_channels, 1)\n    \n    def forward(self, x):\n        e1 = self.enc1(x)\n        e2 = self.enc2(self.pool1(e1))\n        e3 = self.enc3(self.pool2(e2))\n        b = self.bottleneck(self.pool3(e3))\n        d3 = self.upconv3(b)\n        d3 = torch.cat([d3, e3], 1)\n        d3 = self.dec3(d3)\n        d2 = self.upconv2(d3)\n        d2 = torch.cat([d2, e2], 1)\n        d2 = self.dec2(d2)\n        d1 = self.upconv1(d2)\n        d1 = torch.cat([d1, e1], 1)\n        d1 = self.dec1(d1)\n        return torch.sigmoid(self.out(d1))\n\ndef load_tiff(path):\n    with Image.open(path) as tif:\n        return np.stack([np.array(frame) for frame in ImageSequence.Iterator(tif)])\n\ndef save_to_zip(volume, zip_file, filename):\n    volume = (np.clip(volume, 0, 1) * 255).astype(np.uint8)\n    pages = [Image.fromarray(v) for v in volume]\n    buffer = BytesIO()\n    pages[0].save(buffer, format='TIFF', save_all=True, append_images=pages[1:])\n    zip_file.writestr(filename, buffer.getvalue())\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-20T10:04:32.111287Z","iopub.execute_input":"2025-11-20T10:04:32.111580Z","iopub.status.idle":"2025-11-20T10:04:32.122989Z","shell.execute_reply.started":"2025-11-20T10:04:32.111560Z","shell.execute_reply":"2025-11-20T10:04:32.122095Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"##  3: ⚡ Our Coffee Grinding & Brewing Tools (Core Functions)\n\nTime to grind those coffee beans! These functions are our coffee grinder, blender, and strainer, turning raw CT scans into perfect predictions!","metadata":{}},{"cell_type":"code","source":"def dl_inference_fast(volume, model, device, stride=72):\n    vol = (volume - volume.mean()) / (volume.std() + 1e-8)\n    Z,H,W = vol.shape\n    ps = 96\n    out = np.zeros((Z,H,W), dtype=np.float32)\n    cnt = np.zeros((Z,H,W), dtype=np.float32)\n    \n    with torch.no_grad():\n        for z in range(0, max(1,Z-ps+1), stride):\n            for y in range(0, max(1,H-ps+1), stride):\n                for x in range(0, max(1,W-ps+1), stride):\n                    ze,ye,xe = min(z+ps,Z), min(y+ps,H), min(x+ps,W)\n                    patch = vol[z:ze,y:ye,x:xe]\n                    if patch.shape != (ps,ps,ps):\n                        pad = np.zeros((ps,ps,ps))\n                        pad[:patch.shape[0],:patch.shape[1],:patch.shape[2]] = patch\n                        patch = pad\n                    pred = model(torch.from_numpy(patch[None,None]).float().to(device))[0,0].cpu().numpy()\n                    out[z:ze,y:ye,x:xe] += pred[:ze-z,:ye-y,:xe-x]\n                    cnt[z:ze,y:ye,x:xe] += 1\n    \n    return out / np.maximum(cnt, 1)\n\ndef classical_cv_slice_adaptive(volume):\n    filtered = ndimage.median_filter(volume, size=(1, 3, 3))\n    masks = []\n    percentile_cutoff = 100.0 - TOP_PERCENT\n    for sl in filtered:\n        t = np.percentile(sl, percentile_cutoff)\n        masks.append((sl > t).astype(np.float32))\n    return np.stack(masks)\n\ndef slice_adaptive_ensemble(classical, dl):\n    ensemble = np.zeros_like(classical)\n    for i in range(classical.shape[0]):\n        slice_classical = classical[i]\n        slice_dl = dl[i]\n        weight = 0.7 + 0.2 * np.abs(slice_classical - 0.5)\n        ensemble[i] = weight * slice_classical + (1 - weight) * slice_dl\n    return ensemble\n\ndef topology_aware_clean(mask):\n    cc_structure = ndimage.generate_binary_structure(rank=3, connectivity=1)\n    labeled, num = ndimage.label(mask, structure=cc_structure)\n    \n    if num == 0:\n        return mask.astype(np.uint8)\n    \n    component_sizes = ndimage.sum(mask, labeled, index=np.arange(1, num + 1))\n    keep_labels = np.where(component_sizes >= MIN_COMPONENT_VOXELS)[0] + 1\n    cleaned = np.isin(labeled, keep_labels).astype(np.uint8)\n    \n    if CLOSING_ITERS > 0:\n        closing_structure = np.zeros((3, 3, 3), dtype=np.uint8)\n        closing_structure[1, :, :] = 1\n        cleaned = ndimage.binary_closing(cleaned, structure=closing_structure, iterations=CLOSING_ITERS).astype(np.uint8)\n    \n    return cleaned.astype(np.float32)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-20T10:04:57.129272Z","iopub.execute_input":"2025-11-20T10:04:57.130000Z","iopub.status.idle":"2025-11-20T10:04:57.140414Z","shell.execute_reply.started":"2025-11-20T10:04:57.129970Z","shell.execute_reply":"2025-11-20T10:04:57.139739Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"##  4: 🚀 Brewing Our Perfect Cold Coffee (Main Pipeline)\n\nTime to brew! Load the model, grab the test data, and run our coffee making factory to create the submission.zip!","metadata":{}},{"cell_type":"code","source":"print(\"\\n 1️⃣/3️⃣Loading model\")\ndevice = torch.device('cuda' if torch.cuda.is_available() else 'cpu')#mill jaye to gpu na mille to cpu\nmodel = UNet3D(in_channels=1, out_channels=1, init_features=24).to(device)\n\nmodel_path = Path('/kaggle/input/model11/best_model_fast.pth')\nif model_path.exists():\n    model.load_state_dict(torch.load(model_path, map_location=device))\n    print(f\"Loaded: {model_path}\")\n\nmodel.eval()\n\nprint(\"\\n[2️⃣/3️⃣] Now test data\")\nROOT = Path('/kaggle/input/vesuvius-challenge-surface-detection')\ntest_df = pd.read_csv(ROOT / 'test.csv')\ntest_imgs = ROOT / 'test_images'\nprint(f\"Test volumes: {len(test_df)}\")\n\nprint(\"\\n[3️⃣/3️⃣] Just 30 sec everthing completed --> Slice-adaptive + DL + Topology\")\n\nsubmission_zip = Path('/kaggle/working/submission.zip')\n\nwith zipfile.ZipFile(submission_zip, 'w', compression=zipfile.ZIP_DEFLATED, compresslevel=9) as zf:\n    for idx, row in tqdm(test_df.iterrows(), total=len(test_df), desc='Processing'):\n        image_id = row['id']\n        img_path = test_imgs / f'{image_id}.tif'\n        \n        if not img_path.exists():\n            continue\n        \n        try:\n            volume = load_tiff(img_path).astype(np.float32)\n            \n            classical_pred = classical_cv_slice_adaptive(volume)\n            \n            dl_pred = dl_inference_fast(volume, model, device, stride=72)\n            \n            ensemble = slice_adaptive_ensemble(classical_pred, dl_pred)\n            \n            binary_mask = (ensemble > 0.5).astype(np.uint8)\n            final = topology_aware_clean(binary_mask)\n            \n            save_to_zip(final, zf, f'{image_id}.tif')\n            \n        except Exception as e:\n            print(f\"ERROR {image_id}: {str(e)[:80]}\")\n\nzip_size = submission_zip.stat().st_size / 1e6\n\n\nprint(f\"\\n{'_-'*60}\")\nprint(f\"DONE! submission.zip ({zip_size:.1f} MB)\")\nprint(f\"{'_-'*60}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-20T10:09:34.337875Z","iopub.execute_input":"2025-11-20T10:09:34.338722Z","iopub.status.idle":"2025-11-20T10:09:59.457516Z","shell.execute_reply.started":"2025-11-20T10:09:34.338691Z","shell.execute_reply":"2025-11-20T10:09:59.456827Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}