{"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,"isSourceIdPinned":false,"sourceType":"competition"}],"dockerImageVersionId":31192,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\n#for dirname, _, filenames in os.walk('/kaggle/input'):\n#    for filename in filenames:\n#        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-11-15T12:56:07.171901Z","iopub.execute_input":"2025-11-15T12:56:07.172125Z","iopub.status.idle":"2025-11-15T12:56:07.885142Z","shell.execute_reply.started":"2025-11-15T12:56:07.172102Z","shell.execute_reply":"2025-11-15T12:56:07.884578Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Install packages","metadata":{}},{"cell_type":"code","source":"!pip install --quiet imagecodecs\n!pip install segmentation-models-pytorch","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-15T12:56:07.885780Z","iopub.execute_input":"2025-11-15T12:56:07.886018Z","iopub.status.idle":"2025-11-15T12:56:15.383905Z","shell.execute_reply.started":"2025-11-15T12:56:07.886003Z","shell.execute_reply":"2025-11-15T12:56:15.382965Z"},"_kg_hide-output":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Load datasets","metadata":{}},{"cell_type":"code","source":"%%time\nfrom pathlib import Path\nDATA_ROOT = Path('/kaggle/input/vesuvius-challenge-surface-detection')\ntrain_csv = '/kaggle/input/vesuvius-challenge-surface-detection/train.csv'\ntest_csv = '/kaggle/input/vesuvius-challenge-surface-detection/test.csv'\ntest_image_dir = '/kaggle/input/vesuvius-challenge-surface-detection/test_images'\n\ntrain_images_dir = DATA_ROOT / 'train_images'\ntrain_labels_dir = DATA_ROOT / 'train_labels'","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-15T12:56:15.385136Z","iopub.execute_input":"2025-11-15T12:56:15.385428Z","iopub.status.idle":"2025-11-15T12:56:15.391179Z","shell.execute_reply.started":"2025-11-15T12:56:15.385389Z","shell.execute_reply":"2025-11-15T12:56:15.390469Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# For inline plots\n%matplotlib inline\n\n# Install required packages quietly\n!pip install -q tifffile imagecodecs ipywidgets","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-15T12:56:15.392053Z","iopub.execute_input":"2025-11-15T12:56:15.392277Z","iopub.status.idle":"2025-11-15T12:56:18.766828Z","shell.execute_reply.started":"2025-11-15T12:56:15.392252Z","shell.execute_reply":"2025-11-15T12:56:18.765984Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile patch_dataset.py\n# patch_dataset.py\nimport os\nfrom pathlib import Path\nimport random\nimport numpy as np\nimport tifffile\nimport scipy.ndimage as ndi\n\nimport torch\nfrom torch.utils.data import Dataset, DataLoader\n\n# -------------------------\n# Utility helpers\n# -------------------------\ndef load_tiff(path):\n    \"\"\"Load a single 3D tif as numpy array (Z, H, W).\"\"\"\n    arr = tifffile.imread(path)\n    # ensure numpy contiguous and standard dtype (float32 for images)\n    return np.asarray(arr)\n\ndef ensure_dir(p):\n    os.makedirs(p, exist_ok=True)\n\n# -------------------------\n# Surface-focused sampler\n# -------------------------\ndef compute_surface_weights(label_volume, cache_path=None, sigma=5.0):\n    \"\"\"\n    Compute a per-voxel sampling weight that favors voxels near the surface.\n    - label_volume: 3D numpy int array with values {0,1,2} (2 = unlabeled)\n    - returns weights normalized to sum=1 over labeled voxels (where label != 2).\n    \"\"\"\n    # mask only over labeled voxels (0 and 1)\n    labeled_mask = (label_volume != 2)\n    if not labeled_mask.any():\n        # fallback: uniform weights\n        weights = np.ones_like(label_volume, dtype=np.float32)\n        weights[:] = 1.0\n        if cache_path:\n            np.save(cache_path, weights)\n        return weights.astype(np.float32)\n\n    # distance transform from foreground boundary:\n    # compute distance to nearest boundary between foreground and background\n    fg = (label_volume == 1)\n    bg = (label_volume == 0)\n\n    # compute distance to nearest foreground voxel for every voxel (only meaningful for labeled voxels)\n    # use signed distance: distance to nearest foreground - distance to nearest background\n    dist_to_fg = ndi.distance_transform_edt(~fg)\n    dist_to_bg = ndi.distance_transform_edt(~bg)\n    signed_dist = dist_to_fg - dist_to_bg  # positive inside fg, negative in bg\n\n    # absolute distance to surface\n    abs_dist = np.abs(signed_dist)\n\n    # weight = exp(-abs_dist / sigma) but only for labeled voxels\n    weights = np.zeros_like(abs_dist, dtype=np.float32)\n    weights[labeled_mask] = np.exp(-abs_dist[labeled_mask] / float(sigma))\n\n    # small smoothing + ensure not all zeros\n    if weights.sum() == 0:\n        weights[labeled_mask] = 1.0\n\n    # normalize to sum=1 over full volume (so sampling using probabilities is straightforward)\n    weights = weights / weights.sum()\n    if cache_path:\n        np.save(cache_path, weights)\n    return weights.astype(np.float32)\n\n# -------------------------\n# Vesuvius Patch Dataset\n# -------------------------\nclass VesuviusPatchDataset(Dataset):\n    \"\"\"\n    Dataset that returns 3D patches sampled from volumes with a bias towards surface voxels.\n    - images: list of image file paths (.tif)\n    - labels: list of label file paths (.tif) or None (for unlabeled/test)\n    - patch_size: tuple (pz, ph, pw)\n    - samples_per_epoch: approximate number of patches returned per epoch (dataset length)\n    - cache_dir: where to store precomputed sampling weight arrays (npz / .npy)\n    - transforms: optional callable(image, label) -> (image, label)\n    \"\"\"\n    def __init__(self, images, labels, patch_size=(16, 256, 256),\n                 samples_per_epoch=4000, cache_dir='./cache',\n                 sigma=3.0, transforms=None, min_labeled_fraction=0.01):\n        super().__init__()\n        assert len(images) == len(labels)\n        self.images = [str(p) for p in images]\n        self.labels = [str(p) for p in labels]\n        self.patch_size = tuple(int(x) for x in patch_size)\n        self.samples_per_epoch = int(samples_per_epoch)\n        self.transforms = transforms\n        self.cache_dir = Path(cache_dir)\n        ensure_dir(self.cache_dir)\n        self.sigma = sigma\n        self.min_labeled_fraction = min_labeled_fraction\n\n        # lazy metadata cache\n        self._meta = [None] * len(self.images)\n        self._weights_paths = [self.cache_dir / f\"weights_{Path(p).stem}.npy\" for p in self.images]\n\n        # Precompute shapes (quick)\n        for i, img_path in enumerate(self.images):\n            arr = tifffile.imread(img_path, memmap=True)\n            self._meta[i] = {'shape': arr.shape, 'dtype': arr.dtype}\n            del arr\n\n    def __len__(self):\n        return self.samples_per_epoch\n\n    def _get_weights(self, idx):\n        wpath = self._weights_paths[idx]\n        if wpath.exists():\n            try:\n                weights = np.load(wpath)\n                return weights\n            except Exception:\n                pass\n        # compute\n        label = load_tiff(self.labels[idx])\n        weights = compute_surface_weights(label, cache_path=str(wpath), sigma=self.sigma)\n        return weights\n\n    def _sample_center(self, weights, shape):\n        \"\"\"\n        Sample a center coordinate for a patch based on weights probability map.\n        Ensures the patch fits inside the volume.\n        \"\"\"\n        pz, ph, pw = self.patch_size\n        Z, H, W = shape\n        # flatten and sample index\n        flat = weights.ravel()\n        # in case of nan or zeros, fallback to uniform over labeled area\n        if not np.isfinite(flat).all() or flat.sum() == 0:\n            flat = None\n\n        if flat is not None:\n            flat = flat / flat.sum()\n            idx = np.random.choice(flat.size, p=flat)\n            z0 = idx // (H * W)\n            rem = idx % (H * W)\n            y0 = rem // W\n            x0 = rem % W\n        else:\n            # uniform\n            z0 = np.random.randint(0, Z)\n            y0 = np.random.randint(0, H)\n            x0 = np.random.randint(0, W)\n\n        # convert to patch top-left coordinates while clamping to volume borders\n        z = int(np.clip(z0 - pz // 2, 0, max(0, Z - pz)))\n        y = int(np.clip(y0 - ph // 2, 0, max(0, H - ph)))\n        x = int(np.clip(x0 - pw // 2, 0, max(0, W - pw)))\n        return z, y, x\n\n    def __getitem__(self, idx):\n        # pick a random volume (uniform across volumes). Alternative: choose volumes proportional to labeled volume size.\n        vol_idx = random.randint(0, len(self.images) - 1)\n\n        img_path = self.images[vol_idx]\n        lbl_path = self.labels[vol_idx]\n\n        # load memmap / full arrays\n        image = load_tiff(img_path).astype(np.float32)\n        label = load_tiff(lbl_path).astype(np.int16)\n\n        # ensure shape matches meta\n        Z, H, W = image.shape\n\n        # get / compute weights\n        weights = self._get_weights(vol_idx)\n\n        # ensure min labeled fraction: if volume has almost no labeled voxels, fall back to uniform sampling\n        if (label != 2).sum() / label.size < self.min_labeled_fraction:\n            weights = np.ones_like(weights, dtype=np.float32)\n            weights = weights / weights.sum()\n\n        # sample center then crop\n        z, y, x = self._sample_center(weights, (Z, H, W))\n        pz, ph, pw = self.patch_size\n        img_patch = image[z:z + pz, y:y + ph, x:x + pw]\n        lbl_patch = label[z:z + pz, y:y + ph, x:x + pw]\n\n        # optional transforms (numpy in -> numpy out)\n        if self.transforms is not None:\n            img_patch, lbl_patch = self.transforms(img_patch, lbl_patch)\n\n        # convert to tensors: channel-first convention for 3D -> (C=1, Z, H, W)\n        img_patch = torch.from_numpy(img_patch).unsqueeze(0)  # 1, Z, H, W\n        lbl_patch = torch.from_numpy(lbl_patch).long()\n\n        # map unlabeled (2) to ignore index e.g., 255 for loss functions\n        lbl_patch[lbl_patch == 2] = 255  # recommended ignore index for CrossEntropy / custom losses\n\n        return img_patch, lbl_patch\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-15T12:56:18.770329Z","iopub.execute_input":"2025-11-15T12:56:18.770621Z","iopub.status.idle":"2025-11-15T12:56:18.779061Z","shell.execute_reply.started":"2025-11-15T12:56:18.770591Z","shell.execute_reply":"2025-11-15T12:56:18.778468Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\nimport tifffile\nimport numpy as np\nfrom torch.utils.data import Dataset\n\nclass VesuviusPatchDataset(Dataset):\n    def __init__(self, volume_path, mask_path=None, patch_size=(64,128,128)):\n        self.volume_path = volume_path\n        self.mask_path = mask_path\n        self.patch_size = patch_size\n\n        # load as memmap (safe for huge tif)\n        self.volume = tifffile.memmap(self.volume_path)\n\n        if mask_path is not None:\n            self.mask = tifffile.memmap(self.mask_path)\n        else:\n            self.mask = None\n\n        self.D, self.H, self.W = self.volume.shape\n        self.pd, self.ph, self.pw = patch_size\n\n    def __len__(self):\n        return 2000  # random sampling\n\n    def __getitem__(self, idx):\n        # random patch\n        z = np.random.randint(0, self.D - self.pd)\n        y = np.random.randint(0, self.H - self.ph)\n        x = np.random.randint(0, self.W - self.pw)\n\n        v = self.volume[z:z+self.pd, y:y+self.ph, x:x+self.pw].astype(np.float32)\n        v = (v - v.mean()) / (v.std() + 1e-6)\n        v = np.expand_dims(v, 0)  # [C=1, D, H, W]\n\n        if self.mask is not None:\n            m = self.mask[z:z+self.pd, y:y+self.ph, x:x+self.pw]\n            m = (m == 1).astype(np.float32)\n            m = np.expand_dims(m, 0)  # [1, D, H, W]\n        else:\n            m = np.zeros_like(v)\n\n        return v, m\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-15T12:56:18.779819Z","iopub.execute_input":"2025-11-15T12:56:18.780207Z","iopub.status.idle":"2025-11-15T12:56:20.493619Z","shell.execute_reply.started":"2025-11-15T12:56:18.780179Z","shell.execute_reply":"2025-11-15T12:56:20.492774Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\nfrom patch_dataset import VesuviusPatchDataset","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-15T12:56:20.494487Z","iopub.execute_input":"2025-11-15T12:56:20.494904Z","iopub.status.idle":"2025-11-15T12:56:20.830906Z","shell.execute_reply.started":"2025-11-15T12:56:20.494878Z","shell.execute_reply":"2025-11-15T12:56:20.830199Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# EDA","metadata":{}},{"cell_type":"markdown","source":"EDA snippet for Vesuvius volumes:","metadata":{}},{"cell_type":"code","source":"%%time\nimport tifffile\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom pathlib import Path\n\n#DATA_ROOT = Path('/kaggle/input/vesuvius-challenge-surface-detection')\n#train_images_dir = DATA_ROOT / 'train_images'\n#train_labels_dir = DATA_ROOT / 'train_labels'\n\n# pick a sample volume\nsample_img_path = sorted(train_images_dir.glob('*.tif'))[0]\nsample_lbl_path = train_labels_dir / sample_img_path.name\n\n# load\nimg = tifffile.imread(sample_img_path).astype(np.float32)\nlbl = tifffile.imread(sample_lbl_path).astype(np.int16)\n\nprint(f\"Volume shape: {img.shape}, dtype: {img.dtype}\")\nprint(f\"Label unique values: {np.unique(lbl)}\")\nprint(f\"Foreground voxels: {(lbl==1).sum()}, Background: {(lbl==0).sum()}, Unlabeled: {(lbl==2).sum()}\")\n\n# normalize image for display\nimg_norm = (img - img.min()) / (img.max() - img.min())\n\n# show few axial slices\nfig, axes = plt.subplots(2, 5, figsize=(20,6))\nslices = np.linspace(0, img.shape[0]-1, 10, dtype=int)\nfor ax, sl in zip(axes.flatten(), slices):\n    ax.imshow(img_norm[sl], cmap='gray')\n    ax.imshow(lbl[sl], cmap='jet', alpha=0.3)  # overlay mask\n    ax.set_title(f\"Slice {sl}\")\n    ax.axis('off')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-15T12:56:20.832007Z","iopub.execute_input":"2025-11-15T12:56:20.832344Z","iopub.status.idle":"2025-11-15T12:56:22.890031Z","shell.execute_reply.started":"2025-11-15T12:56:20.832324Z","shell.execute_reply":"2025-11-15T12:56:22.889195Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Here’s a ready-to-run surface distribution visualization:","metadata":{}},{"cell_type":"code","source":"%%time\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom pathlib import Path\nimport tifffile\n\n#DATA_ROOT = Path('/kaggle/input/vesuvius-challenge-surface-detection')\n#train_labels_dir = DATA_ROOT / 'train_labels'\n\nlabel_files = sorted(list(train_labels_dir.glob('*.tif')))\n\n# Compute fraction of surface voxels per volume\nsurface_fractions = []\n\nfor f in label_files[:10]:  # sample first 10 for quick check\n    #lbl = tifffile.memmap(f)\n    lbl = tifffile.imread(f).astype(np.int16)\n\n    fg = (lbl == 1)\n    bg = (lbl == 0)\n    \n    # compute approximate surface voxels: voxels with foreground neighbor in 6-connectivity\n    from scipy.ndimage import convolve\n    kernel = np.zeros((3,3,3), dtype=int)\n    kernel[1,1,0] = kernel[1,1,2] = 1\n    kernel[1,0,1] = kernel[1,2,1] = 1\n    kernel[0,1,1] = kernel[2,1,1] = 1\n    \n    neighbor_count = convolve(fg.astype(int), kernel, mode='constant', cval=0)\n    surface = fg & (neighbor_count < 6)  # voxels at boundary\n    surface_frac = surface.sum() / max(fg.sum(),1)\n    surface_fractions.append(surface_frac)\n\n# Plot distribution\nplt.figure(figsize=(8,4))\nplt.bar(range(len(surface_fractions)), surface_fractions)\nplt.xlabel(\"Volume index\")\nplt.ylabel(\"Surface fraction (approx.)\")\nplt.title(\"Approximate fraction of surface voxels per volume\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-15T12:56:22.891076Z","iopub.execute_input":"2025-11-15T12:56:22.891609Z","iopub.status.idle":"2025-11-15T12:56:31.579198Z","shell.execute_reply.started":"2025-11-15T12:56:22.891581Z","shell.execute_reply":"2025-11-15T12:56:31.578588Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\nimport os\nimport pandas as pd\nimport tifffile as tiff\nimport matplotlib.pyplot as plt\nimport numpy as np\n\n\ndf = pd.read_csv(train_csv)\nprint(df.head())\n\n# Pick 1-2 volumes for quick exploration\nsample_ids = df['id'].tolist()[:2]\n\nfor image_id in sample_ids:\n    image_path = os.path.join(train_images_dir, f\"{image_id}.tif\")\n    label_path = os.path.join(train_labels_dir, f\"{image_id}.tif\")\n    \n    vol = tiff.imread(image_path)\n    mask = tiff.imread(label_path)\n    \n    print(f\"Volume {image_id} shape: {vol.shape}, mask shape: {mask.shape}\")\n    print(f\"Mask unique values: {np.unique(mask)}\")\n    \n    # Show middle slice\n    mid_slice = vol.shape[0] // 2\n    plt.figure(figsize=(10,4))\n    plt.subplot(1,2,1)\n    plt.imshow(vol[mid_slice], cmap='gray')\n    plt.title(f'Volume {image_id} mid slice')\n    plt.subplot(1,2,2)\n    plt.imshow(mask[mid_slice], cmap='gray')\n    plt.title(f'Mask {image_id} mid slice')\n    plt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-15T12:56:31.579934Z","iopub.execute_input":"2025-11-15T12:56:31.580167Z","iopub.status.idle":"2025-11-15T12:56:34.119381Z","shell.execute_reply.started":"2025-11-15T12:56:31.580145Z","shell.execute_reply":"2025-11-15T12:56:34.118580Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Build a lightweight 3D UNet suitable for your patch-based training.","metadata":{}},{"cell_type":"code","source":"%%time\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\n\nclass ConvBlock3D(nn.Module):\n    def __init__(self, in_ch, out_ch, kernel_size=3, padding=1):\n        super().__init__()\n        self.conv = nn.Sequential(\n            nn.Conv3d(in_ch, out_ch, kernel_size, padding=padding),\n            nn.BatchNorm3d(out_ch),\n            nn.ReLU(inplace=True),\n            nn.Conv3d(out_ch, out_ch, kernel_size, padding=padding),\n            nn.BatchNorm3d(out_ch),\n            nn.ReLU(inplace=True),\n        )\n    def forward(self, x):\n        return self.conv(x)\n\nclass DownBlock3D(nn.Module):\n    def __init__(self, in_ch, out_ch):\n        super().__init__()\n        self.pool = nn.MaxPool3d(2)\n        self.conv = ConvBlock3D(in_ch, out_ch)\n    def forward(self, x):\n        x = self.pool(x)\n        x = self.conv(x)\n        return x\n\nclass UpBlock3D(nn.Module):\n    def __init__(self, in_ch, out_ch):\n        super().__init__()\n        self.up = nn.ConvTranspose3d(in_ch, out_ch, kernel_size=2, stride=2)\n        self.conv = ConvBlock3D(in_ch, out_ch)\n    def forward(self, x, skip):\n        x = self.up(x)\n        # pad if needed\n        if x.shape != skip.shape:\n            diffZ = skip.size(2) - x.size(2)\n            diffY = skip.size(3) - x.size(3)\n            diffX = skip.size(4) - x.size(4)\n            x = F.pad(x, [diffX//2, diffX - diffX//2,\n                          diffY//2, diffY - diffY//2,\n                          diffZ//2, diffZ - diffZ//2])\n        x = torch.cat([skip, x], dim=1)\n        x = self.conv(x)\n        return x\n\nclass UNet3D(nn.Module):\n    def __init__(self, in_channels=1, out_channels=1, features=[16, 32, 64]):\n        super().__init__()\n        \n        # Encoder\n        self.enc1 = self._block(in_channels, features[0])\n        self.pool1 = nn.MaxPool3d(2, 2)\n        \n        self.enc2 = self._block(features[0], features[1])\n        self.pool2 = nn.MaxPool3d(2, 2)\n        \n        self.enc3 = self._block(features[1], features[2])\n        self.pool3 = nn.MaxPool3d(2, 2)\n        \n        # Bottleneck\n        self.bottleneck = self._block(features[2], features[2]*2)\n        \n        # Decoder\n        self.upconv3 = nn.ConvTranspose3d(features[2]*2, features[2], 2, 2)\n        self.dec3 = self._block(features[2]*2, features[2])\n        \n        self.upconv2 = nn.ConvTranspose3d(features[2], features[1], 2, 2)\n        self.dec2 = self._block(features[1]*2, features[1])\n        \n        self.upconv1 = nn.ConvTranspose3d(features[1], features[0], 2, 2)\n        self.dec1 = self._block(features[0]*2, features[0])\n        \n        self.final = nn.Conv3d(features[0], out_channels, 1)\n    \n    def _block(self, in_channels, out_channels):\n        return nn.Sequential(\n            nn.Conv3d(in_channels, out_channels, 3, padding=1, bias=False),\n            nn.BatchNorm3d(out_channels),\n            nn.ReLU(inplace=True),\n            nn.Conv3d(out_channels, out_channels, 3, padding=1, bias=False),\n            nn.BatchNorm3d(out_channels),\n            nn.ReLU(inplace=True)\n        )\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        \n        b = self.bottleneck(self.pool3(e3))\n        \n        d3 = self.upconv3(b)\n        d3 = torch.cat((d3, e3), dim=1)\n        d3 = self.dec3(d3)\n        \n        d2 = self.upconv2(d3)\n        d2 = torch.cat((d2, e2), dim=1)\n        d2 = self.dec2(d2)\n        \n        d1 = self.upconv1(d2)\n        d1 = torch.cat((d1, e1), dim=1)\n        d1 = self.dec1(d1)\n        \n        return d1  # logits only\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-15T12:56:34.120438Z","iopub.execute_input":"2025-11-15T12:56:34.120718Z","iopub.status.idle":"2025-11-15T12:56:34.131879Z","shell.execute_reply.started":"2025-11-15T12:56:34.120699Z","shell.execute_reply":"2025-11-15T12:56:34.131100Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Minimal example to test the UNet3D on a patch from your dataset:","metadata":{}},{"cell_type":"code","source":"%%time\nimport tifffile\nimport numpy as np\nfrom torch.utils.data import Dataset\n\nclass VesuviusPatchDataset(Dataset):\n    def __init__(self, volume_path, mask_path=None, patch_size=(16,128,128)):\n        \"\"\"\n        volume_path: path to 3D tif volume\n        mask_path: path to 3D tif mask\n        patch_size: (D,H,W) size of 3D patch\n        \"\"\"\n        self.volume_path = volume_path\n        self.mask_path = mask_path\n        self.patch_size = patch_size\n\n        # load volume & mask fully (or use memmap if uncompressed)\n        self.volume = tifffile.imread(self.volume_path).astype(np.float32)\n        if self.mask_path is not None:\n            self.mask = tifffile.imread(self.mask_path).astype(np.int16)\n        else:\n            self.mask = None\n\n        self.D, self.H, self.W = self.volume.shape\n        self.pd, self.ph, self.pw = patch_size\n\n    def __len__(self):\n        return 2000  # number of patches per volume\n\n    def __getitem__(self, idx):\n        # random patch\n        z = np.random.randint(0, self.D - self.pd)\n        y = np.random.randint(0, self.H - self.ph)\n        x = np.random.randint(0, self.W - self.pw)\n\n        vol_patch = self.volume[z:z+self.pd, y:y+self.ph, x:x+self.pw].astype(np.float32)\n        vol_patch = (vol_patch - vol_patch.mean()) / (vol_patch.std() + 1e-6)\n        vol_patch = np.expand_dims(vol_patch, 0)  # [C=1, D,H,W]\n\n        if self.mask is not None:\n            mask_patch = self.mask[z:z+self.pd, y:y+self.ph, x:x+self.pw]\n            mask_patch = (mask_patch == 1).astype(np.float32)\n            mask_patch = np.expand_dims(mask_patch, 0)  # [1,D,H,W]\n        else:\n            mask_patch = np.zeros_like(vol_patch)\n\n        return vol_patch, mask_patch","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-15T12:56:34.132969Z","iopub.execute_input":"2025-11-15T12:56:34.133240Z","iopub.status.idle":"2025-11-15T12:56:34.148644Z","shell.execute_reply.started":"2025-11-15T12:56:34.133216Z","shell.execute_reply":"2025-11-15T12:56:34.147999Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Model Training / Prediction (Patch-based, quick test)\n\nTrain on only 1-2 volumes with few patches\n\nTest that the training loop runs without crashing\n\nPredict on the same volumes → visualize predictions\n\nSub function","metadata":{}},{"cell_type":"code","source":"%%time\nclass ScrollPatchDataset(Dataset):\n    \"\"\"\n    Patch-based dataset for 3D scroll volumes.\n    Extracts random patches from each volume for training.\n    \"\"\"\n    def __init__(self, image_dir, label_dir, csv_file, patch_size=(16,128,128), max_patches_per_volume=10):\n        self.image_dir = image_dir\n        self.label_dir = label_dir\n        self.df = pd.read_csv(csv_file)\n        self.patch_size = patch_size\n        self.max_patches_per_volume = max_patches_per_volume\n\n    def __len__(self):\n        return len(self.df) * self.max_patches_per_volume\n\n    def load_tiff(self, path):\n        try:\n            return tiff.imread(path)\n        except:\n            import imagecodecs\n            return tiff.imread(path)\n\n    def __getitem__(self, idx):\n        # idx -> volume_idx + patch_idx\n        volume_idx = idx // self.max_patches_per_volume\n        image_id = str(self.df.iloc[volume_idx]['id'])\n\n        image_path = os.path.join(self.image_dir, f\"{image_id}.tif\")\n        label_path = os.path.join(self.label_dir, f\"{image_id}.tif\")\n\n        vol = self.load_tiff(image_path)\n        mask = self.load_tiff(label_path)\n\n        # Add channel dim if 2D\n        if len(vol.shape) == 2:\n            vol = vol[np.newaxis, :, :]\n            mask = mask[np.newaxis, :, :]\n\n        # Normalize volume\n        vol = vol.astype(np.float32)\n        vol = (vol - vol.mean()) / (vol.std() + 1e-8)\n\n        # Normalize mask to 0-1\n        mask = mask.astype(np.float32)\n        if mask.max() > 1.0:\n            mask = mask / 255.0\n        mask = np.clip(mask, 0.0, 1.0)\n\n        # Random crop patch\n        D, H, W = vol.shape\n        pd, ph, pw = self.patch_size\n\n        d_start = np.random.randint(0, max(1, D - pd + 1))\n        h_start = np.random.randint(0, max(1, H - ph + 1))\n        w_start = np.random.randint(0, max(1, W - pw + 1))\n\n        vol_patch = vol[d_start:d_start+pd, h_start:h_start+ph, w_start:w_start+pw]\n        mask_patch = mask[d_start:d_start+pd, h_start:h_start+ph, w_start:w_start+pw]\n\n        # Convert to torch [C,D,H,W]\n        vol_patch = torch.from_numpy(vol_patch).float().unsqueeze(0)\n        mask_patch = torch.from_numpy(mask_patch).float().unsqueeze(0)\n\n        return vol_patch, mask_patch","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-15T12:56:34.149357Z","iopub.execute_input":"2025-11-15T12:56:34.149643Z","iopub.status.idle":"2025-11-15T12:56:34.165835Z","shell.execute_reply.started":"2025-11-15T12:56:34.149616Z","shell.execute_reply":"2025-11-15T12:56:34.165122Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## ✅ Changes for fast CPU run:\n\nPatch size reduced: (8, 64, 64) → 8 slices, 64×64 pixels per patch.\n\nMax patches per volume: 1 → fewer samples per volume.\n\nUNet features: [2, 4, 8] → smaller model, faster convs.\n\nEpochs: 1 → test end-to-end pipeline.\n","metadata":{}},{"cell_type":"code","source":"%%time\n# =======================\n# CPU-friendly debug settings\n# =======================\nimport torch.optim as optim\nfrom torch.utils.data import Dataset, DataLoader\nfrom tqdm import tqdm\n\npatch_size = (8, 64, 64)\nmax_patches_per_volume = 1\nfeatures = [1, 4, 8] # [2, 4, 8]\nepochs = 1\n\n\n# =======================\n# Device setup\n# =======================\ndevice = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\nprint(f\"Using device: {device}\")\n\n# =======================\n# Loss\n# =======================\ndef dice_loss(pred, target, smooth=1.):\n    pred = torch.sigmoid(pred)\n    intersection = (pred * target).sum()\n    return 1 - (2. * intersection + smooth) / (pred.sum() + target.sum() + smooth)\n\ndef combined_loss(pred, target):\n    bce = F.binary_cross_entropy_with_logits(pred, target)\n    dice = dice_loss(pred, target)\n    return bce + dice\n    \n# Dataset & DataLoader\ndataset = ScrollPatchDataset(\n    train_images_dir, train_labels_dir, train_csv,\n    patch_size=patch_size, max_patches_per_volume=max_patches_per_volume\n)\ntrain_loader = DataLoader(dataset, batch_size=1, shuffle=True, num_workers=2, pin_memory=True)\nprint(f\"Debug training on {len(dataset)} patches\")\n\n# Smaller UNet3D\nmodel = UNet3D(in_channels=1, out_channels=1, features=features).to(device)\n\n\noptimizer = optim.Adam(model.parameters(), lr=1e-4, weight_decay=1e-5)\nscheduler = optim.lr_scheduler.ReduceLROnPlateau(optimizer, 'min', patience=2, factor=0.5)\n\n# =======================\n# Quick training loop\n# =======================\nmodel.train()\nbest_loss = float('inf')\nfor epoch in range(epochs):\n    epoch_loss = 0\n    progress = tqdm(train_loader, desc=f\"Epoch {epoch+1}/{epochs}\")\n    for data, target in progress:\n        data, target = data.to(device), target.to(device)\n        optimizer.zero_grad()\n        output = model(data)\n        loss = combined_loss(output, target)\n        loss.backward()\n        torch.nn.utils.clip_grad_norm_(model.parameters(), 1.0)\n        optimizer.step()\n        epoch_loss += loss.item()\n        progress.set_postfix({\"loss\": f\"{loss.item():.4f}\"})\n    avg_loss = epoch_loss / len(train_loader)\n    scheduler.step(avg_loss)\n    print(f\"Epoch {epoch+1} finished, avg loss: {avg_loss:.4f}\")\n\n    if avg_loss < best_loss:\n        best_loss = avg_loss\n        torch.save(model.state_dict(), '/kaggle/working/best_model_debug.pth')\n        print(f\"Saved best debug model with loss {best_loss:.4f}\")\n\nprint(\"CPU debug training completed!\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-15T12:56:34.166629Z","iopub.execute_input":"2025-11-15T12:56:34.167095Z","iopub.status.idle":"2025-11-15T13:02:51.239618Z","shell.execute_reply.started":"2025-11-15T12:56:34.167075Z","shell.execute_reply":"2025-11-15T13:02:51.238840Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"##  Prediction & Submission (Debug Mode)","metadata":{}},{"cell_type":"code","source":"%%time\n# this funciton error, pending\nimport torch\nimport numpy as np\nimport pandas as pd\nimport tifffile as tiff\nfrom torch.utils.data import Dataset\nfrom skimage.transform import resize\n\nclass ScrollDataset(Dataset):\n    def __init__(self, image_dir, csv_file, max_depth=16, target_size=128):\n        \"\"\"\n        image_dir: path to test_images/\n        csv_file: path to test.csv\n        max_depth: number of slices to keep\n        target_size: H/W resize to this value\n        \"\"\"\n        self.image_dir = image_dir\n        self.df = pd.read_csv(csv_file)\n        self.max_depth = max_depth\n        self.target_size = target_size\n\n    def __len__(self):\n        return len(self.df)\n\n    def __getitem__(self, idx):\n        \n        image_id = str(self.df.iloc[idx]['id'])\n        image_path = os.path.join(self.image_dir, f\"{image_id}.tif\")\n\n        volume = self.load_tiff(image_path)\n\n        # Handle 2D images → convert to 3D\n        if len(volume.shape) == 2:\n            volume = volume[np.newaxis, :, :]\n\n        # Resize\n        if volume.shape != (self.max_depth, self.target_size, self.target_size):\n            volume = self.resize_volume(volume, self.max_depth, self.target_size, self.target_size)\n\n        volume = volume.astype(np.float32)\n        volume = (volume - volume.mean()) / (volume.std() + 1e-8)\n\n        volume = torch.from_numpy(volume).float().unsqueeze(0)\n        return volume, image_id","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-15T13:02:51.240672Z","iopub.execute_input":"2025-11-15T13:02:51.241102Z","iopub.status.idle":"2025-11-15T13:02:51.321326Z","shell.execute_reply.started":"2025-11-15T13:02:51.241070Z","shell.execute_reply":"2025-11-15T13:02:51.320764Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\nclass ScrollDataset(Dataset):\n    def __init__(self, image_dir, csv_file, max_depth=16, target_size=128):\n        self.image_dir = image_dir\n        self.df = pd.read_csv(csv_file)\n        self.max_depth = max_depth\n        self.target_size = target_size\n\n    def load_tiff(self, path):\n        try:\n            return tiff.imread(path)\n        except:\n            import imagecodecs\n            return tiff.imread(path)\n\n    def resize_volume(self, volume, target_depth, target_height, target_width):\n        resized = np.zeros((target_depth, target_height, target_width), dtype=volume.dtype)\n\n        depth_ratio = volume.shape[0] / target_depth\n        for d in range(target_depth):\n            orig_d = min(int(d * depth_ratio), volume.shape[0] - 1)\n            slice_2d = volume[orig_d]\n            resized_slice = cv2.resize(\n                slice_2d,\n                (target_width, target_height),\n                interpolation=cv2.INTER_AREA\n            )\n            resized[d] = resized_slice\n        return resized\n\n    def __getitem__(self, idx):\n        image_id = str(self.df.iloc[idx]['id'])\n        image_path = os.path.join(self.image_dir, f\"{image_id}.tif\")\n\n        volume = self.load_tiff(image_path)\n\n        # Ensure 3D\n        if len(volume.shape) == 2:\n            volume = volume[np.newaxis, :, :]\n\n        # Resize\n        if volume.shape != (self.max_depth, self.target_size, self.target_size):\n            volume = self.resize_volume(volume, self.max_depth, self.target_size, self.target_size)\n\n        # Normalize\n        volume = volume.astype(np.float32)\n        volume = (volume - volume.mean()) / (volume.std() + 1e-8)\n\n        volume = torch.from_numpy(volume).float().unsqueeze(0)\n\n        return volume, image_id\n\n    def __len__(self):\n        return len(self.df)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-15T13:02:51.322171Z","iopub.execute_input":"2025-11-15T13:02:51.322869Z","iopub.status.idle":"2025-11-15T13:02:51.330059Z","shell.execute_reply.started":"2025-11-15T13:02:51.322849Z","shell.execute_reply":"2025-11-15T13:02:51.329331Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## ✅ Prediction and Submission testing\n1. Sliding Window Predictor","metadata":{}},{"cell_type":"code","source":"%%time\nimport torch\nimport numpy as np\nimport tifffile as tiff\nfrom tqdm import tqdm\n\ndef predict_full_volume(model, volume, patch_size=128, stride=96):\n    \"\"\"\n    volume: numpy array [D, H, W]\n    returns: mask_full [H, W]\n    \"\"\"\n    model.eval()\n    D, H, W = volume.shape\n\n    # output mask accumulation\n    output = np.zeros((H, W), dtype=np.float32)\n    count  = np.zeros((H, W), dtype=np.float32)\n\n    # loop patches\n    for y in tqdm(range(0, H - patch_size + 1, stride), desc=\"Sliding Y\"):\n        for x in range(0, W - patch_size + 1, stride):\n            patch = volume[:, y:y+patch_size, x:x+patch_size].astype(np.float32)\n\n            # normalize patch\n            p = (patch - patch.mean()) / (patch.std() + 1e-8)\n            p = torch.from_numpy(p).float().unsqueeze(0).unsqueeze(0).to(device)  # [1,1,D,H,W]\n\n            with torch.no_grad():\n                pred = model(p).cpu().numpy()[0,0]  # [D,H,W]\n\n            # collapse depth by max\n            pred2d = pred.max(axis=0)\n\n            output[y:y+patch_size, x:x+patch_size] += pred2d\n            count[y:y+patch_size, x:x+patch_size] += 1\n\n    # avoid division by zero\n    count[count == 0] = 1\n    final_mask = output / count\n    final_mask = (final_mask * 255).clip(0,255).astype(np.uint8)\n\n    return final_mask","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-15T13:02:51.330845Z","iopub.execute_input":"2025-11-15T13:02:51.331179Z","iopub.status.idle":"2025-11-15T13:02:51.345993Z","shell.execute_reply.started":"2025-11-15T13:02:51.331155Z","shell.execute_reply":"2025-11-15T13:02:51.345470Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"2. Test Set Volume Loader","metadata":{}},{"cell_type":"code","source":"%%time\ndef load_test_volume(test_id):\n    path = f\"/kaggle/input/vesuvius-challenge-surface-detection/test_images/{test_id}.tif\"\n    vol = tiff.imread(path)\n\n    # ensure 3D\n    if len(vol.shape) == 2:\n        vol = vol[np.newaxis, :, :]\n\n    # reduce depth to match training depth\n    target_depth = 16\n    if vol.shape[0] > target_depth:\n        idx = np.linspace(0, vol.shape[0]-1, target_depth).astype(int)\n        vol = vol[idx]\n\n    return vol","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-15T13:02:51.346749Z","iopub.execute_input":"2025-11-15T13:02:51.347025Z","iopub.status.idle":"2025-11-15T13:02:51.361575Z","shell.execute_reply.started":"2025-11-15T13:02:51.347001Z","shell.execute_reply":"2025-11-15T13:02:51.360903Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"3. Run Prediction on All Test Samples","metadata":{}},{"cell_type":"code","source":"%%time\nimport zipfile\nimport os\n\nmodel_path = \"/kaggle/working/best_model_debug.pth\"\ncheckpoint = torch.load(model_path, map_location=device)\nmodel = UNet3D().to(device)\n#model.load_state_dict(checkpoint[\"model_state_dict\"])\n\nsubmission_dir = \"/kaggle/working/submission_full\"\nos.makedirs(submission_dir, exist_ok=True)\n\ntest_ids = [1407735]  # update if more test samples appear\nprint(\"Predicting test IDs:\", test_ids)\n\nfor tid in test_ids:\n    print(f\"\\n=== Predicting test_volume {tid} ===\")\n    vol = load_test_volume(tid)\n\n    mask = predict_full_volume(model, vol, patch_size=128, stride=96)\n\n    out_path = f\"{submission_dir}/{tid}.tif\"\n    tiff.imwrite(out_path, mask)\n    print(\"Saved:\", out_path)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-15T13:02:51.362386Z","iopub.execute_input":"2025-11-15T13:02:51.362746Z","iopub.status.idle":"2025-11-15T13:02:52.009188Z","shell.execute_reply.started":"2025-11-15T13:02:51.362701Z","shell.execute_reply":"2025-11-15T13:02:52.008528Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Create the final submission.zip and upload to Kaggle.","metadata":{}},{"cell_type":"code","source":"%%time\nimport zipfile\n\nzip_path = \"/kaggle/working/submission.zip\"\n\nwith zipfile.ZipFile(zip_path, 'w') as z:\n    z.write(\"/kaggle/working/submission_full/1407735.tif\", arcname=\"1407735.tif\")\n\nprint(\"Final submission ready:\", zip_path)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-15T13:02:52.010036Z","iopub.execute_input":"2025-11-15T13:02:52.010277Z","iopub.status.idle":"2025-11-15T13:02:52.015892Z","shell.execute_reply.started":"2025-11-15T13:02:52.010259Z","shell.execute_reply":"2025-11-15T13:02:52.015309Z"}},"outputs":[],"execution_count":null}]}