{"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":"none","dataSources":[{"sourceId":117682,"databundleVersionId":15062069,"sourceType":"competition"},{"sourceId":14261875,"sourceType":"datasetVersion","datasetId":8670484},{"sourceId":14364531,"sourceType":"datasetVersion","datasetId":9172693},{"sourceId":14407441,"sourceType":"datasetVersion","datasetId":9181226},{"sourceId":285233081,"sourceType":"kernelVersion"},{"sourceId":286093322,"sourceType":"kernelVersion"},{"sourceId":289436856,"sourceType":"kernelVersion"},{"sourceId":290272223,"sourceType":"kernelVersion"}],"dockerImageVersionId":31192,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"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\nimport matplotlib.pyplot as plt\nimport os\nfrom PIL import Image\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\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:20.065691Z","iopub.execute_input":"2026-01-06T15:33:20.065999Z","iopub.status.idle":"2026-01-06T15:33:20.432939Z","shell.execute_reply.started":"2026-01-06T15:33:20.065967Z","shell.execute_reply":"2026-01-06T15:33:20.431951Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# ADDED NEW FEATURES","metadata":{}},{"cell_type":"markdown","source":"1. Find Normal to flow of sheets.\n2. Find Seeds for Geodesic Veronoi or other similar methods.\n3. Geodesic Voronoi and a similar Euclidian method \n4. Estimating Number of Sheets\n5. EDF Visualization","metadata":{}},{"cell_type":"markdown","source":"# IMPORTANT NOTE :\nThe \"prediction_volume.tif\" will automatically be saved when you run compute button in the mini-app, and prediction will be saved based on the settings selected. Which then you are supposed to use throughout the notebook to try various post processing steps, also update the paths for your probability predictions, I am too lazy to provide some samples and currently can't give my predictions as I still am competing.","metadata":{}},{"cell_type":"markdown","source":"# Summary\n- This is a basic notebook, to  observe probability outputs of your models via various standard techniques, like quantiles, std of predictions, Histogram.\n- You can compute your competition metric score.\n- It also has an interactive widget, to find best probability thresholds for inference, if you are not using the given topological refinements, simply use thigh=tlow in the slider to find the best threshold for you.\n- I made this to fine tune some hyperparamters for post processing.\n- Currently two train-cases are hardcoded, 1004283650 and 1006462223. 1006462223 is quite tough.","metadata":{}},{"cell_type":"code","source":"def load_tif_volume(path):\n    \"\"\"\n    Memory optimized loader. Pre-allocates array to avoid \n    memory spike during stacking.\n    \"\"\"\n    \n    img = Image.open(path)\n    W, H = img.size\n    Z = getattr(img, 'n_frames', 1)\n    \n    print(f\"Loading {path} ({Z}x{H}x{W})...\")\n    \n    # Pre-allocate buffer (Standard RAM)\n    vol = np.zeros((Z, H, W), dtype=np.uint8)\n    \n    for i in range(Z):\n        try:\n            img.seek(i)\n            # Write directly to buffer\n            vol[i] = np.array(img, dtype=np.uint8)\n        except EOFError:\n            break\n            \n    return vol\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:30.510114Z","iopub.execute_input":"2026-01-06T15:33:30.511201Z","iopub.status.idle":"2026-01-06T15:33:30.518501Z","shell.execute_reply.started":"2026-01-06T15:33:30.511170Z","shell.execute_reply":"2026-01-06T15:33:30.517092Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"**Paths of your predictions**\n- hardcoded for .npy can be changed very simply","metadata":{}},{"cell_type":"code","source":"pred1p = \"/kaggle/input/50-epochs-occlusion/probabilities_npy/1004283650_prob.npy\"\npred2p = \"/kaggle/input/50-epochs-occlusion/probabilities_npy/1006462223_prob.npy\"\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:13.916061Z","iopub.execute_input":"2026-01-06T15:33:13.916416Z","iopub.status.idle":"2026-01-06T15:33:13.922178Z","shell.execute_reply.started":"2026-01-06T15:33:13.916384Z","shell.execute_reply":"2026-01-06T15:33:13.920988Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pred3p = \"/kaggle/input/50-epochs-occlusion/probabilities_npy/1182039038_prob.npy\"\n# pred4p = \"/kaggle/input/250-epochs-highres/nnUNet_results/Dataset100_VesuviusSurface/nnUNetTrainer_500epochs__nnUNetResEncUNetMPlans__3d_fullres/fold_all/validation/1370134939.tif\"","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:14.154930Z","iopub.execute_input":"2026-01-06T15:33:14.155869Z","iopub.status.idle":"2026-01-06T15:33:14.160416Z","shell.execute_reply.started":"2026-01-06T15:33:14.155784Z","shell.execute_reply":"2026-01-06T15:33:14.159354Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"gt1p = \"/kaggle/input/vesuvius-challenge-surface-detection/train_labels/1004283650.tif\"\ngt2p = \"/kaggle/input/vesuvius-challenge-surface-detection/train_labels/1006462223.tif\"","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:14.517942Z","iopub.execute_input":"2026-01-06T15:33:14.518317Z","iopub.status.idle":"2026-01-06T15:33:14.523342Z","shell.execute_reply.started":"2026-01-06T15:33:14.518266Z","shell.execute_reply":"2026-01-06T15:33:14.522326Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"gt3p = \"/kaggle/input/vesuvius-challenge-surface-detection/train_labels/1182039038.tif\"\ngt4p = \"/kaggle/input/vesuvius-challenge-surface-detection/train_labels/1370134939.tif\"","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:14.857936Z","iopub.execute_input":"2026-01-06T15:33:14.858337Z","iopub.status.idle":"2026-01-06T15:33:14.863730Z","shell.execute_reply.started":"2026-01-06T15:33:14.858275Z","shell.execute_reply":"2026-01-06T15:33:14.862554Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Loading files\n- vol1 and vol2 are our probability volumes\n- gt1 and gt2 are ground truths/ labels\n- inp_1 and inp_2 are the inputs","metadata":{}},{"cell_type":"code","source":"vol1= np.load(pred1p)\nvol2= np.load(pred2p)\n# vol1 = load_tif_volume(pred1p)\n# vol2 = load_tif_volume(pred2p)\ngt1 =load_tif_volume(gt1p)\ngt2 =load_tif_volume(gt2p)\ninp_1 =load_tif_volume(\"/kaggle/input/vesuvius-challenge-surface-detection/train_images/1004283650.tif\")\ninp_2 =load_tif_volume(\"/kaggle/input/vesuvius-challenge-surface-detection/train_images/1006462223.tif\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:33.543034Z","iopub.execute_input":"2026-01-06T15:33:33.543399Z","iopub.status.idle":"2026-01-06T15:33:39.752605Z","shell.execute_reply.started":"2026-01-06T15:33:33.543374Z","shell.execute_reply":"2026-01-06T15:33:39.751580Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"vol3= np.load(pred3p)\n# vol4= np.load(pred4p)\n# vol3 = load_tif_volume(pred3p)\n# vol4 = load_tif_volume(pred4p)\n\ngt3 =load_tif_volume(gt3p)\n# gt4 =load_tif_volume(gt4p)\ninp_3 =load_tif_volume(\"/kaggle/input/vesuvius-challenge-surface-detection/train_images/1182039038.tif\")\n# inp_4 =load_tif_volume(\"/kaggle/input/vesuvius-challenge-surface-detection/train_images/975031774.tif\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:39.754383Z","iopub.execute_input":"2026-01-06T15:33:39.754738Z","iopub.status.idle":"2026-01-06T15:33:43.003568Z","shell.execute_reply.started":"2026-01-06T15:33:39.754707Z","shell.execute_reply":"2026-01-06T15:33:43.002324Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# vol1, vol2 = vol1.reshape(320,320,320,2), vol2.reshape(320,320,320,2)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:43.004585Z","iopub.execute_input":"2026-01-06T15:33:43.004822Z","iopub.status.idle":"2026-01-06T15:33:43.009822Z","shell.execute_reply.started":"2026-01-06T15:33:43.004799Z","shell.execute_reply":"2026-01-06T15:33:43.008810Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# vol1 = vol1[:,:,:,1:2].reshape(320,320,320)\n# vol2 = vol2[:,:,:,1:2].reshape(320,320,320)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:43.012098Z","iopub.execute_input":"2026-01-06T15:33:43.012517Z","iopub.status.idle":"2026-01-06T15:33:43.037542Z","shell.execute_reply.started":"2026-01-06T15:33:43.012487Z","shell.execute_reply":"2026-01-06T15:33:43.036386Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"vol1.shape","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:43.038489Z","iopub.execute_input":"2026-01-06T15:33:43.038736Z","iopub.status.idle":"2026-01-06T15:33:43.060210Z","shell.execute_reply.started":"2026-01-06T15:33:43.038715Z","shell.execute_reply":"2026-01-06T15:33:43.059132Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from ipywidgets import interact\nimport ipywidgets as widgets\nimport panel as pn","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:43.061229Z","iopub.execute_input":"2026-01-06T15:33:43.061597Z","iopub.status.idle":"2026-01-06T15:33:44.750236Z","shell.execute_reply.started":"2026-01-06T15:33:43.061574Z","shell.execute_reply":"2026-01-06T15:33:44.748986Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Simple Interactive window\nto visualize input, pred, and the ground truth side by side","metadata":{}},{"cell_type":"code","source":"def explore_side_by_side(gt, pred,inp):\n    \"\"\"\n    Creates an interactive slider to scroll through the Z-axis.\n    \"\"\"\n    # 1. Determine the maximum depth based on your array shape (H, W, Depth)\n    max_z = gt.shape[2] - 1\n    \n    # 2. Define the callback function that runs every time the slider moves\n    def plot_slice(z):\n        fig, axs = plt.subplots(1, 3, figsize=(12, 4))\n        \n        # Ground Truth\n        axs[0].imshow(gt[z,:,:], cmap='gray')\n        axs[0].set_title(f\"GT (z={z})\")\n        \n        # Prediction\n        axs[1].imshow(pred[z,:,:], cmap='gray')\n        axs[1].set_title(f\"Pred (z={z})\")\n\n\n        axs[2].imshow(inp[z,:,:],cmap = \"gray\")\n        axs[2].set_title(\"input (z={z})\")\n        \n        \n        for ax in axs: ax.axis('off')\n        plt.show()\n\n    # 3. Create the widget\n    interact(plot_slice, \n             z=widgets.IntSlider(min=0, max=max_z, step=1, value=max_z//2, description='Slice Z:'))\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:44.751349Z","iopub.execute_input":"2026-01-06T15:33:44.751692Z","iopub.status.idle":"2026-01-06T15:33:44.760607Z","shell.execute_reply.started":"2026-01-06T15:33:44.751668Z","shell.execute_reply":"2026-01-06T15:33:44.759398Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"explore_side_by_side(gt1,vol1,inp_1) ##1004283650\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:44.761780Z","iopub.execute_input":"2026-01-06T15:33:44.762162Z","iopub.status.idle":"2026-01-06T15:33:45.200130Z","shell.execute_reply.started":"2026-01-06T15:33:44.762128Z","shell.execute_reply":"2026-01-06T15:33:45.199091Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"explore_side_by_side(gt2,vol2,inp_2) ##1006462223\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:45.201361Z","iopub.execute_input":"2026-01-06T15:33:45.201763Z","iopub.status.idle":"2026-01-06T15:33:45.567241Z","shell.execute_reply.started":"2026-01-06T15:33:45.201732Z","shell.execute_reply":"2026-01-06T15:33:45.566163Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"explore_side_by_side(gt3,vol3,inp_3) \n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:45.569581Z","iopub.execute_input":"2026-01-06T15:33:45.569899Z","iopub.status.idle":"2026-01-06T15:33:45.901159Z","shell.execute_reply.started":"2026-01-06T15:33:45.569876Z","shell.execute_reply":"2026-01-06T15:33:45.900081Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# explore_side_by_side(gt4,vol4,inp_4) \n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:45.902214Z","iopub.execute_input":"2026-01-06T15:33:45.902608Z","iopub.status.idle":"2026-01-06T15:33:45.908375Z","shell.execute_reply.started":"2026-01-06T15:33:45.902581Z","shell.execute_reply":"2026-01-06T15:33:45.906757Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!pip install connected-components-3d \nimport cc3d\n\n\nimport numpy as np\n\nfrom scipy.spatial.distance import cdist\nfrom scipy.optimize import linear_sum_assignment\nfrom scipy.ndimage import convolve\nfrom scipy.ndimage import gaussian_filter\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:45.909673Z","iopub.execute_input":"2026-01-06T15:33:45.910136Z","iopub.status.idle":"2026-01-06T15:33:50.669489Z","shell.execute_reply.started":"2026-01-06T15:33:45.910101Z","shell.execute_reply":"2026-01-06T15:33:50.668172Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Widget to fine tune\n- This is the widget to fine tune a few hyperparameters for some basic topological post processing","metadata":{}},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Quantiles and histograms","metadata":{}},{"cell_type":"code","source":"fg_probs = vol1[gt1 == 1].astype(np.float64)\nmean_fg = fg_probs.mean()\nstd_fg  = fg_probs.std()\nq25_fg = np.percentile(fg_probs, 25)\nq50_fg = np.percentile(fg_probs, 50)\nq75_fg = np.percentile(fg_probs, 75)\nprint(std_fg,q25_fg,q50_fg,q75_fg)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:54.170819Z","iopub.execute_input":"2026-01-06T15:33:54.171147Z","iopub.status.idle":"2026-01-06T15:33:54.468699Z","shell.execute_reply.started":"2026-01-06T15:33:54.171119Z","shell.execute_reply":"2026-01-06T15:33:54.467547Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"fg_probs = vol1[gt1 == 0].astype(np.float64)\nmean_fg = fg_probs.mean()\nstd_fg  = fg_probs.std()\nq25_fg = np.percentile(fg_probs, 25)\nq50_fg = np.percentile(fg_probs, 50)\nq75_fg = np.percentile(fg_probs, 90)\nprint(std_fg,q25_fg,q50_fg,q75_fg)\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:54.470390Z","iopub.execute_input":"2026-01-06T15:33:54.471153Z","iopub.status.idle":"2026-01-06T15:33:55.028115Z","shell.execute_reply.started":"2026-01-06T15:33:54.471126Z","shell.execute_reply":"2026-01-06T15:33:55.027136Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"fg_probs = vol2[gt2 == 1].astype(np.float64)\nmean_fg = fg_probs.mean()\nstd_fg  = fg_probs.std()\nq25_fg = np.percentile(fg_probs, 25)\nq50_fg = np.percentile(fg_probs, 50)\nq75_fg = np.percentile(fg_probs, 75)\nprint(std_fg,q25_fg,q50_fg,q75_fg)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:55.029627Z","iopub.execute_input":"2026-01-06T15:33:55.030124Z","iopub.status.idle":"2026-01-06T15:33:55.417843Z","shell.execute_reply.started":"2026-01-06T15:33:55.030094Z","shell.execute_reply":"2026-01-06T15:33:55.416647Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"fg_probs = vol2[gt2 == 0].astype(np.float64)\nmean_fg = fg_probs.mean()\nstd_fg  = fg_probs.std()\nq25_fg = np.percentile(fg_probs, 25)\nq50_fg = np.percentile(fg_probs, 50)\nq75_fg = np.percentile(fg_probs, 75)\nprint(std_fg,q25_fg,q50_fg,q75_fg)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:55.419499Z","iopub.execute_input":"2026-01-06T15:33:55.419858Z","iopub.status.idle":"2026-01-06T15:33:57.035282Z","shell.execute_reply.started":"2026-01-06T15:33:55.419829Z","shell.execute_reply":"2026-01-06T15:33:57.034258Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"fg_probs = vol3[gt3 == 1].astype(np.float64)\nmean_fg = fg_probs.mean()\nstd_fg  = fg_probs.std()\nq25_fg = np.percentile(fg_probs, 25)\nq50_fg = np.percentile(fg_probs, 50)\nq75_fg = np.percentile(fg_probs, 75)\nprint(std_fg,q25_fg,q50_fg,q75_fg)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:33:59.368526Z","iopub.execute_input":"2026-01-06T15:33:59.368841Z","iopub.status.idle":"2026-01-06T15:33:59.507012Z","shell.execute_reply.started":"2026-01-06T15:33:59.368818Z","shell.execute_reply":"2026-01-06T15:33:59.506065Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"fg_probs = vol3[gt3 == 0].astype(np.float64)\nmean_fg = fg_probs.mean()\nstd_fg  = fg_probs.std()\nq25_fg = np.percentile(fg_probs, 25)\nq50_fg = np.percentile(fg_probs, 50)\nq75_fg = np.percentile(fg_probs, 75)\nprint(std_fg,q25_fg,q50_fg,q75_fg)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T12:54:17.542124Z","iopub.execute_input":"2026-01-06T12:54:17.542516Z","iopub.status.idle":"2026-01-06T12:54:17.971518Z","shell.execute_reply.started":"2026-01-06T12:54:17.542485Z","shell.execute_reply":"2026-01-06T12:54:17.970269Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# fg_probs = vol4[gt4 == 1].astype(np.float64)\n# mean_fg = fg_probs.mean()\n# std_fg  = fg_probs.std()\n# q25_fg = np.percentile(fg_probs, 25)\n# q50_fg = np.percentile(fg_probs, 50)\n# q75_fg = np.percentile(fg_probs, 75)\n# print(std_fg,q25_fg,q50_fg,q75_fg)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T12:54:17.972373Z","iopub.execute_input":"2026-01-06T12:54:17.972683Z","iopub.status.idle":"2026-01-06T12:54:17.979138Z","shell.execute_reply.started":"2026-01-06T12:54:17.972662Z","shell.execute_reply":"2026-01-06T12:54:17.978121Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# fg_probs = vol4[gt4 == 0].astype(np.float64)\n# mean_fg = fg_probs.mean()\n# std_fg  = fg_probs.std()\n# q25_fg = np.percentile(fg_probs, 25)\n# q50_fg = np.percentile(fg_probs, 50)\n# q75_fg = np.percentile(fg_probs, 75)\n# print(std_fg,q25_fg,q50_fg,q75_fg)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T12:54:17.980194Z","iopub.execute_input":"2026-01-06T12:54:17.980559Z","iopub.status.idle":"2026-01-06T12:54:18.001615Z","shell.execute_reply.started":"2026-01-06T12:54:17.980528Z","shell.execute_reply":"2026-01-06T12:54:18.000336Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# print(vol1.shape,vol1.max(),vol1.min())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T12:54:18.002934Z","iopub.execute_input":"2026-01-06T12:54:18.003338Z","iopub.status.idle":"2026-01-06T12:54:18.023381Z","shell.execute_reply.started":"2026-01-06T12:54:18.003307Z","shell.execute_reply":"2026-01-06T12:54:18.022311Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\n\ndef plot_prob_hist(probs, title, bins=200):\n    plt.figure(figsize=(6,4))\n    plt.hist(probs, bins=bins, density=False)\n    plt.yscale('log')\n    plt.title(title)\n    plt.xlabel(\"Probability\")\n    plt.ylabel(\"Count (log)\")\n    plt.grid(alpha=0.3)\n    plt.show()\n\n# ---------- VOL1 ----------\nfg1 = vol1[gt1 == 1].astype(np.float64)\nbg1 = vol1[gt1 == 0].astype(np.float64)\n\nplot_prob_hist(fg1, \"VOL1 – FG Probabilities\")\nplot_prob_hist(bg1, \"VOL1 – BG Probabilities\")\n\n\n# ---------- VOL2 ----------\nfg2 = vol2[gt2 == 1].astype(np.float64)\nbg2 = vol2[gt2 == 0].astype(np.float64)\n\nplot_prob_hist(fg2, \"VOL2 – FG Probabilities\")\nplot_prob_hist(bg2, \"VOL2 – BG Probabilities\")\n\n# ---------- VOL3 ----------\nfg3 = vol3[gt3 == 1].astype(np.float64)\nbg3 = vol3[gt3 == 0].astype(np.float64)\n\nplot_prob_hist(fg3, \"VOL3 – FG Probabilities\")\nplot_prob_hist(bg3, \"VOL3 – BG Probabilities\")\n\n# # ---------- VOL4----------\n# fg4 = vol4[gt4 == 1].astype(np.float64)\n# bg4 = vol4[gt4 == 0].astype(np.float64)\n\n# plot_prob_hist(fg4, \"VOL4 – FG Probabilities\")\n# plot_prob_hist(bg4, \"VOL4 – BG Probabilities\")\n\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T12:54:18.030859Z","iopub.execute_input":"2026-01-06T12:54:18.031441Z","iopub.status.idle":"2026-01-06T12:54:24.021729Z","shell.execute_reply.started":"2026-01-06T12:54:18.031400Z","shell.execute_reply":"2026-01-06T12:54:24.020281Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# import numpy as np\n# import matplotlib.pyplot as plt\n# from skimage.filters import meijering\n# from skimage import exposure\n\n# # 1. Select a slice\n# z_idx = 30\n# prob_slice = vol2[z_idx, :, :]\n# gt_slice = gt2[z_idx, :, :]\n\n# # 2. Parameters\n# # sigmas: range(1, 4) usually works well for ink lines (1-3px radius)\n# sigmas = range(1, 4)\n# pad_size = 20\n\n# # 3. Apply Meijering with Padding (Crucial to avoid edge artifacts)\n# # Reflect padding prevents the border from looking like a \"cliff\"\n# padded_slice = np.pad(prob_slice, pad_size, mode='reflect')\n# ridge_map_padded = meijering(padded_slice, sigmas=sigmas, black_ridges=False)\n# ridge_map = ridge_map_padded[pad_size:-pad_size, pad_size:-pad_size]\n\n# # 4. Normalize and Threshold\n# # Normalizing helps make the threshold (0.05 - 0.2) consistent across slices\n# if ridge_map.max() > 0:\n#     ridge_norm = ridge_map / ridge_map.max()\n# else:\n#     ridge_norm = ridge_map\n\n# final_mask = ridge_norm > 0.15  # Tune this threshold\n\n# # 5. Visualization\n# fig, ax = plt.subplots(1, 4, figsize=(20, 5))\n\n# # A. Original Threshold (likely merged)\n# ax[0].imshow(prob_slice > 0.8, cmap='gray')\n# ax[0].set_title(f\"Original Threshold (>0.8)\")\n# ax[0].axis('off')\n\n# # B. Meijering Output (Heatmap)\n# im1 = ax[1].imshow(ridge_norm, cmap='inferno')\n# ax[1].set_title(\"Meijering (Tubeness)\")\n# ax[1].axis('off')\n# plt.colorbar(im1, ax=ax[1], fraction=0.046, pad=0.04)\n\n# # C. Meijering Thresholded\n# ax[2].imshow(final_mask, cmap='gray')\n# ax[2].set_title(f\"Meijering Mask (>0.15)\")\n# ax[2].axis('off')\n\n# # D. Ground Truth\n# ax[3].imshow(gt_slice, cmap='gray')\n# ax[3].set_title(\"Ground Truth\")\n# ax[3].axis('off')\n\n# plt.tight_layout()\n# plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T12:54:24.023143Z","iopub.execute_input":"2026-01-06T12:54:24.023571Z","iopub.status.idle":"2026-01-06T12:54:24.029646Z","shell.execute_reply.started":"2026-01-06T12:54:24.023539Z","shell.execute_reply":"2026-01-06T12:54:24.028312Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# The widget","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\nimport ipywidgets as widgets\nfrom PIL import Image\nfrom IPython.display import display, clear_output\n\nfrom scipy.ndimage import (\n    label,\n    binary_closing,\n    generate_binary_structure,\n    binary_propagation,\n)\nfrom skimage.morphology import remove_small_objects\nfrom skimage.filters import meijering\nfrom skimage import exposure\n\n%matplotlib inline\n\n# ==========================================\n# 1. HELPERS\n# ==========================================\n\ndef apply_hysteresis_3d(prob_vol, t_low, t_high):\n    strong = prob_vol >= t_high\n    weak   = prob_vol >= t_low\n    struct3d = generate_binary_structure(3, 3)\n    return binary_propagation(strong, mask=weak, structure=struct3d)\n\ndef apply_hysteresis_2d(prob_slice, t_low, t_high):\n    strong = prob_slice >= t_high\n    weak   = prob_slice >= t_low\n    struct2d = generate_binary_structure(2, 2)\n    return binary_propagation(strong, mask=weak, structure=struct2d)\n\ndef build_anisotropic_struct(z_radius: int, xy_radius: int):\n    z, r = z_radius, xy_radius\n    if z == 0 and r == 0: return None\n    if z == 0 and r > 0:\n        size = 2 * r + 1\n        struct = np.zeros((1, size, size), dtype=bool)\n        y, x = np.ogrid[-r:r+1, -r:r+1]\n        mask = x**2 + y**2 <= r**2\n        struct[0] = mask\n        return struct\n    depth = 2 * z + 1\n    size = 2 * r + 1\n    struct = np.zeros((depth, size, size), dtype=bool)\n    y, x = np.ogrid[-r:r+1, -r:r+1]\n    mask = x**2 + y**2 <= r**2\n    for dz in range(depth):\n        struct[dz] = mask\n    return struct\n\ndef apply_anisotropic_closing_3d(mask3d, z_radius: int, xy_radius: int):\n    struct = build_anisotropic_struct(z_radius, xy_radius)\n    if struct is None: return mask3d\n    return binary_closing(mask3d, structure=struct)\n\ndef apply_meijering_slice_safe(img_slice, sigma_range, filter_padding=20):\n    padded = np.pad(img_slice, filter_padding, mode='reflect')\n    filtered = meijering(padded, sigmas=sigma_range, black_ridges=False)\n    if filter_padding > 0:\n        return filtered[filter_padding:-filter_padding, filter_padding:-filter_padding]\n    return filtered\n\n# ==========================================\n# 2. INTERACTIVE TOPOLOGY TOOL (FIXED)\n# ==========================================\n\nclass KaggleTopologyOptimizer:\n    def __init__(self, prob_volume, gt_volume, spacing=(1, 1, 1)):\n        self.vol = prob_volume.astype(np.float32)\n        self.gt  = gt_volume\n        self.spacing = spacing\n\n        self.out_plot = widgets.Output()\n        self.out_log  = widgets.Output()\n\n        style = {'description_width': 'initial'}\n        Z = self.vol.shape[0]\n\n        # --- SLIDERS (continuous_update=False to prevent lag) ---\n        self.s_z      = widgets.IntSlider(min=0, max=Z - 1, value=Z // 2, description='Z Slice:', style=style, continuous_update=False)\n        \n        # Hysteresis / Simple Threshold\n        self.s_tlow   = widgets.FloatSlider(min=0.0, max=1.0, value=0.4, step=0.01, description='Hyst Low:',  style=style, continuous_update=False)\n        self.s_thigh  = widgets.FloatSlider(min=0.0, max=1.0, value=0.8, step=0.01, description='Hyst High:', style=style, continuous_update=False)\n        \n        # Meijering\n        self.chk_meij   = widgets.Checkbox(value=False, description='Enable Meijering (Center Only)')\n        self.s_margin   = widgets.IntSlider(min=0, max=200, value=50, description='Central Margin:', style=style, continuous_update=False)\n        self.s_mthresh  = widgets.FloatSlider(min=0.0, max=1.0, value=0.1, step=0.01, description='Tube Thresh:', style=style, continuous_update=False)\n        \n        # Meijering Tuning\n        self.s_msig_min = widgets.IntSlider(min=1, max=5, value=1, description='Sigma Min:', style=style, continuous_update=False)\n        self.s_msig_max = widgets.IntSlider(min=2, max=8, value=4, description='Sigma Max:', style=style, continuous_update=False)\n        self.s_mpad     = widgets.IntSlider(min=0, max=50, value=20, description='Filter Pad:', style=style, continuous_update=False)\n\n        # Morphology\n        self.s_zrad   = widgets.IntSlider(min=0, max=5, value=1, description='Close Z Rad:',  style=style, continuous_update=False)\n        self.s_xyrad  = widgets.IntSlider(min=0, max=5, value=0, description='Close XY Rad:', style=style, continuous_update=False)\n        self.s_minsize = widgets.IntSlider(min=0, max=2000, value=100, step=50, description='Dust Min:', style=style, continuous_update=False)\n\n        self.btn_compute = widgets.Button(description='Compute Full 3D Metrics', button_style='info')\n        self.btn_compute.on_click(self.compute_3d_metrics)\n\n        # Explicitly observe all widgets\n        widgets_list = [\n            self.s_z, self.s_tlow, self.s_thigh, \n            self.s_zrad, self.s_xyrad, self.s_minsize, \n            self.chk_meij, self.s_msig_min, self.s_msig_max, \n            self.s_mpad, self.s_mthresh, self.s_margin\n        ]\n        \n        for w in widgets_list:\n            w.observe(self.update_view, names='value')\n\n        ui = widgets.VBox([\n            self.s_z,\n            widgets.HTML(\"<b>1. Base Thresholds (Used for Edges & Hysteresis)</b>\"),\n            widgets.HBox([self.s_tlow, self.s_thigh]),\n            widgets.HTML(\"<b>2. Meijering Filter (Applied to Central Patch)</b>\"),\n            widgets.HBox([self.chk_meij, self.s_margin]),\n            widgets.HBox([self.s_mthresh, self.s_msig_min, self.s_msig_max, self.s_mpad]),\n            widgets.HTML(\"<i>Note: Outside the 'Central Margin', we use a simple threshold (Hyst High).</i>\"),\n            widgets.HTML(\"<b>3. Morphology & Cleaning</b>\"),\n            widgets.HBox([self.s_zrad, self.s_xyrad, self.s_minsize]),\n            widgets.HTML(\"<hr>\"),\n            self.btn_compute\n        ])\n\n        display(ui, self.out_plot, self.out_log)\n        self.update_view(None)\n\n    def get_processed_slice(self, z):\n        prob_slice = self.vol[z]\n        h, w = prob_slice.shape\n        margin = self.s_margin.value\n        \n        # --- STRATEGY SELECTION ---\n        if self.chk_meij.value:\n            # 1. Define Central Patch Mask (Fixed Zero Margin Logic)\n            center_mask = np.zeros_like(prob_slice, dtype=bool)\n            \n            if margin == 0:\n                center_mask[:] = True\n            elif margin < h//2 and margin < w//2:\n                center_mask[margin:-margin, margin:-margin] = True\n            else:\n                center_mask[:] = True \n            \n            # 2. Compute Meijering (Inner)\n            s_min = self.s_msig_min.value\n            s_max = max(s_min + 1, self.s_msig_max.value)\n            f_pad = self.s_mpad.value\n            \n            tubeness = apply_meijering_slice_safe(prob_slice, range(s_min, s_max), filter_padding=f_pad)\n            if tubeness.max() > 0: tubeness /= tubeness.max()\n            \n            mask_inner = tubeness > self.s_mthresh.value\n            \n            # 3. Compute Simple Threshold (Outer)\n            mask_outer = prob_slice > self.s_thigh.value\n            \n            # 4. Combine: Inner logic in center, Outer logic at borders\n            mask2d = np.where(center_mask, mask_inner, mask_outer)\n            \n        else:\n            # Fallback to standard Hysteresis everywhere\n            t_low, t_high = self.s_tlow.value, self.s_thigh.value\n            mask2d = apply_hysteresis_2d(prob_slice, t_low, t_high)\n\n        # --- MORPHOLOGY ---\n        xy_r = self.s_xyrad.value\n        if xy_r > 0:\n            struct = build_anisotropic_struct(0, xy_r)\n            from scipy.ndimage import binary_closing as closing_2d\n            struct_2d = struct[0] \n            mask2d = closing_2d(mask2d, structure=struct_2d)\n\n        # --- DUST REMOVAL ---\n        if self.s_minsize.value > 0:\n            mask2d = remove_small_objects(mask2d.astype(bool), min_size=self.s_minsize.value)\n\n        return mask2d\n\n    def update_view(self, change):\n        z = self.s_z.value\n        m = self.s_margin.value\n        \n        pred_mask = self.get_processed_slice(z)\n        gt_slice  = self.gt[z]\n\n        gt_viz = np.zeros((*gt_slice.shape, 3), dtype=np.float32)\n        gt_viz[gt_slice == 1] = [1, 1, 1]  # FG\n        gt_viz[gt_slice == 2] = [1, 0, 0]  # Ignore\n\n        with self.out_plot:\n            clear_output(wait=True)\n            \n            fig, axes = plt.subplots(1, 2, figsize=(12, 6))\n\n            axes[0].imshow(gt_viz)\n            axes[0].set_title(f\"GT (z={z})\\nWhite=FG, Red=Ignore\", fontsize=13)\n            axes[0].axis('off')\n            \n            # Draw Prediction Image\n            axes[1].imshow(pred_mask, cmap=\"gray\")\n            \n            # --- FIX: Draw Margin Box ALWAY ---\n            # We draw this regardless of self.chk_meij.value so you can set it up first\n            h, w = pred_mask.shape\n            \n            if m > 0 and m < h//2 and m < w//2:\n                # x, y, width, height\n                rect = plt.Rectangle((m, m), w - 2*m, h - 2*m, \n                                     linewidth=2, edgecolor='cyan', facecolor='none', linestyle='--')\n                axes[1].add_patch(rect)\n\n            title_txt = \"Meijering (Center) + Thresh (Edge)\" if self.chk_meij.value else \"Hysteresis (Full)\"\n            axes[1].set_title(f\"{title_txt} (z={z})\", fontsize=13)\n            axes[1].axis('off')\n\n            plt.tight_layout()\n            plt.show()\n\n    def compute_3d_metrics(self, b):\n        with self.out_log:\n            clear_output()\n            print(\"⏳ Running full 3D pipeline...\")\n            \n            # --- 1. Generate Base Mask ---\n            if self.chk_meij.value:\n                print(\"- Strategy: Meijering in Center + Simple Threshold at Edge\")\n                s_min = self.s_msig_min.value\n                s_max = max(s_min + 1, self.s_msig_max.value)\n                f_pad = self.s_mpad.value\n                thresh_tube = self.s_mthresh.value\n                thresh_simple = self.s_thigh.value\n                margin = self.s_margin.value\n                \n                mask3d = np.zeros_like(self.vol, dtype=bool)\n                h, w = self.vol.shape[1], self.vol.shape[2]\n                \n                # Pre-calculate center slice indices\n                if margin == 0:\n                     sl_y, sl_x = slice(None), slice(None)\n                     has_center = True\n                elif margin < h//2 and margin < w//2:\n                     sl_y = slice(margin, h - margin)\n                     sl_x = slice(margin, w - margin)\n                     has_center = True\n                else:\n                     has_center = False\n                \n                for z in range(self.vol.shape[0]):\n                    if z % 10 == 0: print(f\"  Processing slice {z}/{self.vol.shape[0]}...\", end='\\r')\n                    \n                    # A. Default to Simple Threshold\n                    slice_mask = self.vol[z] > thresh_simple\n                    \n                    # B. Apply Meijering to Center\n                    if has_center:\n                        tubeness = apply_meijering_slice_safe(self.vol[z], range(s_min, s_max), filter_padding=f_pad)\n                        if tubeness.max() > 0: tubeness /= tubeness.max()\n                        slice_mask[sl_y, sl_x] = tubeness[sl_y, sl_x] > thresh_tube\n                    \n                    mask3d[z] = slice_mask\n                print(\"\\n  Hybrid Filtering Complete.\")\n            else:\n                print(\"- Applying 3D Hysteresis...\")\n                mask3d = apply_hysteresis_3d(self.vol, self.s_tlow.value, self.s_thigh.value)\n\n            # --- 2. 3D Anisotropic Closing ---\n            z_r, xy_r = self.s_zrad.value, self.s_xyrad.value\n            if z_r > 0 or xy_r > 0:\n                print(f\"- 3D anisotropic closing: z_radius={z_r}, xy_radius={xy_r}\")\n                mask3d = apply_anisotropic_closing_3d(mask3d, z_r, xy_r)\n\n            # --- 3. Dust Removal ---\n            mins = self.s_minsize.value\n            if mins > 0:\n                print(f\"- Removing small components (< {mins} voxels)\")\n                mask3d = remove_small_objects(mask3d, min_size=mins)\n\n            # --- 4. Evaluate ---\n            print(\"- Applying ignore mask (GT == 2)\")\n            valid = (self.gt != 2)\n            pred_eval = mask3d & valid\n            gt_eval   = (self.gt == 1)\n\n            # Save TIFF\n            print(\"- Saving prediction_volume.tif...\")\n            pred_visual = (pred_eval.astype(np.uint8) * 255)\n            imgs = [Image.fromarray(pred_visual[i]) for i in range(pred_visual.shape[0])]\n            imgs[0].save(\"prediction_volume.tif\", compression=\"tiff_deflate\", save_all=True, append_images=imgs[1:])\n\n            # Metrics\n            inter = np.logical_and(pred_eval, gt_eval).sum()\n            s_pred = pred_eval.sum()\n            s_gt   = gt_eval.sum()\n            dice = 2.0 * inter / max(s_pred + s_gt, 1)\n            \n            print(\"- Counting Connected Components (26-connectivity)...\")\n            struct26 = generate_binary_structure(3, 3)\n            _, n_pred = label(pred_eval, structure=struct26)\n            _, n_gt   = label(gt_eval,   structure=struct26)\n            \n            print(\"\\n✅ 3D Metrics (valid region):\")\n            print(f\"Dice           : {dice:.4f}\")\n            print(f\"Components     : Pred={n_pred}, GT={n_gt}\")\n            print(f\"|Δ components| : {abs(n_pred - n_gt)}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:34:06.726120Z","iopub.execute_input":"2026-01-06T15:34:06.726497Z","iopub.status.idle":"2026-01-06T15:34:07.278300Z","shell.execute_reply.started":"2026-01-06T15:34:06.726470Z","shell.execute_reply":"2026-01-06T15:34:07.277126Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"app = KaggleTopologyOptimizer((vol2*(gt2!=2))[:,:,:].astype(np.float32),gt2[:,:,:])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:05:43.563676Z","iopub.execute_input":"2026-01-06T16:05:43.563992Z","iopub.status.idle":"2026-01-06T16:05:44.896157Z","shell.execute_reply.started":"2026-01-06T16:05:43.563970Z","shell.execute_reply":"2026-01-06T16:05:44.894698Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# The competition metric\nwhen you click compute in the above widget with some specifications, the prediction is saved in prediction_volume, then run the below code, with appropriate ground truth path to get you competition metric evaluation","metadata":{}},{"cell_type":"code","source":"import vesuvius_2025_metric_demo as metric\nmetric.score_single_tif(gt2p,\"/kaggle/working/prediction_volume.tif\",2.0)\n# 0.749\n#0.575 (0.35-0.6)\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from topometrics import compute_leaderboard_score\npred = load_tif_volume(\"/kaggle/working/prediction_volume.tif\")\nrep = compute_leaderboard_score(\n    predictions=pred,\n    labels=gt1,\n    dims=(0, 1, 2),\n    spacing=(1.0, 1.0, 1.0),  # (z, y, x)\n    surface_tolerance=2.0,  # in spacing units\n    voi_connectivity=26,\n    voi_transform=\"one_over_one_plus\",\n    voi_alpha=0.3,\n    combine_weights=(0.3, 0.35, 0.35),  # (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\nprint(\"Leaderboard score:\", rep.score)  # scalar in [0,1]\nprint(\"Topo score:\", rep.topo.toposcore)  # [0,1]\nprint(\"Surface Dice:\", rep.surface_dice)  # [0,1]\nprint(\"VOI score:\", rep.voi.voi_score)  # (0,1]\nprint(\"VOI split/merge:\", rep.voi.voi_split, rep.voi.voi_merge)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:07:13.479575Z","iopub.status.idle":"2026-01-06T16:07:13.479904Z","shell.execute_reply.started":"2026-01-06T16:07:13.479752Z","shell.execute_reply":"2026-01-06T16:07:13.479767Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Connected Components","metadata":{}},{"cell_type":"code","source":"P = vol2","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:19:50.181993Z","iopub.execute_input":"2026-01-06T16:19:50.182379Z","iopub.status.idle":"2026-01-06T16:19:50.187419Z","shell.execute_reply.started":"2026-01-06T16:19:50.182343Z","shell.execute_reply":"2026-01-06T16:19:50.186400Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pip install connected-components-3d\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:19:50.719792Z","iopub.execute_input":"2026-01-06T16:19:50.720239Z","iopub.status.idle":"2026-01-06T16:19:55.096528Z","shell.execute_reply.started":"2026-01-06T16:19:50.720212Z","shell.execute_reply":"2026-01-06T16:19:55.095194Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import cc3d\nimport cv2 # You can still import cv2 for other operations\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:19:55.098551Z","iopub.execute_input":"2026-01-06T16:19:55.098858Z","iopub.status.idle":"2026-01-06T16:19:55.104646Z","shell.execute_reply.started":"2026-01-06T16:19:55.098825Z","shell.execute_reply":"2026-01-06T16:19:55.103568Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"P = load_tif_volume(\"/kaggle/working/prediction_volume.tif\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:19:55.105608Z","iopub.execute_input":"2026-01-06T16:19:55.105999Z","iopub.status.idle":"2026-01-06T16:19:55.582660Z","shell.execute_reply.started":"2026-01-06T16:19:55.105974Z","shell.execute_reply":"2026-01-06T16:19:55.581275Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import cc3d\n\nC3D = cc3d.connected_components(\n    P.astype(np.uint8),\n    connectivity=26\n)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:20:02.696088Z","iopub.execute_input":"2026-01-06T16:20:02.696487Z","iopub.status.idle":"2026-01-06T16:20:02.931063Z","shell.execute_reply.started":"2026-01-06T16:20:02.696462Z","shell.execute_reply":"2026-01-06T16:20:02.929993Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"np.unique(C3D)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:20:03.349608Z","iopub.execute_input":"2026-01-06T16:20:03.350399Z","iopub.status.idle":"2026-01-06T16:20:04.180550Z","shell.execute_reply.started":"2026-01-06T16:20:03.350366Z","shell.execute_reply":"2026-01-06T16:20:04.179716Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"labels_out = cc3d.connected_components(P)\n\n# 'labels_out' is a new array where all connected voxels have a unique ID.\n# You can then extract individual components:\nN = np.max(labels_out)\nprint(N)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:20:04.819858Z","iopub.execute_input":"2026-01-06T16:20:04.820188Z","iopub.status.idle":"2026-01-06T16:20:05.070075Z","shell.execute_reply.started":"2026-01-06T16:20:04.820166Z","shell.execute_reply":"2026-01-06T16:20:05.068793Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"zslice = 140","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:31:49.441832Z","iopub.execute_input":"2026-01-06T16:31:49.442214Z","iopub.status.idle":"2026-01-06T16:31:49.446930Z","shell.execute_reply.started":"2026-01-06T16:31:49.442187Z","shell.execute_reply":"2026-01-06T16:31:49.445802Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"img = (labels_out == 7)[zslice]\nplt.imshow(img)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:31:49.873122Z","iopub.execute_input":"2026-01-06T16:31:49.873468Z","iopub.status.idle":"2026-01-06T16:31:50.111697Z","shell.execute_reply.started":"2026-01-06T16:31:49.873446Z","shell.execute_reply":"2026-01-06T16:31:50.110383Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Finding Normal To the flow of sheets","metadata":{}},{"cell_type":"code","source":"import numpy as np\nfrom scipy import ndimage, signal\nfrom skimage.feature import structure_tensor\n\n\ndef dominant_sheet_normal_2d(slice2d, sigma=2.0, thresh=0.5):\n    mask = slice2d > thresh\n    if mask.sum() < 100:\n        raise ValueError(\"Not enough signal\")\n\n    Axx, Axy, Ayy = structure_tensor(slice2d.astype(np.float32), sigma=sigma)\n\n    Jxx = Axx[mask].sum()\n    Jxy = Axy[mask].sum()\n    Jyy = Ayy[mask].sum()\n\n    J = np.array([[Jxx, Jxy],\n                  [Jxy, Jyy]], dtype=np.float64)\n\n    vals, vecs = np.linalg.eigh(J)\n    normal = vecs[:, np.argmax(vals)]\n    normal /= np.linalg.norm(normal)\n    return normal","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nslice2d = vol2[160,:,:]*(gt2[160,:,:]!=2)\nflow = dominant_sheet_normal_2d(slice2d)  \nplt.imshow(slice2d, cmap='gray')\ncx, cy = slice2d.shape[0]//2, slice2d.shape[1]//2\nplt.arrow(cy, cx, flow[1]*50, flow[0]*50, color='red', width=2)\nplt.title(\"Dominant sheet flow direction\")\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:31:53.745706Z","iopub.execute_input":"2026-01-06T16:31:53.746727Z","iopub.status.idle":"2026-01-06T16:31:54.024025Z","shell.execute_reply.started":"2026-01-06T16:31:53.746691Z","shell.execute_reply":"2026-01-06T16:31:54.022666Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(flow)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:31:54.025528Z","iopub.execute_input":"2026-01-06T16:31:54.025897Z","iopub.status.idle":"2026-01-06T16:31:54.031734Z","shell.execute_reply.started":"2026-01-06T16:31:54.025868Z","shell.execute_reply":"2026-01-06T16:31:54.030672Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Finding Seeds","metadata":{}},{"cell_type":"code","source":"import numpy as np\n\ndef cast_ray(img, start, normal, max_len=50, ds=1.0):\n    \"\"\"\n    Cast a ray from start along ±normal.\n\n    Parameters\n    ----------\n    img : 2D ndarray (binary or prob)\n    start : (r0, c0)\n    normal : (ny, nx)  -- SHOULD be unit length\n    max_len : maximum distance (pixels)\n    ds : step size\n\n    Returns\n    -------\n    samples : list of (r, c, value)\n    \"\"\"\n\n    h, w = img.shape\n    r0, c0 = start\n    ny, nx = normal\n\n    samples = []\n\n    for sign in (-1, +1):\n        t = 0.0\n        while t <= max_len:\n            r = int(round(r0 + sign * t * ny))\n            c = int(round(c0 + sign * t * nx))\n\n            if r < 0 or r >= h or c < 0 or c >= w:\n                break\n\n            samples.append((r, c, img[r, c]))\n            t += ds\n\n    return samples\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:31:54.304081Z","iopub.execute_input":"2026-01-06T16:31:54.304917Z","iopub.status.idle":"2026-01-06T16:31:54.313181Z","shell.execute_reply.started":"2026-01-06T16:31:54.304876Z","shell.execute_reply":"2026-01-06T16:31:54.311779Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"cnt = 0\nseeds = []\nstart = (160,160)\nnormal = flow\nray = cast_ray(img , start, normal, max_len=160)\nlast = ()\nfor r, c, v in ray:\n    if v == 1:\n        cnt += 1\n        last = (r,c)\n    else:\n        if cnt > 3:\n            seeds.append((zslice,last[0],last[1]))\n        cnt = 0\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:31:56.652693Z","iopub.execute_input":"2026-01-06T16:31:56.653048Z","iopub.status.idle":"2026-01-06T16:31:56.663095Z","shell.execute_reply.started":"2026-01-06T16:31:56.653023Z","shell.execute_reply":"2026-01-06T16:31:56.662046Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"fig, ax = plt.subplots()\n\nax.imshow(img)\n\n# 5. Overlay the points using scatter plot\n# 'c' sets the color (e.g., red), 's' sets the size of the points\nax.scatter([x[2] for x in seeds], [y[1] for y in seeds], c='red', s=100, marker='o', edgecolors='white', linewidths=1.5, label='Points of Interest')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:31:56.924233Z","iopub.execute_input":"2026-01-06T16:31:56.924627Z","iopub.status.idle":"2026-01-06T16:31:57.159968Z","shell.execute_reply.started":"2026-01-06T16:31:56.924603Z","shell.execute_reply":"2026-01-06T16:31:57.159065Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Geodesic Vornoi\ncredits -> hengchk23","metadata":{}},{"cell_type":"code","source":"#geodesic Voronoi\n\nimport numpy as np\nimport math\n\n\n# import pyvista as pv\n# from matplotlib.colors import ListedColormap\n\n# def show_in_3d(data, is_label=True,  cmap=\"tab20\", title=None):\n#     vol = data.astype(np.float32)  # image  #mask  predict overlay_3d\n#     vol = vol / vol.max()  # scale to [0, 1]\n#     #vol = vol ** 1.15\n\n#     # Wrap as a pyvista UniformGrid\n#     grid = pv.wrap(vol)\n\n#     # ----------------------------------------------------------\n#     # 3. Wrap volume as PyVista grid and attach label scalars\n#     if is_label:\n#         labels = data.astype(np.int32)\n#         grid[\"labels\"] = labels.flatten(order=\"F\")  # PyVista expects Fortran order\n#     else:\n#         pass\n\n#     if cmap=='my_color':\n#         my_color = np.full((256,4),fill_value=64,dtype=np.uint8)\n#         my_color[0]   = [  0,   0,   0, 255]\n#         my_color[1]   = [255, 255,   0, 255]\n#         my_color[2]   = [255,   0, 255, 255]\n#         my_color[3]   = [  0, 255, 255, 255]\n#         my_color[4]   = [  0,   0, 255, 255]\n#         my_color[5]   = [  0, 255,   0, 255]\n#         my_color[6]   = [  0,   0,   0, 255]\n#         my_color[254] = [  0,   0, 255, 255]\n#         my_color[255] = [255,   0,   0, 255]\n#         cmap = ListedColormap(my_color/255)\n\n#     opacity = np.ones(255) * 0.5\n#     opacity[0] = 0.0\n#     opacity[254] = 1.0\n\n    # surface = grid.extract_surface()  # marching cubes surface\n    # surface[\"label\"] = surface.point_data[\"values\"]\n\n    # ----------------------------------------------------------\n    # plotter = pv.Plotter(title=title)\n    # if is_label:\n    #     plotter = pv.Plotter(off_screen=True)\n    #     plotter.add_volume(\n    #         grid,\n    #         shade=True,\n    #         opacity=opacity,\n    #         scalars=\"labels\",\n    #         cmap=cmap,\n    #     )\n        \n    #     plotter.export_html(\"volume.html\")\n\n    # else:\n    #     #not implemented\n    #     pass\n    # plotter.show()\n\n\n#############################################################\ndef chamfer_geodesic_voronoi_3d(mask, seeds):\n    \"\"\"\n    mask: 3D bool array, True = allowed/foreground, False = blocked\n    seeds: list of (z,y,x) integer coords (must lie in mask)\n\n    Returns:\n      dist: float32 array, inf outside mask\n      lab : int32 array, -1 outside mask, otherwise index of nearest seed\n    \"\"\"\n    D, H, W = mask.shape\n    INF = np.float32(1e9)\n\n    dist = np.full((D, H, W), INF, dtype=np.float32)\n    lab  = np.full((D, H, W), -1,  dtype=np.int32)\n\n    # init seeds\n    for i, (z, y, x) in enumerate(seeds):\n        if not mask[z, y, x]:\n            raise ValueError(f\"Seed {i} at {(z,y,x)} is outside mask.\")\n        dist[z, y, x] = 0.0\n        lab[z, y, x] = i\n\n    # 26-neighborhood offsets + chamfer costs\n    # axial: 1, face-diagonal: sqrt(2), body-diagonal: sqrt(3)\n    nbrs = []\n    for dz in (-1, 0, 1):\n        for dy in (-1, 0, 1):\n            for dx in (-1, 0, 1):\n                if dz == 0 and dy == 0 and dx == 0:\n                    continue\n                k = abs(dz) + abs(dy) + abs(dx)\n                if k == 1:\n                    w = 1.0\n                elif k == 2:\n                    w = math.sqrt(2.0)\n                else:\n                    w = math.sqrt(3.0)\n                nbrs.append((dz, dy, dx, np.float32(w)))\n\n    # Split into forward/backward neighbor sets based on scan order.\n    # Forward scan order: z increasing, then y, then x.\n    # So \"already visited\" neighbors are those lexicographically smaller:\n    # (dz < 0) or (dz==0 and dy < 0) or (dz==0 and dy==0 and dx < 0)\n    fwd = [(dz,dy,dx,w) for (dz,dy,dx,w) in nbrs\n           if (dz < 0) or (dz == 0 and dy < 0) or (dz == 0 and dy == 0 and dx < 0)]\n    bwd = [(dz,dy,dx,w) for (dz,dy,dx,w) in nbrs\n           if (dz > 0) or (dz == 0 and dy > 0) or (dz == 0 and dy == 0 and dx > 0)]\n\n    # One relaxation pass\n    def relax(scan_z, scan_y, scan_x, neigh_list):\n        for z in scan_z:\n            for y in scan_y:\n                for x in scan_x:\n                    if not mask[z, y, x]:\n                        continue\n                    best_d = dist[z, y, x]\n                    best_l = lab[z, y, x]\n\n                    for dz, dy, dx, w in neigh_list:\n                        zz, yy, xx = z + dz, y + dy, x + dx\n                        if 0 <= zz < D and 0 <= yy < H and 0 <= xx < W and mask[zz, yy, xx]:\n                            cand = dist[zz, yy, xx] + w\n                            if cand < best_d:\n                                best_d = cand\n                                best_l = lab[zz, yy, xx]\n\n                    dist[z, y, x] = best_d\n                    lab[z, y, x]  = best_l\n\n    # Two-pass chamfer propagation\n    relax(range(D), range(H), range(W), fwd)                 # forward pass\n    relax(range(D-1, -1, -1), range(H-1, -1, -1), range(W-1, -1, -1), bwd)  # backward pass\n\n    return dist, lab\n\n\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:32:10.668424Z","iopub.execute_input":"2026-01-06T16:32:10.669225Z","iopub.status.idle":"2026-01-06T16:32:10.689925Z","shell.execute_reply.started":"2026-01-06T16:32:10.669194Z","shell.execute_reply":"2026-01-06T16:32:10.688463Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\ngg = chamfer_geodesic_voronoi_3d(labels_out == 7,seeds)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:32:11.090381Z","iopub.execute_input":"2026-01-06T16:32:11.091472Z","iopub.status.idle":"2026-01-06T16:32:48.043594Z","shell.execute_reply.started":"2026-01-06T16:32:11.091438Z","shell.execute_reply":"2026-01-06T16:32:48.042376Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.imshow(gg[1][zslice])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:33:27.736968Z","iopub.execute_input":"2026-01-06T16:33:27.737322Z","iopub.status.idle":"2026-01-06T16:33:27.973196Z","shell.execute_reply.started":"2026-01-06T16:33:27.737273Z","shell.execute_reply":"2026-01-06T16:33:27.971391Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.imshow(pred1[zslice])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:32:48.045489Z","iopub.execute_input":"2026-01-06T16:32:48.046197Z","iopub.status.idle":"2026-01-06T16:32:48.283346Z","shell.execute_reply.started":"2026-01-06T16:32:48.046161Z","shell.execute_reply":"2026-01-06T16:32:48.282074Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Euclidean distance based merge removal","metadata":{}},{"cell_type":"code","source":"def euclidean_based(mask, seeds):\n    \"\"\"\n    mask: 3D bool array, True = allowed/foreground, False = blocked\n    seeds: list of (z,y,x) integer coords (must lie in mask)\n\n    Returns:\n      dist: float32 array, inf outside mask\n      lab : int32 array, -1 outside mask, otherwise index of nearest seed\n    \"\"\"\n    D, H, W = mask.shape\n    INF = np.float32(1e9)\n\n    dist = np.full((D, H, W), INF, dtype=np.float32)\n    lab  = np.full((D, H, W), -1,  dtype=np.int32)\n\n    # sanity check seeds\n    for i, (z, y, x) in enumerate(seeds):\n        if not mask[z, y, x]:\n            raise ValueError(f\"Seed {i} at {(z,y,x)} is outside mask.\")\n\n    # pre-pack seeds for speed\n    seeds_arr = np.array(seeds, dtype=np.int32)  # shape (N,3)\n\n    for z in range(D):\n        for y in range(H):\n            for x in range(W):\n                if not mask[z, y, x]:\n                    continue\n\n                best_d = INF\n                best_l = -1\n\n                for i, (sz, sy, sx) in enumerate(seeds_arr):\n                    dz = z - sz\n                    dy = y - sy\n                    dx = x - sx\n                    d = math.sqrt(dz*dz + dy*dy + dx*dx)\n\n                    if d < best_d:\n                        best_d = d\n                        best_l = i\n\n                dist[z, y, x] = best_d\n                lab[z, y, x]  = best_l\n\n    return dist, lab\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:32:48.284666Z","iopub.execute_input":"2026-01-06T16:32:48.284949Z","iopub.status.idle":"2026-01-06T16:32:48.294550Z","shell.execute_reply.started":"2026-01-06T16:32:48.284927Z","shell.execute_reply":"2026-01-06T16:32:48.293260Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\ngg2 = eulcidian_distance(labels_out == 7,seeds)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:46:30.802346Z","iopub.execute_input":"2026-01-06T16:46:30.802714Z","iopub.status.idle":"2026-01-06T16:47:26.588844Z","shell.execute_reply.started":"2026-01-06T16:46:30.802689Z","shell.execute_reply":"2026-01-06T16:47:26.587680Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.imshow(gg2[1][zslice])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T16:47:26.590263Z","iopub.execute_input":"2026-01-06T16:47:26.590743Z","iopub.status.idle":"2026-01-06T16:47:26.826241Z","shell.execute_reply.started":"2026-01-06T16:47:26.590717Z","shell.execute_reply":"2026-01-06T16:47:26.825242Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Find number of layers in the prediction\nsimply load the volume, and pass in the function given below to get the numnber of layers, make sure to not have sparse thin broken predictions, it will work even with merges and thick but sometimes not otherwise.","metadata":{}},{"cell_type":"code","source":"\n\n# ============================================================\n# Ray-box intersection (unchanged)\n# ============================================================\n\ndef get_ray_bounds(y0, x0, dy, dx, Y, X):\n    t_min = -np.inf\n    t_max = np.inf\n\n    if abs(dy) < 1e-9:\n        if not (0 <= y0 < Y):\n            return None\n    else:\n        t1 = (0 - y0) / dy\n        t2 = (Y - 1 - y0) / dy\n        t_min = max(t_min, min(t1, t2))\n        t_max = min(t_max, max(t1, t2))\n\n    if abs(dx) < 1e-9:\n        if not (0 <= x0 < X):\n            return None\n    else:\n        t1 = (0 - x0) / dx\n        t2 = (X - 1 - x0) / dx\n        t_min = max(t_min, min(t1, t2))\n        t_max = min(t_max, max(t1, t2))\n\n    if t_min > t_max:\n        return None\n\n    return t_min, t_max\n\n\n# ============================================================\n# Single-slice ray test (segmentation only)\n# ============================================================\n\ndef count_layers_global_perpendicular_rays(\n    volume,\n    z0,\n    n_rays=300,\n    step=0.5,\n    thresh=0.5,\n    min_prominence=0.01,\n    seed=0,\n):\n    rng = np.random.default_rng(seed)\n\n    Z, Y, X = volume.shape\n    slice2d = volume[z0]\n\n    try:\n        dy, dx = dominant_sheet_normal_2d(slice2d, sigma=1.0, thresh=thresh)\n    except ValueError:\n        return []\n\n    mask = slice2d > thresh\n    coords = np.argwhere(mask)\n    if len(coords) == 0:\n        return []\n\n    counts = []\n    random_indices = rng.integers(len(coords), size=n_rays)\n\n    for idx in random_indices:\n        y0, x0 = coords[idx]\n\n        bounds = get_ray_bounds(y0, x0, dy, dx, Y, X)\n        if bounds is None:\n            continue\n\n        t_start, t_end = bounds\n        t_vals = np.arange(t_start, t_end, step)\n\n        if len(t_vals) < 2:\n            continue\n\n        ys = y0 + t_vals * dy\n        xs = x0 + t_vals * dx\n        zs = np.full_like(t_vals, z0)\n\n        sample_coords = np.stack([zs, ys, xs])\n\n        profile = ndimage.map_coordinates(\n            volume,\n            sample_coords,\n            order=1,\n            mode=\"constant\",\n            cval=0.0,\n        )\n\n        peaks, _ = signal.find_peaks(\n            profile,\n            height=thresh,\n            prominence=min_prominence,\n        )\n\n        if len(peaks) > 0:\n            counts.append(len(peaks))\n\n    return counts\n\n\n# ============================================================\n# Multi-slice aggregation (FINAL ENTRY POINT)\n# ============================================================\n\ndef count_layers_multislice_segmentation(\n    volume,\n    n_slices=20,\n    n_rays_per_slice=300,\n    step=0.5,\n    thresh=0.5,\n    min_prominence=0.01,\n    seed=0,\n):\n    rng = np.random.default_rng(seed)\n    Z, _, _ = volume.shape\n\n    z_choices = rng.choice(\n        np.arange(5, Z - 5),\n        size=min(n_slices, Z - 10),\n        replace=False,\n    )\n\n    all_counts = []\n\n    for z0 in z_choices:\n        counts = count_layers_global_perpendicular_rays(\n            volume=volume,\n            z0=z0,\n            n_rays=n_rays_per_slice,\n            step=step,\n            thresh=thresh,\n            min_prominence=min_prominence,\n            seed=int(rng.integers(1e9)),\n        )\n        all_counts.extend(counts)\n\n    return all_counts\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:46:24.046035Z","iopub.execute_input":"2026-01-06T15:46:24.046387Z","iopub.status.idle":"2026-01-06T15:46:24.064206Z","shell.execute_reply.started":"2026-01-06T15:46:24.046363Z","shell.execute_reply":"2026-01-06T15:46:24.063160Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pred1 = load_tif_volume(\"/kaggle/working/prediction_volume.tif\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:46:24.217490Z","iopub.execute_input":"2026-01-06T15:46:24.217811Z","iopub.status.idle":"2026-01-06T15:46:24.619563Z","shell.execute_reply.started":"2026-01-06T15:46:24.217789Z","shell.execute_reply":"2026-01-06T15:46:24.618554Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\ncounts = count_layers_multislice_segmentation(\n    pred1,\n    n_slices=100,\n    n_rays_per_slice=50,\n    thresh=0.6,\n    \n)\n\nhist = np.bincount(counts)\nK = hist.argmax()\n\nprint(\"Histogram:\", hist)\nprint(\"Predicted number of sheets:\", K)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T15:46:24.696894Z","iopub.execute_input":"2026-01-06T15:46:24.697220Z","iopub.status.idle":"2026-01-06T15:46:26.993512Z","shell.execute_reply.started":"2026-01-06T15:46:24.697196Z","shell.execute_reply":"2026-01-06T15:46:26.992362Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T13:48:53.373180Z","iopub.execute_input":"2026-01-06T13:48:53.373771Z","iopub.status.idle":"2026-01-06T13:48:53.656298Z","shell.execute_reply.started":"2026-01-06T13:48:53.373723Z","shell.execute_reply":"2026-01-06T13:48:53.654865Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T13:53:27.938288Z","iopub.execute_input":"2026-01-06T13:53:27.939439Z","iopub.status.idle":"2026-01-06T13:53:27.947449Z","shell.execute_reply.started":"2026-01-06T13:53:27.939393Z","shell.execute_reply":"2026-01-06T13:53:27.946404Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T13:58:06.323140Z","iopub.execute_input":"2026-01-06T13:58:06.323496Z","iopub.status.idle":"2026-01-06T13:58:06.333212Z","shell.execute_reply.started":"2026-01-06T13:58:06.323474Z","shell.execute_reply":"2026-01-06T13:58:06.331952Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pred1[seeds[0]]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T13:58:07.358465Z","iopub.execute_input":"2026-01-06T13:58:07.359022Z","iopub.status.idle":"2026-01-06T13:58:07.366690Z","shell.execute_reply.started":"2026-01-06T13:58:07.358984Z","shell.execute_reply":"2026-01-06T13:58:07.365422Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"seeds","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T13:58:09.268203Z","iopub.execute_input":"2026-01-06T13:58:09.269358Z","iopub.status.idle":"2026-01-06T13:58:09.275934Z","shell.execute_reply.started":"2026-01-06T13:58:09.269314Z","shell.execute_reply":"2026-01-06T13:58:09.274882Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"flow","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T13:49:23.120075Z","iopub.execute_input":"2026-01-06T13:49:23.120470Z","iopub.status.idle":"2026-01-06T13:49:23.128605Z","shell.execute_reply.started":"2026-01-06T13:49:23.120441Z","shell.execute_reply":"2026-01-06T13:49:23.127447Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T12:55:11.091841Z","iopub.status.idle":"2026-01-06T12:55:11.092108Z","shell.execute_reply.started":"2026-01-06T12:55:11.091965Z","shell.execute_reply":"2026-01-06T12:55:11.091978Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# SDF","metadata":{}},{"cell_type":"code","source":"import numpy as np\nfrom scipy import ndimage\n# shape: (D, H, W)\n# foreground (inside surface) = 1\n# background (outside) = 0\ndef signed_distance_transform(mask):\n    \"\"\"\n    Computes signed distance field for a 3D binary mask.\n\n    Parameters\n    ----------\n    mask : np.ndarray (D, H, W), binary\n        1 = inside object\n        0 = outside object\n\n    Returns\n    -------\n    sdf : np.ndarray (D, H, W), float\n        Negative inside, positive outside, 0 on boundary\n    \"\"\"\n    # Distance to background (outside)\n    dist_out = ndimage.distance_transform_edt(mask == 0)\n\n    # Distance to foreground (inside)\n    dist_in = ndimage.distance_transform_edt(mask == 1)\n\n    # Signed distance\n    sdf = dist_out - dist_in\n\n    return sdf\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T12:55:11.104366Z","iopub.status.idle":"2026-01-06T12:55:11.105082Z","shell.execute_reply.started":"2026-01-06T12:55:11.104595Z","shell.execute_reply":"2026-01-06T12:55:11.104612Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"sdfpred1 = signed_distance_transform(pred1)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-06T12:55:11.107288Z","iopub.status.idle":"2026-01-06T12:55:11.107660Z","shell.execute_reply.started":"2026-01-06T12:55:11.107494Z","shell.execute_reply":"2026-01-06T12:55:11.107509Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.imshow(sdfpred1[zslice])","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}