{"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":[{"sourceType":"competition","sourceId":117682,"databundleVersionId":15062069},{"sourceType":"datasetVersion","sourceId":14407441,"datasetId":9181226,"databundleVersionId":15222839},{"sourceType":"datasetVersion","sourceId":14364531,"datasetId":9172693,"databundleVersionId":15175109},{"sourceType":"datasetVersion","sourceId":14261875,"datasetId":8670484,"databundleVersionId":15061995},{"sourceType":"kernelVersion","sourceId":285233081},{"sourceType":"kernelVersion","sourceId":286093322},{"sourceType":"kernelVersion","sourceId":289436856},{"sourceType":"kernelVersion","sourceId":292281173}],"dockerImageVersionId":31192,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Forward Field Estimation\nI currently can't make this notebook pretty so that it's easy to understand, nor can I try this idea for next 7-10 days, So I'm making this public and will give brief info about the idea, You can ask gpt for explanation, and try similar ideas.\n1. Main aim is to estimate where a particular pixel in slice z will move in slice z+1.\n2. To do so I simply take the 3x3 grid in the next slice, and allocate the current pixel's field as weighted(by inferred probability) direction (2d vector denoting change in x and y), This way we get an approx but noisy field initialization.\n3. Now I iteratively improve this field, by allocating current pixel by weighted average of fields of the 3x3 pixel region in the front slice and previous slice.\n4. After we have an estimated field, We can use it to fill holes, I tried this method for 20-30 mins a few days ago and in that 20-30 mins it looked like it could workout.\nFew directions to try -->\n1. Try updating by taking context of different combinations of nearby sheets.\n2. Maybe predicting a field and then iteratively improving it like this.\n\nI coded the function using gpt and have not checked since so there might be error, Also I'm not aware if this is some well known method.\n                           ","metadata":{}},{"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-14T06:47:07.223864Z","iopub.execute_input":"2026-01-14T06:47:07.224156Z","iopub.status.idle":"2026-01-14T06:47:09.443017Z","shell.execute_reply.started":"2026-01-14T06:47:07.224132Z","shell.execute_reply":"2026-01-14T06:47:09.442008Z"}},"outputs":[],"execution_count":null},{"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-14T06:47:09.445300Z","iopub.execute_input":"2026-01-14T06:47:09.445758Z","iopub.status.idle":"2026-01-14T06:47:09.452359Z","shell.execute_reply.started":"2026-01-14T06:47:09.445732Z","shell.execute_reply":"2026-01-14T06:47:09.451301Z"}},"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-14T06:47:09.453542Z","iopub.execute_input":"2026-01-14T06:47:09.453869Z","iopub.status.idle":"2026-01-14T06:47:09.475458Z","shell.execute_reply.started":"2026-01-14T06:47:09.453847Z","shell.execute_reply":"2026-01-14T06:47:09.474560Z"}},"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-14T06:47:09.476438Z","iopub.execute_input":"2026-01-14T06:47:09.476779Z","iopub.status.idle":"2026-01-14T06:47:09.495147Z","shell.execute_reply.started":"2026-01-14T06:47:09.476752Z","shell.execute_reply":"2026-01-14T06:47:09.494065Z"}},"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-14T06:47:09.496202Z","iopub.execute_input":"2026-01-14T06:47:09.496919Z","iopub.status.idle":"2026-01-14T06:47:09.513090Z","shell.execute_reply.started":"2026-01-14T06:47:09.496883Z","shell.execute_reply":"2026-01-14T06:47:09.512139Z"}},"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-14T06:47:09.513996Z","iopub.execute_input":"2026-01-14T06:47:09.514235Z","iopub.status.idle":"2026-01-14T06:47:09.533669Z","shell.execute_reply.started":"2026-01-14T06:47:09.514214Z","shell.execute_reply":"2026-01-14T06:47:09.532680Z"}},"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-14T06:47:09.536197Z","iopub.execute_input":"2026-01-14T06:47:09.536596Z","iopub.status.idle":"2026-01-14T06:47:16.306282Z","shell.execute_reply.started":"2026-01-14T06:47:09.536572Z","shell.execute_reply":"2026-01-14T06:47:16.305524Z"}},"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-14T06:47:16.307190Z","iopub.execute_input":"2026-01-14T06:47:16.307459Z","iopub.status.idle":"2026-01-14T06:47:19.062616Z","shell.execute_reply.started":"2026-01-14T06:47:16.307438Z","shell.execute_reply":"2026-01-14T06:47:19.061836Z"}},"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-14T06:47:19.063390Z","iopub.execute_input":"2026-01-14T06:47:19.063623Z","iopub.status.idle":"2026-01-14T06:47:19.067583Z","shell.execute_reply.started":"2026-01-14T06:47:19.063604Z","shell.execute_reply":"2026-01-14T06:47:19.066537Z"}},"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-14T06:47:19.068419Z","iopub.execute_input":"2026-01-14T06:47:19.068745Z","iopub.status.idle":"2026-01-14T06:47:19.095546Z","shell.execute_reply.started":"2026-01-14T06:47:19.068716Z","shell.execute_reply":"2026-01-14T06:47:19.094623Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"vol1.shape","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-14T06:47:19.096580Z","iopub.execute_input":"2026-01-14T06:47:19.096861Z","iopub.status.idle":"2026-01-14T06:47:19.113934Z","shell.execute_reply.started":"2026-01-14T06:47:19.096837Z","shell.execute_reply":"2026-01-14T06:47:19.112903Z"}},"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-14T06:47:19.115025Z","iopub.execute_input":"2026-01-14T06:47:19.116006Z","iopub.status.idle":"2026-01-14T06:47:20.597177Z","shell.execute_reply.started":"2026-01-14T06:47:19.115975Z","shell.execute_reply":"2026-01-14T06:47:20.596330Z"}},"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-14T06:47:20.598116Z","iopub.execute_input":"2026-01-14T06:47:20.598592Z","iopub.status.idle":"2026-01-14T06:47:20.607720Z","shell.execute_reply.started":"2026-01-14T06:47:20.598567Z","shell.execute_reply":"2026-01-14T06:47:20.606453Z"}},"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-14T06:47:20.608953Z","iopub.execute_input":"2026-01-14T06:47:20.609376Z","iopub.status.idle":"2026-01-14T06:47:20.998866Z","shell.execute_reply.started":"2026-01-14T06:47:20.609353Z","shell.execute_reply":"2026-01-14T06:47:20.998129Z"}},"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-14T06:47:20.999667Z","iopub.execute_input":"2026-01-14T06:47:20.999938Z","iopub.status.idle":"2026-01-14T06:47:21.315345Z","shell.execute_reply.started":"2026-01-14T06:47:20.999911Z","shell.execute_reply":"2026-01-14T06:47:21.314537Z"}},"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-14T06:47:21.316192Z","iopub.execute_input":"2026-01-14T06:47:21.316446Z","iopub.status.idle":"2026-01-14T06:47:21.611157Z","shell.execute_reply.started":"2026-01-14T06:47:21.316427Z","shell.execute_reply":"2026-01-14T06:47:21.610446Z"}},"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-14T06:47:21.611947Z","iopub.execute_input":"2026-01-14T06:47:21.612216Z","iopub.status.idle":"2026-01-14T06:47:21.616184Z","shell.execute_reply.started":"2026-01-14T06:47:21.612196Z","shell.execute_reply":"2026-01-14T06:47:21.615394Z"}},"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-14T06:47:21.617032Z","iopub.execute_input":"2026-01-14T06:47:21.617689Z","iopub.status.idle":"2026-01-14T06:47:28.195347Z","shell.execute_reply.started":"2026-01-14T06:47:21.617667Z","shell.execute_reply":"2026-01-14T06:47:28.194530Z"}},"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-14T06:47:28.196473Z","iopub.execute_input":"2026-01-14T06:47:28.196919Z","iopub.status.idle":"2026-01-14T06:47:28.656986Z","shell.execute_reply.started":"2026-01-14T06:47:28.196885Z","shell.execute_reply":"2026-01-14T06:47:28.656207Z"}},"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-14T06:47:28.657871Z","iopub.execute_input":"2026-01-14T06:47:28.658347Z","iopub.status.idle":"2026-01-14T06:47:29.938669Z","shell.execute_reply.started":"2026-01-14T06:47:28.658326Z","shell.execute_reply":"2026-01-14T06:47:29.937919Z"}},"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,"execution":{"iopub.status.busy":"2026-01-14T07:28:14.329639Z","iopub.execute_input":"2026-01-14T07:28:14.330356Z","iopub.status.idle":"2026-01-14T07:30:21.777667Z","shell.execute_reply.started":"2026-01-14T07:28:14.330324Z","shell.execute_reply":"2026-01-14T07:30:21.776756Z"}},"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-12T16:21:43.823265Z","iopub.execute_input":"2026-01-12T16:21:43.823597Z","iopub.status.idle":"2026-01-12T16:21:43.828022Z","shell.execute_reply.started":"2026-01-12T16:21:43.823571Z","shell.execute_reply":"2026-01-12T16:21:43.827155Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pip install connected-components-3d\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-12T16:21:44.116467Z","iopub.execute_input":"2026-01-12T16:21:44.117152Z","iopub.status.idle":"2026-01-12T16:21:48.428244Z","shell.execute_reply.started":"2026-01-12T16:21:44.117094Z","shell.execute_reply":"2026-01-12T16:21:48.426766Z"}},"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-12T16:21:48.430191Z","iopub.execute_input":"2026-01-12T16:21:48.430499Z","iopub.status.idle":"2026-01-12T16:21:48.718303Z","shell.execute_reply.started":"2026-01-12T16:21:48.430469Z","shell.execute_reply":"2026-01-12T16:21:48.717325Z"}},"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-14T07:08:03.627902Z","iopub.execute_input":"2026-01-14T07:08:03.628222Z","iopub.status.idle":"2026-01-14T07:08:03.960558Z","shell.execute_reply.started":"2026-01-14T07:08:03.628198Z","shell.execute_reply":"2026-01-14T07:08:03.959754Z"}},"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-14T07:08:03.961983Z","iopub.execute_input":"2026-01-14T07:08:03.962208Z","iopub.status.idle":"2026-01-14T07:08:04.160409Z","shell.execute_reply.started":"2026-01-14T07:08:03.962190Z","shell.execute_reply":"2026-01-14T07:08:04.159669Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"np.unique(C3D)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-14T07:08:04.197957Z","iopub.execute_input":"2026-01-14T07:08:04.198235Z","iopub.status.idle":"2026-01-14T07:08:05.047822Z","shell.execute_reply.started":"2026-01-14T07:08:04.198217Z","shell.execute_reply":"2026-01-14T07:08:05.046996Z"}},"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-14T07:08:05.049260Z","iopub.execute_input":"2026-01-14T07:08:05.050142Z","iopub.status.idle":"2026-01-14T07:08:05.255634Z","shell.execute_reply.started":"2026-01-14T07:08:05.050118Z","shell.execute_reply":"2026-01-14T07:08:05.254572Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"zslice = 140","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-14T07:08:07.769159Z","iopub.execute_input":"2026-01-14T07:08:07.769565Z","iopub.status.idle":"2026-01-14T07:08:07.774281Z","shell.execute_reply.started":"2026-01-14T07:08:07.769480Z","shell.execute_reply":"2026-01-14T07:08:07.773201Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"img = (labels_out)[zslice]\nplt.imshow(img)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-14T07:08:11.920772Z","iopub.execute_input":"2026-01-14T07:08:11.921366Z","iopub.status.idle":"2026-01-14T07:08:12.131036Z","shell.execute_reply.started":"2026-01-14T07:08:11.921337Z","shell.execute_reply":"2026-01-14T07:08:12.130129Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"sheet = labels_out==8","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-12T16:22:36.821232Z","iopub.execute_input":"2026-01-12T16:22:36.821971Z","iopub.status.idle":"2026-01-12T16:22:36.840423Z","shell.execute_reply.started":"2026-01-12T16:22:36.821939Z","shell.execute_reply":"2026-01-12T16:22:36.839528Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.imshow(sheet[40])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-12T16:22:58.695009Z","iopub.execute_input":"2026-01-12T16:22:58.695380Z","iopub.status.idle":"2026-01-12T16:22:58.905089Z","shell.execute_reply.started":"2026-01-12T16:22:58.695355Z","shell.execute_reply":"2026-01-12T16:22:58.904067Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nimport plotly.graph_objects as go\nfrom skimage.measure import marching_cubes\n\ndef binary_volume_to_html(volume, out_file=\"mesh.html\"):\n    verts, faces, _, _ = marching_cubes(volume, level=0.5)\n\n    fig = go.Figure(data=[\n        go.Mesh3d(\n            x=verts[:, 0],\n            y=verts[:, 1],\n            z=verts[:, 2],\n            i=faces[:, 0],\n            j=faces[:, 1],\n            k=faces[:, 2],\n            opacity=1.0\n        )\n    ])\n\n    fig.update_layout(scene=dict(aspectmode=\"data\"))\n    fig.write_html(out_file)\n    print(f\"Saved to {out_file}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-14T06:52:09.119224Z","iopub.execute_input":"2026-01-14T06:52:09.119599Z","iopub.status.idle":"2026-01-14T06:52:09.149704Z","shell.execute_reply.started":"2026-01-14T06:52:09.119564Z","shell.execute_reply":"2026-01-14T06:52:09.148864Z"}},"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-14T07:14:46.159259Z","iopub.execute_input":"2026-01-14T07:14:46.159661Z","iopub.status.idle":"2026-01-14T07:14:46.502828Z","shell.execute_reply.started":"2026-01-14T07:14:46.159632Z","shell.execute_reply":"2026-01-14T07:14:46.501975Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"P[0:2,:,:]=0\nP[:,0:2,:] = 0\nP[:,:,0:2] =0 ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-14T06:56:35.623884Z","iopub.execute_input":"2026-01-14T06:56:35.624186Z","iopub.status.idle":"2026-01-14T06:56:35.631932Z","shell.execute_reply.started":"2026-01-14T06:56:35.624164Z","shell.execute_reply":"2026-01-14T06:56:35.631134Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"binary_volume_to_html(P)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-14T06:56:36.368393Z","iopub.execute_input":"2026-01-14T06:56:36.368814Z","iopub.status.idle":"2026-01-14T06:56:45.891991Z","shell.execute_reply.started":"2026-01-14T06:56:36.368787Z","shell.execute_reply":"2026-01-14T06:56:45.891059Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\n\n# Directions (dy, dx)\nDIRS = np.array([\n    [-1, -1], [0, -1], [1, -1],\n    [-1,  0], [0,  0], [1,  0],\n    [-1,  1], [0,  1], [1,  1],\n], dtype=np.int32)\n\ndef forward_initialize(P):\n    \"\"\"\n    P: (Z, H, W)\n    Returns V: (Z, H, W, 2)\n    \"\"\"\n    Z, H, W = P.shape\n    V = np.zeros((Z, H, W, 2), dtype=np.float32)\n\n    # We compute only for interior pixels\n    ys = slice(1, H-1)\n    xs = slice(1, W-1)\n\n    for z in range(Z - 1):\n        P_next = P[z+1]\n\n        weighted_sum = np.zeros((H, W, 2), dtype=np.float32)\n        weight_sum   = np.zeros((H, W), dtype=np.float32)\n\n        for dy, dx in DIRS:\n            shifted = P_next[1+dy:H-1+dy, 1+dx:W-1+dx]   # (H-2, W-2)\n\n            weight_sum[ys, xs] += shifted\n            weighted_sum[ys, xs, 0] += shifted * dx\n            weighted_sum[ys, xs, 1] += shifted * dy\n\n        mask = weight_sum > 1e-6\n        V[z, mask] = weighted_sum[mask] / weight_sum[mask, None]\n\n    # Copy second-last to last slice\n    V[-1] = V[-2]\n\n    return V\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-12T17:36:27.314843Z","iopub.execute_input":"2026-01-12T17:36:27.315218Z","iopub.status.idle":"2026-01-12T17:36:27.324385Z","shell.execute_reply.started":"2026-01-12T17:36:27.315192Z","shell.execute_reply":"2026-01-12T17:36:27.323341Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\ninit = forward_initialize(vol2)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-12T17:36:27.559593Z","iopub.execute_input":"2026-01-12T17:36:27.560550Z","iopub.status.idle":"2026-01-12T17:36:32.385956Z","shell.execute_reply.started":"2026-01-12T17:36:27.560510Z","shell.execute_reply":"2026-01-12T17:36:32.384922Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import plotly.graph_objects as go\nimport numpy as np\n\ndef show_3d_vector_field(V, P=None, stride=10, max_points=50000):\n    \"\"\"\n    V: (Z, H, W, 2)\n    P: optional (Z, H, W) probability volume\n    \"\"\"\n    Z, H, W, _ = V.shape\n\n    zs, ys, xs = np.mgrid[0:Z, 0:H, 0:W]\n\n    # Subsample\n    zs = zs[::stride, ::stride, ::stride]\n    ys = ys[::stride, ::stride, ::stride]\n    xs = xs[::stride, ::stride, ::stride]\n\n    U = V[::stride, ::stride, ::stride, 0]\n    Vv = V[::stride, ::stride, ::stride, 1]\n    Wv = np.ones_like(U)  # dz = 1 → forward flow\n\n    # Flatten\n    xs = xs.flatten()\n    ys = ys.flatten()\n    zs = zs.flatten()\n    U = U.flatten()\n    Vv = Vv.flatten()\n    Wv = Wv.flatten()\n\n    if P is not None:\n        probs = P[::stride, ::stride, ::stride].flatten()\n        mask = probs > 0.2   # show only meaningful areas\n        xs, ys, zs = xs[mask], ys[mask], zs[mask]\n        U, Vv, Wv = U[mask], Vv[mask], Wv[mask]\n\n    # Limit number of arrows (important for performance)\n    if len(xs) > max_points:\n        idx = np.random.choice(len(xs), max_points, replace=False)\n        xs, ys, zs = xs[idx], ys[idx], zs[idx]\n        U, Vv, Wv = U[idx], Vv[idx], Wv[idx]\n\n    fig = go.Figure(data=go.Cone(\n        x=xs, y=ys, z=zs,\n        u=U, v=Vv, w=Wv,\n        sizemode=\"absolute\",\n        sizeref=2,\n        anchor=\"tail\"\n    ))\n\n    fig.update_layout(\n        scene=dict(\n            xaxis_title=\"X\",\n            yaxis_title=\"Y\",\n            zaxis_title=\"Z (slice)\",\n        ),\n        title=\"3D Vector Field\"\n    )\n\n    fig.show()\n    return fig","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-12T16:55:53.266025Z","iopub.execute_input":"2026-01-12T16:55:53.266415Z","iopub.status.idle":"2026-01-12T16:55:53.277834Z","shell.execute_reply.started":"2026-01-12T16:55:53.266388Z","shell.execute_reply":"2026-01-12T16:55:53.276669Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"fig = show_3d_vector_field(init[:100,:100,:100,:], stride=8)\nfig.write_html(\"vector_field.html\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-12T16:55:54.178041Z","iopub.execute_input":"2026-01-12T16:55:54.179253Z","iopub.status.idle":"2026-01-12T16:55:54.227667Z","shell.execute_reply.started":"2026-01-12T16:55:54.179215Z","shell.execute_reply":"2026-01-12T16:55:54.226804Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"V min:\", refined.min())\nprint(\"V max:\", refined.max())\nprint(\"Mean magnitude:\", np.mean(np.linalg.norm(refined, axis=-1)))\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-12T17:12:19.984807Z","iopub.execute_input":"2026-01-12T17:12:19.985162Z","iopub.status.idle":"2026-01-12T17:12:20.988587Z","shell.execute_reply.started":"2026-01-12T17:12:19.985128Z","shell.execute_reply":"2026-01-12T17:12:20.987502Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nfrom numba import njit, prange\n\nNEIGHBORS = np.array([\n    (0, -1, -1), (0, -1, 0), (0, -1, 1),\n    (0,  0, -1), (0,  0, 0), (0,  0, 1),\n    (0,  1, -1), (0,  1, 0), (0,  1, 1),\n    (-1, 0, 0),\n    (1, 0, 0),\n], dtype=np.int32)\n\n@njit(parallel=True, fastmath=True)\ndef iterative_refine_parallel(V, P, n_iters, alpha=0.2):\n    \"\"\"\n    V: (Z, H, W, 2) initial vector field (float32)\n    P: (Z, H, W) probability volume (float32)\n    alpha: how strongly we anchor to the original initialization\n           0.1–0.3 usually works well\n    \"\"\"\n    Z, H, W, _ = V.shape\n\n    V_init = V.copy()   # anchor (never changes)\n    V = V.copy()\n\n    for _ in range(n_iters):\n        V_new = np.zeros_like(V)\n\n        # Parallelize over voxels\n        for z in prange(1, Z-1):\n            for y in range(1, H-1):\n                for x in range(1, W-1):\n\n                    wx = 0.0\n                    wy = 0.0\n                    wsum = 0.0\n\n                    for k in range(NEIGHBORS.shape[0]):\n                        dz, dy, dx = NEIGHBORS[k]\n\n                        zz = z + dz\n                        yy = y + dy\n                        xx = x + dx\n\n                        w = P[zz, yy, xx]\n\n                        wx += w * V[zz, yy, xx, 0]\n                        wy += w * V[zz, yy, xx, 1]\n                        wsum += w\n\n                    if wsum > 1e-8:\n                        # Diffused update\n                        vx = wx / wsum\n                        vy = wy / wsum\n\n                        # 🔒 Anchor to original initialization\n                        V_new[z, y, x, 0] = alpha * V_init[z, y, x, 0] + (1 - alpha) * vx\n                        V_new[z, y, x, 1] = alpha * V_init[z, y, x, 1] + (1 - alpha) * vy\n                    else:\n                        # fallback: keep previous\n                        V_new[z, y, x, 0] = V[z, y, x, 0]\n                        V_new[z, y, x, 1] = V[z, y, x, 1]\n\n        # Preserve boundary slices\n        V_new[0]  = V[0]\n        V_new[-1] = V[-1]\n\n        V = V_new\n\n    return V\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-12T17:36:36.440784Z","iopub.execute_input":"2026-01-12T17:36:36.441088Z","iopub.status.idle":"2026-01-12T17:36:36.454403Z","shell.execute_reply.started":"2026-01-12T17:36:36.441064Z","shell.execute_reply":"2026-01-12T17:36:36.453391Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\ninit = init.astype(np.float32)\nvol2 = vol2.astype(np.float32)\nrefined = iterative_refine_parallel(init, vol2, n_iters=5)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-12T17:36:36.655024Z","iopub.execute_input":"2026-01-12T17:36:36.655896Z","iopub.status.idle":"2026-01-12T17:36:45.431228Z","shell.execute_reply.started":"2026-01-12T17:36:36.655864Z","shell.execute_reply":"2026-01-12T17:36:45.429945Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"fig = show_3d_vector_field(refined[:100,:100,:100,:], stride=8)\nfig.write_html(\"refined.html\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-12T17:36:45.432552Z","iopub.execute_input":"2026-01-12T17:36:45.432833Z","iopub.status.idle":"2026-01-12T17:36:45.495536Z","shell.execute_reply.started":"2026-01-12T17:36:45.432809Z","shell.execute_reply":"2026-01-12T17:36:45.494390Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"curr = sheet[:100,:100,:100]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-12T17:15:07.408824Z","iopub.execute_input":"2026-01-12T17:15:07.409203Z","iopub.status.idle":"2026-01-12T17:15:07.413826Z","shell.execute_reply.started":"2026-01-12T17:15:07.409171Z","shell.execute_reply":"2026-01-12T17:15:07.412899Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"binary_volume_to_html(curr)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-12T17:15:21.774886Z","iopub.execute_input":"2026-01-12T17:15:21.775255Z","iopub.status.idle":"2026-01-12T17:15:21.844844Z","shell.execute_reply.started":"2026-01-12T17:15:21.775227Z","shell.execute_reply":"2026-01-12T17:15:21.843505Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nfrom numba import njit, prange\n\n@njit(parallel=True)\ndef extract_surface_parallel(P, V, seed_thresh=0.5, n_steps=50):\n    \"\"\"\n    P: (Z, H, W) float32\n    V: (Z, H, W, 2) float32\n    Returns:\n        S: (Z, H, W) uint8 mask (0/1)\n    \"\"\"\n    Z, H, W = P.shape\n\n    # Seeds\n    S = np.zeros(P.shape)\n    S = (P > seed_thresh).astype(np.uint8)\n    seed = S.copy()\n    # Active frontier\n    active = S.copy()\n\n    for _ in range(n_steps):\n        new_active = np.zeros_like(S)\n\n        # Parallel over volume\n        for z in prange(Z - 1):\n            for y in range(H):\n                for x in range(W):\n\n                    if active[z, y, x] == 0:\n                        continue\n\n                    dx = V[z, y, x, 0]\n                    dy = V[z, y, x, 1]\n\n                    ny = int(y + dy + 0.5)\n                    nx = int(x + dx + 0.5)\n                    nz = z + 1\n\n                    if 0 <= ny < H and 0 <= nx < W:\n                        if S[nz, ny, nx] == 0:\n                            S[nz, ny, nx] = 1\n                            new_active[nz, ny, nx] = 1\n\n        # Move frontier\n        active = new_active\n        # Early stop\n        if np.sum(active) == 0:\n            break\n\n    return S, seed","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-12T17:39:55.926524Z","iopub.execute_input":"2026-01-12T17:39:55.926871Z","iopub.status.idle":"2026-01-12T17:39:55.939029Z","shell.execute_reply.started":"2026-01-12T17:39:55.926840Z","shell.execute_reply":"2026-01-12T17:39:55.937927Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"P = P.astype(np.float32)\nV = V.astype(np.float32)\n\nS, seed = extract_surface_parallel(vol2, refined, seed_thresh=0.5, n_steps=10)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-12T17:43:41.621449Z","iopub.execute_input":"2026-01-12T17:43:41.622187Z","iopub.status.idle":"2026-01-12T17:43:42.187869Z","shell.execute_reply.started":"2026-01-12T17:43:41.622150Z","shell.execute_reply":"2026-01-12T17:43:42.186897Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"Seeds:\", (seed).sum())\nprint(\"Extracted:\", S.sum())\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-12T17:43:42.189133Z","iopub.execute_input":"2026-01-12T17:43:42.189548Z","iopub.status.idle":"2026-01-12T17:43:42.255612Z","shell.execute_reply.started":"2026-01-12T17:43:42.189520Z","shell.execute_reply":"2026-01-12T17:43:42.254779Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.imshow(seed[150])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-12T17:43:53.417546Z","iopub.execute_input":"2026-01-12T17:43:53.417955Z","iopub.status.idle":"2026-01-12T17:43:53.639932Z","shell.execute_reply.started":"2026-01-12T17:43:53.417922Z","shell.execute_reply":"2026-01-12T17:43:53.638904Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.imshow(S[150])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-12T17:43:57.054550Z","iopub.execute_input":"2026-01-12T17:43:57.055389Z","iopub.status.idle":"2026-01-12T17:43:57.267829Z","shell.execute_reply.started":"2026-01-12T17:43:57.055353Z","shell.execute_reply":"2026-01-12T17:43:57.267011Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"S.shape","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-12T17:28:29.483850Z","iopub.execute_input":"2026-01-12T17:28:29.484240Z","iopub.status.idle":"2026-01-12T17:28:29.491061Z","shell.execute_reply.started":"2026-01-12T17:28:29.484211Z","shell.execute_reply":"2026-01-12T17:28:29.490203Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}