{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.12.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":117682,"databundleVersionId":14443416,"sourceType":"competition"},{"sourceId":13698445,"sourceType":"datasetVersion","datasetId":8670484}],"dockerImageVersionId":31234,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"\"\"\"\nVesuvius competition metric.\n\nExpects standard Kaggle paths and Linux in order to manage dependencies.\n\"\"\"\n\nimport glob\nimport importlib\nimport os\nimport subprocess\nimport sys\nimport numpy as np\nimport pandas as pd\nfrom PIL import Image, ImageSequence\nfrom scipy.ndimage import distance_transform_edt, binary_dilation, grey_dilation, gaussian_gradient_magnitude, map_coordinates, gaussian_filter, binary_fill_holes\nimport networkx as nx\nfrom scipy.sparse import lil_matrix\nfrom scipy.sparse.linalg import lsqr\nfrom scipy.sparse import csr_matrix\nfrom scipy.sparse.linalg import spsolve\nfrom pathlib import Path\nfrom scipy.sparse import coo_matrix\nfrom scipy.sparse.linalg import cg\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-12-20T03:19:00.871991Z","iopub.execute_input":"2025-12-20T03:19:00.872381Z","iopub.status.idle":"2025-12-20T03:19:01.842474Z","shell.execute_reply.started":"2025-12-20T03:19:00.872334Z","shell.execute_reply":"2025-12-20T03:19:01.841307Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class ParticipantVisibleError(Exception):\n    pass\n\n\nclass HostVisibleError(Exception):\n    pass\n\n\ndef load_volume(path):\n    im = Image.open(path)\n    slices = []\n    for i, page in enumerate(ImageSequence.Iterator(im)):\n        slice_array = np.array(page)\n        slices.append(slice_array)\n    volume = np.stack(slices, axis=0)\n    return volume\n\n\ndef install_dependencies():\n    \"\"\"On Kaggle, the topometrics library must be installed during the run. This function handles the entire process.\"\"\"\n    try:\n        import topometrics.leaderboard\n\n        return None\n    # The broad exception is necessary as the initial import can fail for multiple reasons.\n    except:\n        pass\n\n    resources_dir = '/kaggle/input/vesuvius-metric-resources'\n    install_dir = '/kaggle/working/topological-metrics-kaggle'\n\n    try:\n        subprocess.run(\n            f'cd {resources_dir} && uv pip install --no-index --find-links=wheels -r topological-metrics-kaggle/requirements.txt',\n            shell=True,\n            check=True,\n        )\n        subprocess.run(f'cd /kaggle/working && cp -r {resources_dir}/topological-metrics-kaggle .', shell=True, check=True)\n        subprocess.run(\n            f'cd {install_dir} && chmod +x scripts/setup_submodules.sh scripts/build_betti.sh && make build-betti',\n            shell=True,\n            check=True,\n        )\n        subprocess.run(\n            f'cd {install_dir} && uv pip install -e . --no-deps --no-index --no-build-isolation -v',\n            shell=True,\n            check=True,\n        )\n        # Add the new library to Python's path and invalidate caches to ensure it's found.\n        sys.path.append('/kaggle/working/topological-metrics-kaggle/src')\n        importlib.invalidate_caches()\n\n    except Exception as err:\n        raise HostVisibleError(f'Failed to install topometrics library: {err}')\n\n\ndef generate_standard_submission(submission_dir: str) -> None:\n    # Dependencies installed here as generate_standard_submission is the first metric function that gets called by the orchestrator.\n    submission_tifs = glob.glob(f'{submission_dir}/**/*.tif', recursive=True)\n    if len(submission_tifs) == 0:\n        submission_tifs = glob.glob('/kaggle/tmp/**/*.tif', recursive=True)\n    if len(submission_tifs) == 0:\n        raise ParticipantVisibleError('No submission files found')\n    df = pd.DataFrame({'tif_paths': submission_tifs})\n    df['id'] = df['tif_paths'].apply(lambda x: x.split('/')[-1].split('.')[0])\n    os.chdir('/kaggle/working')\n    df[['id', 'tif_paths']].to_csv('submission.csv', index=False)\n\n\ndef score_single_tif(\n    gt_path,\n    pred_path,\n    surface_tolerance,\n    voi_connectivity=26,\n    voi_transform='one_over_one_plus',\n    voi_alpha=0.3,\n    topo_weight=0.3,\n    surface_dice_weight=0.35,\n    voi_weight=0.35,\n):\n    gt: np.ndarray = load_volume(gt_path)\n    pr: np.ndarray = load_volume(pred_path)\n\n    install_dependencies()\n    # The import is here to ensure dependencies are loaded first.\n    try:\n        # Use a standard import now that the path is reliably set.\n        import topometrics.leaderboard\n    except Exception as err:\n        raise HostVisibleError(f'Failed to import topometrics after installation: {err}')\n\n    score_report = topometrics.leaderboard.compute_leaderboard_score(\n        predictions=pr,\n        labels=gt,\n        dims=(0, 1, 2),\n        spacing=(1.0, 1.0, 1.0),  # (z, y, x)\n        surface_tolerance=surface_tolerance,  # in spacing units\n        voi_connectivity=voi_connectivity,\n        voi_transform=voi_transform,\n        voi_alpha=voi_alpha,\n        combine_weights=(topo_weight, surface_dice_weight, voi_weight),  # (Topo, SurfaceDice, VOI)\n        fg_threshold=None,  # None => legacy \"!= 0\"; else uses \"x > threshold\"\n        ignore_label=2,  # voxels with this GT label are ignored\n        ignore_mask=None,  # or pass an explicit boolean mask\n    )\n    print(\"topo_score:\", score_report.topo.toposcore,\n        \"voi_score:\", score_report.voi.voi_score,\n        \"surface_dice:\", score_report.surface_dice)\n    return np.clip(score_report.score, a_min=0.0, a_max=1.0)\n\ninstall_dependencies()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-20T03:19:03.441883Z","iopub.execute_input":"2025-12-20T03:19:03.443375Z","iopub.status.idle":"2025-12-20T03:19:04.221581Z","shell.execute_reply.started":"2025-12-20T03:19:03.443321Z","shell.execute_reply":"2025-12-20T03:19:04.220327Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Decomposition A","metadata":{}},{"cell_type":"code","source":"class ExactScrollDecomposition:\n    \"\"\"\n    Exact scroll mask decomposition using:\n    - sdf: signed distance field\n    - normals: gradient of sdf\n    - thickness: distance to mask boundary along normal\n    Fully vectorized inverse transform using sparse accumulation.\n    \"\"\"\n\n    def __init__(self, sdf_truncate=20.0, substeps=5):\n        self.sdf_truncate = sdf_truncate\n        self.substeps = substeps  # sub-voxel sampling along normal\n\n    def transform(self, mask):\n        mask = mask.astype(bool)\n        self.mask_shape = mask.shape\n        sdf = self._compute_sdf(mask)\n        normals = self._normals_from_sdf(sdf)\n        thickness = self._compute_thickness(mask)\n        return {\"sdf\": sdf, \"normals\": normals, \"thickness\": thickness}\n\n    def inverse_transform(self, pred):\n        D, H, W = self.mask_shape\n        sdf = pred[\"sdf\"]\n        normals = pred[\"normals\"]\n        thickness = pred[\"thickness\"]\n\n        # --- Center surface coordinates ---\n        center_mask = sdf <= 0\n        coords = np.argwhere(center_mask)\n        if coords.size == 0:\n            return np.zeros((D, H, W), dtype=np.uint8)\n\n        dz = normals[0][center_mask]\n        dy = normals[1][center_mask]\n        dx = normals[2][center_mask]\n        t_half = thickness[center_mask] / 2.0\n\n        # --- Sub-voxel steps along normals ---\n        num_centers = coords.shape[0]\n        steps = np.linspace(-1, 1, self.substeps)[None, :] * t_half[:, None]\n\n        z_coords = coords[:, 0:1] + dz[:, None] * steps\n        y_coords = coords[:, 1:2] + dy[:, None] * steps\n        x_coords = coords[:, 2:3] + dx[:, None] * steps\n\n        # --- Round and clip voxel indices ---\n        z_idx = np.clip(np.round(z_coords).astype(int), 0, D-1).ravel()\n        y_idx = np.clip(np.round(y_coords).astype(int), 0, H-1).ravel()\n        x_idx = np.clip(np.round(x_coords).astype(int), 0, W-1).ravel()\n\n        # --- Sparse accumulation ---\n        flat_idx = z_idx*H*W + y_idx*W + x_idx\n        ones = np.ones_like(flat_idx, dtype=np.uint16)\n        recon_sparse = coo_matrix((ones, (flat_idx, np.zeros_like(flat_idx))),\n                                  shape=(D*H*W, 1))\n\n        # Convert to dense mask\n        recon_flat = (recon_sparse.toarray().ravel() > 0).astype(np.uint8)\n        recon = recon_flat.reshape(D, H, W)\n\n        return recon\n\n    def _compute_sdf(self, mask):\n        dist_out = distance_transform_edt(~mask)\n        dist_in = distance_transform_edt(mask)\n        sdf = dist_out - dist_in\n        return np.clip(sdf, -self.sdf_truncate, self.sdf_truncate).astype(np.float32)\n\n    def _normals_from_sdf(self, sdf):\n        dz, dy, dx = np.gradient(sdf)\n        grad = np.stack([dz, dy, dx], axis=0)\n        norm = np.linalg.norm(grad, axis=0, keepdims=True) + 1e-8\n        return (grad / norm).astype(np.float32)\n\n    def _compute_thickness(self, mask):\n        dist_in = distance_transform_edt(mask)\n        dist_out = distance_transform_edt(~mask)\n        thickness = dist_in + dist_out\n        return thickness.astype(np.float32)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-20T03:13:55.297532Z","iopub.execute_input":"2025-12-20T03:13:55.297868Z","iopub.status.idle":"2025-12-20T03:13:55.314527Z","shell.execute_reply.started":"2025-12-20T03:13:55.297840Z","shell.execute_reply":"2025-12-20T03:13:55.313535Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"data_path = Path('/kaggle/input/vesuvius-challenge-surface-detection/train_labels/1407735.tif')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-20T03:06:24.122745Z","iopub.execute_input":"2025-12-20T03:06:24.123077Z","iopub.status.idle":"2025-12-20T03:06:24.128825Z","shell.execute_reply.started":"2025-12-20T03:06:24.123048Z","shell.execute_reply":"2025-12-20T03:06:24.127822Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"mask = load_volume(data_path)\nmask = mask * (mask != 2)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-20T03:06:24.282617Z","iopub.execute_input":"2025-12-20T03:06:24.282976Z","iopub.status.idle":"2025-12-20T03:06:24.787317Z","shell.execute_reply.started":"2025-12-20T03:06:24.282944Z","shell.execute_reply":"2025-12-20T03:06:24.786297Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"decomp = ExactScrollDecomposition()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-20T03:13:57.712628Z","iopub.execute_input":"2025-12-20T03:13:57.712975Z","iopub.status.idle":"2025-12-20T03:13:57.717605Z","shell.execute_reply.started":"2025-12-20T03:13:57.712945Z","shell.execute_reply":"2025-12-20T03:13:57.716416Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pred = decomp.transform(mask)\nrecon = decomp.inverse_transform(pred)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-20T03:13:57.857689Z","iopub.execute_input":"2025-12-20T03:13:57.858078Z","iopub.status.idle":"2025-12-20T03:14:45.631600Z","shell.execute_reply.started":"2025-12-20T03:13:57.858049Z","shell.execute_reply":"2025-12-20T03:14:45.630304Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.imshow(recon[100])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-20T03:14:52.771145Z","iopub.execute_input":"2025-12-20T03:14:52.771491Z","iopub.status.idle":"2025-12-20T03:14:52.962668Z","shell.execute_reply.started":"2025-12-20T03:14:52.771460Z","shell.execute_reply":"2025-12-20T03:14:52.961757Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"fig, ax = plt.subplots(ncols = 5, figsize = (24, 24))\nax[0].imshow(pred['thickness'][100])\nax[0].set_title('Thinkness')\nax[1].imshow(pred['sdf'][100])\nax[1].set_title('sign distance function')\nax[2].imshow(pred['normals'][0, 100])\nax[2].set_title('normals dx')\nax[3].imshow(pred['normals'][1, 100])\nax[3].set_title('normals dy')\nax[4].imshow(pred['normals'][2, 100])\nax[4].set_title('normals dz')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-20T03:14:54.287419Z","iopub.execute_input":"2025-12-20T03:14:54.287967Z","iopub.status.idle":"2025-12-20T03:14:55.372873Z","shell.execute_reply.started":"2025-12-20T03:14:54.287932Z","shell.execute_reply":"2025-12-20T03:14:55.371760Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"Mask voxels:\", mask.sum())\nprint(\"Recon voxels:\", recon.sum())\nprint(\"Reconstruction error:\", np.sum(recon != mask))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-20T03:14:55.374538Z","iopub.execute_input":"2025-12-20T03:14:55.374877Z","iopub.status.idle":"2025-12-20T03:14:55.455658Z","shell.execute_reply.started":"2025-12-20T03:14:55.374848Z","shell.execute_reply":"2025-12-20T03:14:55.454562Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Decomposition B","metadata":{}},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}