{"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":14776787,"datasetId":9445668,"databundleVersionId":15629475},{"sourceType":"modelInstanceVersion","sourceId":673516,"databundleVersionId":14789938,"modelInstanceId":499479,"modelId":510647},{"sourceType":"modelInstanceVersion","sourceId":681152,"databundleVersionId":14910215,"modelInstanceId":516822,"modelId":510647},{"sourceType":"modelInstanceVersion","sourceId":665924,"databundleVersionId":14691066,"modelInstanceId":504051,"modelId":510647},{"sourceType":"modelInstanceVersion","sourceId":674747,"databundleVersionId":14811492,"modelInstanceId":503784,"modelId":510647},{"sourceType":"modelInstanceVersion","sourceId":732880,"databundleVersionId":15477237,"modelInstanceId":516822,"modelId":510647},{"sourceType":"modelInstanceVersion","sourceId":660383,"databundleVersionId":14626934,"modelInstanceId":499479,"modelId":510647},{"sourceType":"modelInstanceVersion","sourceId":665589,"databundleVersionId":14687610,"modelInstanceId":503784,"modelId":510647},{"sourceType":"modelInstanceVersion","sourceId":655294,"databundleVersionId":14570138,"modelInstanceId":495238,"modelId":510647},{"sourceType":"modelInstanceVersion","sourceId":672178,"databundleVersionId":14761443,"modelInstanceId":495238,"modelId":510647},{"sourceType":"kernelVersion","sourceId":290917305},{"sourceType":"kernelVersion","sourceId":296543947},{"sourceType":"kernelVersion","sourceId":296804121}],"dockerImageVersionId":31193,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"#### Masking 2.5D dataset before 3d inference - Proof od Concept\n\nThis is a test to check fesability of masking with a 2.5D generated mask dataset before inference of a 3d Model ( simulation of masking of labels) initial result was : 0.545 private score\n\n\nThe ideea is simple: 2.5 ensemble detail volume of sheets with details including exfoliations, 3D model shuld detect inside precalculted mask ( apoximated volume burned sheets) \nwarning the 2.5 models were created early in the competition and the quality may not be optima. \n\n2.5D mask aproximate the volume, and original image is masked with background color for non relevant a areas:\n\nResults with masking : \n* Score with mask : 0.542 **Private score: 0.557**\n* Results without masking: **Private score 0.545 **\n\nVisual result is not much different but private score is improved by 0.012\n\nFull credits for 3D model copied from this notebook: https://www.kaggle.com/code/ibratusmonov/02-02-2026-vesuvius-v3\n\nThe two scripts: masking and 3d inference, were fused using LLM.","metadata":{}},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"code","source":"Copilot=False\nif Copilot:\n    ################################################################################\n    # 0. SETUP & IMPORTS\n    ################################################################################\n    import os\n    import sys\n    import subprocess\n    import time\n    from pathlib import Path\n    from collections import OrderedDict\n    import zipfile\n    import io\n    \n    import numpy as np\n    import torch\n    import torch.nn as nn\n    from PIL import Image, ImageSequence\n    from tqdm import tqdm\n    import tifffile\n    from scipy import ndimage\n    import cv2\n    import matplotlib.pyplot as plt\n    import gc\n    import pandas as pd\n    \n    # Instalare pachete offline pentru 3D (keras, medicai, imagecodecs etc.)\n    WHEEL_DIR = \"/kaggle/input/vsdetection-packages-offline-installer-only/whls\"\n    if os.path.exists(WHEEL_DIR):\n        print(f\"Installing packages from: {WHEEL_DIR}\")\n        subprocess.run(\n            f\"pip install --no-index --find-links {WHEEL_DIR} \"\n            f\"{WHEEL_DIR}/keras_nightly-*.whl \"\n            f\"{WHEEL_DIR}/tifffile-*.whl \"\n            f\"{WHEEL_DIR}/imagecodecs-*.whl \"\n            f\"{WHEEL_DIR}/medicai-*.whl\",\n            shell=True, check=True\n        )\n    else:\n        print(f\"ERROR: Package directory not found at: {WHEEL_DIR}\")\n        sys.exit(1)\n    \n    os.environ[\"KERAS_BACKEND\"] = \"tensorflow\"\n    os.environ[\"TF_CPP_MIN_LOG_LEVEL\"] = \"2\"\n    \n    from keras import ops\n    from medicai.transforms import Compose, NormalizeIntensity\n    from medicai.models import TransUNet\n    from medicai.utils.inference import SlidingWindowInference\n    from skimage.morphology import remove_small_objects\n    import scipy.ndimage as ndi\n    \n    ################################################################################\n    # 1. GLOBAL CONFIG\n    ################################################################################\n    root_dir = \"/kaggle/input/vesuvius-challenge-surface-detection\"\n    TEST_IMG_DIR = Path(root_dir) / \"test_images\"\n    \n    MASK_2D5_DIR = Path(\"/kaggle/working/masks_2_5d\")\n    MASK_2D5_DIR.mkdir(exist_ok=True)\n    \n    MASKED_TEST_DIR = Path(\"/kaggle/working/masked_test_volumes\")\n    MASKED_TEST_DIR.mkdir(exist_ok=True)\n    \n    SUBMISSION_2D5_ZIP = \"/kaggle/working/submission_2_5d.zip\"\n    SUBMISSION_3D_ZIP = \"/kaggle/working/submission_3d.zip\"\n    \n    DEVICE = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n    \n    # ---------------- 2.5D CONFIG ----------------\n    Z_CONTEXT = 3\n    THRESHOLD = 0.055\n    MIN_CC_SIZE = 3\n    TTA_2D5 = True\n    \n    CHECKPOINT_PATHS_2D = [\n        Path(\"/kaggle/input/datasets/crischir/vesuvius-competition-models/results/checkpoints/best_model.pth\"),\n        Path(\"/kaggle/input/datasets/crischir/vesuvius-competition-models/results/checkpoints/model_epoch_20.pth\"),\n        Path(\"/kaggle/input/datasets/crischir/vesuvius-competition-models/results/checkpoints/latest_model.pth\"),\n        Path(\"/kaggle/input/datasets/crischir/vesuvius-competition-models/results/checkpoints/model_epoch_30.pth\"),\n        Path(\"/kaggle/input/datasets/crischir/vesuvius-competition-models/results/checkpoints/model_epoch_40.pth\"),\n        Path(\"/kaggle/input/notebooks/crischir/pix2pix-2-5d-training-discriminator/checkpoints/best_model.pth\"),\n        Path(\"/kaggle/input/notebooks/crischir/pix2pix-2-5d-training-discriminator/checkpoints/model_epoch_40.pth\"),\n        Path(\"/kaggle/input/notebooks/crischir/pix2pix-2-5d-training-discriminator/checkpoints/model_epoch_60.pth\"),\n        Path(\"/kaggle/input/notebooks/crischir/pix2pix-2-5d-training-discriminator/checkpoints/model_epoch_100.pth\"),\n        Path(\"/kaggle/input/notebooks/crischir/pix2pix-2-5d-training-discriminator/checkpoints/model_epoch_100.pth\"),\n    ]\n    \n    # ---------------- 3D CONFIG ----------------\n    NUM_CLASSES = 3\n    PATCH_SIZE = (160, 160, 160)\n    OVERLAP = 0.25\n    BATCH_SIZE = 2\n    \n    T_LOW = 0.45\n    T_HIGH = 0.80\n    Z_RADIUS = 2\n    XY_RADIUS = 0\n    DUST_MIN_SIZE = 300\n    \n    ################################################################################\n    # 2. 2.5D MODEL (Pix2Pix)\n    ################################################################################\n    class Pix2PixGenerator(nn.Module):\n        def __init__(self, in_ch=2 * Z_CONTEXT + 1, out_ch=1):\n            super().__init__()\n            self.e1 = nn.Sequential(\n                nn.Conv2d(in_ch, 64, 4, 2, 1),\n                nn.LeakyReLU(0.2, inplace=True)\n            )\n            self.e2 = nn.Sequential(\n                nn.Conv2d(64, 128, 4, 2, 1),\n                nn.BatchNorm2d(128),\n                nn.LeakyReLU(0.2, inplace=True)\n            )\n            self.e3 = nn.Sequential(\n                nn.Conv2d(128, 256, 4, 2, 1),\n                nn.BatchNorm2d(256),\n                nn.LeakyReLU(0.2, inplace=True)\n            )\n            self.d1 = nn.Sequential(\n                nn.ConvTranspose2d(256, 128, 4, 2, 1),\n                nn.BatchNorm2d(128),\n                nn.ReLU(inplace=True)\n            )\n            self.d2 = nn.Sequential(\n                nn.ConvTranspose2d(256, 64, 4, 2, 1),\n                nn.BatchNorm2d(64),\n                nn.ReLU(inplace=True)\n            )\n            self.final = nn.ConvTranspose2d(128, out_ch, 4, 2, 1)\n    \n        def forward(self, x):\n            en1 = self.e1(x)\n            en2 = self.e2(en1)\n            en3 = self.e3(en2)\n            de1 = self.d1(en3)\n            de2 = self.d2(torch.cat([de1, en2], 1))\n            return self.final(torch.cat([de2, en1], 1))\n    \n    def load_ensemble_2d(paths, device):\n        models = []\n        for path in paths:\n            if not path.exists():\n                print(f\"Warning: Checkpoint {path} not found. Skipping.\")\n                continue\n            model = Pix2PixGenerator().to(device)\n            checkpoint = torch.load(path, map_location=device)\n            state_dict = checkpoint[\"state_dict\"] if \"state_dict\" in checkpoint else checkpoint\n    \n            new_state_dict = OrderedDict()\n            for k, v in state_dict.items():\n                name = k[7:] if k.startswith(\"module.\") else k\n                new_state_dict[name] = v\n    \n            model.load_state_dict(new_state_dict, strict=False)\n            model.eval()\n            models.append(model)\n            print(f\"Loaded 2.5D: {path}\")\n        return models\n    \n    ################################################################################\n    # 3. 2.5D UTILS\n    ################################################################################\n    def read_volume(path: Path):\n        with Image.open(str(path)) as img:\n            frames = [np.array(frame) for frame in ImageSequence.Iterator(img)]\n        vol = np.stack(frames, axis=0)\n        if vol.ndim == 2:\n            vol = vol[None, ...]\n        return vol\n    \n    def apply_cc_filter(mask: np.ndarray, min_size: int = MIN_CC_SIZE):\n        labeled, num_features = ndimage.label(mask)\n        if num_features == 0:\n            return mask\n        sizes = np.bincount(labeled.ravel())\n        mask[sizes[labeled] < min_size] = 0\n        return mask\n    \n    def predict_with_tta_2d5(model, inp, device):\n        logits = []\n    \n        with torch.no_grad():\n            logits.append(torch.sigmoid(model(inp)).cpu().numpy()[0, 0])\n    \n        for dims in [(2,), (3,), (2, 3)]:\n            img_f = torch.flip(inp, dims=dims)\n            with torch.no_grad():\n                p = torch.sigmoid(model(img_f)).cpu().numpy()[0, 0]\n            np_dims = tuple(d - 2 for d in dims)\n            p = np.flip(p, axis=np_dims)\n            logits.append(p)\n    \n        for k in [1, 2, 3]:\n            img_r = torch.rot90(inp, k=k, dims=(2, 3))\n            with torch.no_grad():\n                p = torch.sigmoid(model(img_r)).cpu().numpy()[0, 0]\n            p = np.rot90(p, k=-k, axes=(0, 1))\n            logits.append(p)\n    \n        return np.mean(logits, axis=0)\n    \n    ################################################################################\n    # 4. 2.5D FULL-SLICE INFERENCE\n    ################################################################################\n    def run_full_2d5_inference(models, device, out_dir: Path):\n        test_files = list(TEST_IMG_DIR.glob(\"*.tif\"))\n        out_paths = []\n    \n        for tif_path in test_files:\n            print(f\"--- 2.5D: Processing {tif_path.name} with Ensemble of {len(models)} ---\")\n            vol = read_volume(tif_path)\n            z_max, h, w = vol.shape\n    \n            vol_f = vol.astype(np.float32)\n            v_mean, v_std = vol_f.mean(), vol_f.std() + 1e-6\n            vol_f = (vol_f - v_mean) / v_std\n    \n            vol_padded = np.pad(\n                vol_f,\n                ((Z_CONTEXT, Z_CONTEXT), (0, 0), (0, 0)),\n                mode=\"reflect\"\n            )\n            output_volume = np.zeros((z_max, h, w), dtype=np.uint8)\n    \n            for z in tqdm(range(z_max)):\n                z_start, z_end = z, z + 2 * Z_CONTEXT + 1\n                slice_stack = vol_padded[z_start:z_end]\n    \n                slice_stack = slice_stack.astype(np.float32)\n                s_mean, s_std = slice_stack.mean(), slice_stack.std() + 1e-6\n                slice_stack = (slice_stack - s_mean) / s_std\n    \n                inp = torch.from_numpy(slice_stack).unsqueeze(0).to(device)\n    \n                ensemble_pred = np.zeros((h, w), dtype=np.float32)\n    \n                with torch.no_grad():\n                    for model in models:\n                        if TTA_2D5:\n                            pred = predict_with_tta_2d5(model, inp, device)\n                        else:\n                            pred = torch.sigmoid(model(inp)).cpu().numpy()[0, 0]\n                        ensemble_pred += pred\n    \n                    ensemble_pred /= len(models)\n    \n                mask = (ensemble_pred > THRESHOLD).astype(np.uint8)\n                mask = apply_cc_filter(mask, min_size=MIN_CC_SIZE)\n                output_volume[z] = mask\n    \n            out_name = out_dir / tif_path.name\n            tifffile.imwrite(str(out_name), output_volume, compression=\"deflate\")\n            out_paths.append(out_name)\n    \n        return out_paths\n    \n    ################################################################################\n    # 5. DEBUG: VIZUALIZARE 5 SLICE-URI DIN PRIMUL TIF 2.5D\n    ################################################################################\n    def process_binary_tif_v3_from_tif(tif_path, page_idx=10, min_area=20):\n        vol = tifffile.imread(str(tif_path))\n        if vol.ndim == 2:\n            vol = vol[None, ...]\n        z_max = vol.shape[0]\n        page_idx = max(0, min(page_idx, z_max - 1))\n    \n        raw_frame = vol[page_idx]\n    \n        if raw_frame.max() == 1:\n            frame_8bit = (raw_frame * 255).astype(np.uint8)\n        else:\n            frame_8bit = cv2.normalize(\n                raw_frame, None, 0, 255,\n                cv2.NORM_MINMAX, dtype=cv2.CV_8U\n            )\n    \n        num_labels, labels, stats, _ = cv2.connectedComponentsWithStats(\n            frame_8bit, connectivity=8\n        )\n    \n        keep_labels = np.where(stats[:, cv2.CC_STAT_AREA] >= min_area)[0]\n        keep_labels = keep_labels[keep_labels != 0]\n    \n        mask = np.isin(labels, keep_labels).astype(np.uint8) * 255\n    \n        plt.figure(figsize=(12, 5))\n        plt.subplot(1, 2, 1)\n        plt.imshow(frame_8bit, cmap='gray')\n        plt.title(f\"Original (slice {page_idx})\")\n    \n        plt.subplot(1, 2, 2)\n        plt.imshow(mask, cmap='gray')\n        plt.title(f\"Mask > {min_area}px (slice {page_idx})\")\n        plt.show()\n    \n        return frame_8bit, mask\n    \n    def visualize_first_tif_5_slices(mask_paths):\n        if not mask_paths:\n            print(\"No 2.5D masks found for visualization.\")\n            return\n        first_tif = mask_paths[0]\n        print(f\"Visualizing 5 slices from: {first_tif.name}\")\n    \n        vol = tifffile.imread(str(first_tif))\n        z_max = vol.shape[0]\n        if z_max < 5:\n            indices = list(range(z_max))\n        else:\n            indices = np.linspace(0, z_max - 1, 5, dtype=int)\n    \n        for idx in indices:\n            process_binary_tif_v3_from_tif(first_tif, page_idx=idx, min_area=15)\n    \n    ################################################################################\n    # 6. APLICĂ MĂȘTILE 2.5D PESTE VOLUMELE BRUTE (PATCH LZW CU PIL)\n    ################################################################################\n    def read_tif_any(path):\n        with Image.open(path) as img:\n            frames = [np.array(frame) for frame in ImageSequence.Iterator(img)]\n        return np.stack(frames, axis=0)\n    \n    def apply_2d5_masks_to_raw(mask_dir: Path, raw_dir: Path, out_dir: Path):\n        mask_files = sorted(mask_dir.glob(\"*.tif\"))\n        print(\"Applying 2.5D masks to raw test volumes...\")\n    \n        for mask_path in tqdm(mask_files):\n            raw_path = raw_dir / mask_path.name\n            if not raw_path.exists():\n                print(f\"Raw volume not found for {mask_path.name}, skipping.\")\n                continue\n    \n            raw_vol = read_tif_any(str(raw_path)).astype(np.uint16)\n            mask_vol = tifffile.imread(str(mask_path)).astype(np.uint8)\n    \n            if raw_vol.shape != mask_vol.shape:\n                print(f\"Shape mismatch for {mask_path.name}, skipping.\")\n                continue\n    \n            masked = raw_vol.copy()\n            masked[mask_vol == 0] = 0\n    \n            out_path = out_dir / mask_path.name\n            tifffile.imwrite(str(out_path), masked, compression=\"deflate\")\n    \n        print(\"Masked volumes saved to:\", out_dir)\n    \n    ################################################################################\n    # 7. 3D MODEL (TransUNet)\n    ################################################################################\n    def get_ensemble_models_3d():\n        base_path = \"/kaggle/input/vsd-model/keras/transunet\"\n    \n        path_combo = f\"{base_path}/3/transunet.seresnext50.160px.comboloss.weights.h5\"\n        path_tpu = \"/kaggle/input/train-vesuvius-surface-3d-detection-on-tpu/model.weights.h5\"\n    \n        models = []\n    \n        if os.path.exists(path_combo):\n            print(f\"Loading Model 1 (ComboLoss): {path_combo}\")\n            m = TransUNet(\n                input_shape=(160, 160, 160, 1),\n                encoder_name='seresnext50',\n                classifier_activation=None,\n                num_classes=NUM_CLASSES\n            )\n            m.load_weights(path_combo)\n            models.append(m)\n    \n        if os.path.exists(path_tpu):\n            print(f\"Loading Model 2 (TPU): {path_tpu}\")\n            m = TransUNet(\n                input_shape=(160, 160, 160, 1),\n                encoder_name='seresnext50',\n                classifier_activation=None,\n                num_classes=NUM_CLASSES\n            )\n            m.load_weights(path_tpu)\n            models.append(m)\n    \n        if not models:\n            print(\"WARNING: No weights found. Initializing random model.\")\n            m = TransUNet(\n                input_shape=(160, 160, 160, 1),\n                encoder_name='seresnext50',\n                classifier_activation=None,\n                num_classes=NUM_CLASSES\n            )\n            return [m]\n    \n        return models\n    \n    def val_transformation(image):\n        data = {\"image\": image}\n        pipeline = Compose([\n            NormalizeIntensity(keys=[\"image\"], nonzero=True, channel_wise=False)\n        ])\n        result = pipeline(data)\n        return result[\"image\"]\n    \n    def load_volume_3d(path):\n        vol = tifffile.imread(path).astype(np.float32)\n        vol = vol[None, ..., None]\n        return vol\n    \n    def predict_with_tta_3d(inputs, swi):\n        logits = []\n        logits.append(swi(inputs))\n        for k in [1, 2, 3]:\n            img_r = np.rot90(inputs, k=k, axes=(2, 3))\n            p = swi(img_r)\n            p = np.rot90(p, k=-k, axes=(2, 3))\n            logits.append(p)\n        return np.mean(logits, axis=0)\n    \n    def ensemble_predict_3d(inputs, models):\n        if len(models) == 2:\n            weights = [0.75, 0.25]\n        else:\n            weights = [1.0]\n    \n        ensemble_probs = []\n    \n        for i, model in enumerate(models):\n            swi = SlidingWindowInference(\n                model,\n                num_classes=NUM_CLASSES,\n                roi_size=PATCH_SIZE,\n                sw_batch_size=BATCH_SIZE,\n                mode='gaussian',\n                overlap=OVERLAP\n            )\n    \n            logits = predict_with_tta_3d(inputs, swi)\n            probs = ops.softmax(logits, axis=-1)\n            fg_probs = probs[..., 1]\n            ensemble_probs.append(fg_probs * weights[i])\n    \n        final_probs = np.sum(ensemble_probs, axis=0)\n        return np.squeeze(final_probs)\n    \n    def build_anisotropic_struct(z_radius, xy_radius):\n        z, r = z_radius, xy_radius\n        if z == 0 and r == 0:\n            return None\n        depth = 2 * z + 1\n        size = 2 * r + 1\n        struct = np.zeros((depth, size, size), dtype=bool)\n        cz, cy, cx = z, r, r\n        for dz in range(-z, z + 1):\n            for dy in range(-r, r + 1):\n                for dx in range(-r, r + 1):\n                    if dy * dy + dx * dx <= r * r:\n                        struct[cz + dz, cy + dy, cx + dx] = True\n        return struct\n    \n    def topo_postprocess(probs, T_low, T_high, z_radius, xy_radius, dust_min_size):\n        strong = probs >= T_high\n        weak = probs >= T_low\n    \n        if not strong.any():\n            return np.zeros_like(probs, dtype=np.uint8)\n    \n        struct_hyst = ndi.generate_binary_structure(3, 3)\n        mask = ndi.binary_propagation(strong, mask=weak, structure=struct_hyst)\n    \n        if not mask.any():\n            return np.zeros_like(probs, dtype=np.uint8)\n    \n        struct = build_anisotropic_struct(z_radius, xy_radius)\n        if struct is not None:\n            mask = ndi.binary_closing(mask, structure=struct)\n    \n        if dust_min_size > 0:\n            mask = remove_small_objects(mask.astype(bool), min_size=dust_min_size)\n    \n        return mask.astype(np.uint8)\n    \n    def run_3d_inference_on_masked_volumes(models, masked_dir: Path, zip_path: str):\n        test_df = pd.read_csv(f\"{root_dir}/test.csv\")\n    \n        with zipfile.ZipFile(zip_path, \"w\", compression=zipfile.ZIP_DEFLATED) as z:\n            for idx, row in test_df.iterrows():\n                image_id = row[\"id\"]\n                print(f\"\\n[3D {idx+1}/{len(test_df)}] Processing {image_id}...\")\n    \n                vol_path = masked_dir / f\"{image_id}.tif\"\n                if not vol_path.exists():\n                    print(f\"Masked volume not found for {image_id}, skipping.\")\n                    continue\n    \n                vol = load_volume_3d(str(vol_path))\n                vol = val_transformation(vol)\n    \n                probs = ensemble_predict_3d(vol, models)\n    \n                final_mask = topo_postprocess(\n                    probs,\n                    T_low=T_LOW,\n                    T_high=T_HIGH,\n                    z_radius=Z_RADIUS,\n                    xy_radius=XY_RADIUS,\n                    dust_min_size=DUST_MIN_SIZE\n                )\n    \n                print(f\"    Foreground voxels: {final_mask.sum():,}\")\n    \n                out_path = f\"/kaggle/working/submission_masks/{image_id}.tif\"\n                os.makedirs(\"/kaggle/working/submission_masks\", exist_ok=True)\n                tifffile.imwrite(out_path, final_mask)\n                z.write(out_path, arcname=f\"{image_id}.tif\")\n                os.remove(out_path)\n    \n                del vol, probs, final_mask\n                gc.collect()\n    \n    ################################################################################\n    # 8. MAIN PIPELINE\n    ################################################################################\n    if __name__ == \"__main__\":\n        start_time = time.time()\n    \n        print(\"=\" * 60)\n        print(\"STEP 1: 2.5D ENSEMBLE INFERENCE\")\n        print(\"=\" * 60)\n        ensemble_models_2d5 = load_ensemble_2d(CHECKPOINT_PATHS_2D, DEVICE)\n        if not ensemble_models_2d5:\n            print(\"Error: No 2.5D models loaded. Check your CHECKPOINT_PATHS_2D.\")\n            sys.exit(1)\n    \n        mask_paths = run_full_2d5_inference(ensemble_models_2d5, DEVICE, MASK_2D5_DIR)\n    \n        with zipfile.ZipFile(SUBMISSION_2D5_ZIP, \"w\", compression=zipfile.ZIP_DEFLATED) as z:\n            for p in mask_paths:\n                z.write(str(p), arcname=p.name)\n        print(f\"2.5D masks zipped to: {SUBMISSION_2D5_ZIP}\")\n    \n        print(\"=\" * 60)\n        print(\"STEP 2: DEBUG – 5 SLICES FROM FIRST 2.5D MASK\")\n        print(\"=\" * 60)\n        visualize_first_tif_5_slices(mask_paths)\n    \n        print(\"=\" * 60)\n        print(\"STEP 3: APPLY 2.5D MASKS TO RAW TEST VOLUMES\")\n        print(\"=\" * 60)\n        apply_2d5_masks_to_raw(MASK_2D5_DIR, TEST_IMG_DIR, MASKED_TEST_DIR)\n    \n        print(\"=\" * 60)\n        print(\"STEP 4: 3D ENSEMBLE INFERENCE ON MASKED VOLUMES\")\n        print(\"=\" * 60)\n        models_3d = get_ensemble_models_3d()\n        run_3d_inference_on_masked_volumes(models_3d, MASKED_TEST_DIR, SUBMISSION_3D_ZIP)\n    \n        end_time = time.time()\n        print(\"\\n\" + \"=\" * 60)\n        print(\"PIPELINE COMPLETE\")\n        print(f\"Total execution time: {end_time - start_time:.2f} seconds\")\n        print(f\"2.5D submission: {SUBMISSION_2D5_ZIP}\")\n        print(f\"3D submission:  {SUBMISSION_3D_ZIP}\")\n        print(\"=\" * 60)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-02T20:40:10.631276Z","iopub.execute_input":"2026-03-02T20:40:10.632210Z","iopub.status.idle":"2026-03-02T20:44:59.720346Z","shell.execute_reply.started":"2026-03-02T20:40:10.632156Z","shell.execute_reply":"2026-03-02T20:44:59.719475Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"Copilot=False\nif Copilot:\n    import zipfile\n    import io\n    import numpy as np\n    import tifffile\n    import matplotlib.pyplot as plt\n    import cv2\n    \n    def visualize_tifs_from_zip(zip_path, num_files=1, num_slices=5, min_area=0):\n        \"\"\"\n        Vizualizează TIFF-urile dintr-un ZIP (2.5D sau 3D).\n        \n        Parametri:\n            zip_path   : calea către arhiva ZIP\n            num_files  : câte fișiere .tif să vizualizeze\n            num_slices : câte slice-uri să afișeze din fiecare fișier\n            min_area   : dacă >0, aplică filtrare connected components\n        \"\"\"\n        print(f\"Deschid ZIP: {zip_path}\")\n    \n        with zipfile.ZipFile(zip_path, \"r\") as z:\n            tif_list = [f for f in z.namelist() if f.endswith(\".tif\")]\n            tif_list = sorted(tif_list)\n    \n            if not tif_list:\n                print(\"Nu există fișiere .tif în ZIP.\")\n                return\n    \n            tif_list = tif_list[:num_files]\n    \n            for tif_name in tif_list:\n                print(f\"\\n--- Vizualizare {tif_name} ---\")\n    \n                with z.open(tif_name) as f:\n                    tif_bytes = io.BytesIO(f.read())\n                    vol = tifffile.imread(tif_bytes)\n    \n                if vol.ndim == 2:\n                    vol = vol[None, ...]\n    \n                z_max = vol.shape[0]\n    \n                if z_max <= num_slices:\n                    indices = list(range(z_max))\n                else:\n                    indices = np.linspace(0, z_max - 1, num_slices, dtype=int)\n    \n                for idx in indices:\n                    frame = vol[idx]\n    \n                    if frame.max() <= 1:\n                        frame_8bit = (frame * 255).astype(np.uint8)\n                    else:\n                        frame_8bit = cv2.normalize(frame, None, 0, 255,\n                                                   cv2.NORM_MINMAX, dtype=cv2.CV_8U)\n    \n                    if min_area > 0:\n                        num_labels, labels, stats, _ = cv2.connectedComponentsWithStats(\n                            frame_8bit, connectivity=8\n                        )\n                        keep = np.where(stats[:, cv2.CC_STAT_AREA] >= min_area)[0]\n                        keep = keep[keep != 0]\n                        mask = np.isin(labels, keep).astype(np.uint8) * 255\n                    else:\n                        mask = frame_8bit\n    \n                    plt.figure(figsize=(12, 5))\n                    plt.subplot(1, 2, 1)\n                    plt.imshow(frame_8bit, cmap=\"gray\")\n                    plt.title(f\"{tif_name} – slice {idx} (raw)\")\n    \n                    plt.subplot(1, 2, 2)\n                    plt.imshow(mask, cmap=\"gray\")\n                    plt.title(f\"{tif_name} – slice {idx} (filtered)\")\n                    plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-02T20:59:34.808276Z","iopub.execute_input":"2026-03-02T20:59:34.808602Z","iopub.status.idle":"2026-03-02T20:59:34.818848Z","shell.execute_reply.started":"2026-03-02T20:59:34.808579Z","shell.execute_reply":"2026-03-02T20:59:34.818219Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"Copilot=False\nif Copilot:\n    visualize_tifs_from_zip(\"/kaggle/working/submission_2_5d.zip\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-02T21:00:42.793382Z","iopub.execute_input":"2026-03-02T21:00:42.793693Z","iopub.status.idle":"2026-03-02T21:00:44.544837Z","shell.execute_reply.started":"2026-03-02T21:00:42.793669Z","shell.execute_reply":"2026-03-02T21:00:44.543941Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"Copilot=False\nif Copilot:\n    visualize_tifs_from_zip(\"/kaggle/working/submission_3d.zip\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-02T21:01:24.951898Z","iopub.execute_input":"2026-03-02T21:01:24.952509Z","iopub.status.idle":"2026-03-02T21:01:26.493398Z","shell.execute_reply.started":"2026-03-02T21:01:24.952484Z","shell.execute_reply":"2026-03-02T21:01:26.492616Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"\n================================================================================\n   VESUVIUS COMBINED PIPELINE\n   \n   Flow:\n     1. Run 2.5D Pix2Pix ensemble inference on test TIFs\n     2. Apply process_binary_tif_v3 to clean binary predictions\n        - For the FIRST test file: plot sample slices for visual inspection\n     3. Use cleaned binary masks to zero-out (black) voxels in the original\n        test TIFs before feeding them to the 3D TransUNet model\n     4. Run 3D ensemble inference (V27 Hysteresis) on masked volumes\n================================================================================\n\"\"\"\n\nimport os\nimport subprocess\nimport time\nimport zipfile\nimport io\nimport gc\nfrom pathlib import Path\nfrom collections import OrderedDict\n\nimport numpy as np\nimport torch\nimport torch.nn as nn\nimport tifffile\nimport scipy.ndimage as ndi\nfrom PIL import Image, ImageSequence\nfrom tqdm import tqdm\nimport cv2\nimport matplotlib.pyplot as plt\n\n# ============================================================================\n# INSTALL 3D-MODEL PACKAGES (Kaggle offline)\n# ============================================================================\nvar = \"/kaggle/input/vsdetection-packages-offline-installer-only/whls\"\nif os.path.exists(var):\n    print(f\"Installing packages from: {var}\")\n    subprocess.run(\n        f\"pip install {var}/keras_nightly-*.whl {var}/tifffile-*.whl \"\n        f\"{var}/imagecodecs-*.whl {var}/medicai-*.whl \"\n        \"--no-index --find-links \" + var,\n        shell=True, check=True\n    )\nelse:\n    print(f\"⚠  Package directory not found at {var}. Skipping offline install.\")\n\nos.environ[\"KERAS_BACKEND\"] = \"tensorflow\"\nos.environ[\"TF_CPP_MIN_LOG_LEVEL\"] = \"2\"\n\n# Protobuf patch for older Keras/TF combos\ntry:\n    from google.protobuf import message_factory as _mf\n    if not hasattr(_mf.MessageFactory, \"GetPrototype\"):\n        from google.protobuf.message_factory import GetMessageClass\n        _mf.MessageFactory.GetPrototype = lambda self, d: GetMessageClass(d)\nexcept Exception:\n    pass\n\nimport keras\nfrom keras import ops\nfrom medicai.transforms import Compose, NormalizeIntensity\nfrom medicai.models import TransUNet\nfrom medicai.utils.inference import SlidingWindowInference\nimport pandas as pd\nfrom skimage.morphology import remove_small_objects\n\n# ============================================================================\n# ██████████████████████  CONFIG  ██████████████████████\n# ============================================================================\n\n# ---------- Paths -----------------------------------------------------------\nDATA_PATH        = Path(\"/kaggle/input/vesuvius-challenge-surface-detection\")\nTEST_IMG_DIR     = DATA_PATH / \"test_images\"\nOUTPUT_DIR_2D    = Path(\"/kaggle/working/masks_2d\")        # 2.5D binary predictions\nMASKED_VOL_DIR   = Path(\"/kaggle/working/masked_volumes\")  # original TIFs after masking\nOUTPUT_DIR_3D    = Path(\"/kaggle/working/submission_masks\")\nZIP_PATH_3D      = Path(\"/kaggle/working/submission.zip\")\n\nOUTPUT_DIR_2D.mkdir(parents=True, exist_ok=True)\nMASKED_VOL_DIR.mkdir(parents=True, exist_ok=True)\nOUTPUT_DIR_3D.mkdir(parents=True, exist_ok=True)\n\n# ---------- 2.5D model config -----------------------------------------------\nZ_CONTEXT   = 3\nTHRESHOLD   = 0.055\nMIN_CC_SIZE = 3\nTTA         = True   # Toggle TTA for 2.5D models\n\nCHECKPOINT_PATHS_2D = [\n    Path(\"/kaggle/input/datasets/crischir/vesuvius-competition-models/results/checkpoints/best_model.pth\"),\n    Path(\"/kaggle/input/datasets/crischir/vesuvius-competition-models/results/checkpoints/model_epoch_20.pth\"),\n    Path(\"/kaggle/input/datasets/crischir/vesuvius-competition-models/results/checkpoints/latest_model.pth\"),\n    Path(\"/kaggle/input/datasets/crischir/vesuvius-competition-models/results/checkpoints/model_epoch_30.pth\"),\n    Path(\"/kaggle/input/datasets/crischir/vesuvius-competition-models/results/checkpoints/model_epoch_40.pth\"),\n    Path(\"/kaggle/input/notebooks/crischir/pix2pix-2-5d-training-discriminator/checkpoints/best_model.pth\"),\n    Path(\"/kaggle/input/notebooks/crischir/pix2pix-2-5d-training-discriminator/checkpoints/model_epoch_40.pth\"),\n    Path(\"/kaggle/input/notebooks/crischir/pix2pix-2-5d-training-discriminator/checkpoints/model_epoch_60.pth\"),\n    Path(\"/kaggle/input/notebooks/crischir/pix2pix-2-5d-training-discriminator/checkpoints/model_epoch_100.pth\"),\n    Path(\"/kaggle/input/notebooks/crischir/pix2pix-2-5d-training-discriminator/checkpoints/model_epoch_100.pth\"),\n]\n\nDEVICE = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n\n# ---------- Visualization config (applied to first test file only) -----------\nPLOT_SLICE_INDICES = [10, 25, 50, 75, 100]   # which Z slices to visualize\nPLOT_MIN_AREA      = 15                       # min CC area for visualization\n\n# Dilation applied to 2.5D mask before zeroing original volume (Step 3)\n# Increase these if the 3D model is getting too much of the volume blacked out\nMASK_DILATION_XY   = 8   # pixels — grows each detected blob outward in XY\nMASK_DILATION_Z    = 2   # slices — grows protection up/down through Z\n\n# ---------- 3D model config -------------------------------------------------\nNUM_CLASSES  = 3\nPATCH_SIZE   = (160, 160, 160)\nOVERLAP      = 0.25\nBATCH_SIZE   = 2\n\nT_LOW        = 0.45\nT_HIGH       = 0.80\nZ_RADIUS     = 2\nXY_RADIUS    = 0\nDUST_MIN_SIZE = 300\n\n\n# ============================================================================\n# ██████████████████████  STEP 1 – 2.5D MODELS  ██████████████████████\n# ============================================================================\n\nclass Pix2PixGenerator(nn.Module):\n    def __init__(self, in_ch=7, out_ch=1):\n        super().__init__()\n        self.e1 = nn.Sequential(nn.Conv2d(in_ch, 64, 4, 2, 1), nn.LeakyReLU(0.2, inplace=True))\n        self.e2 = nn.Sequential(nn.Conv2d(64, 128, 4, 2, 1), nn.BatchNorm2d(128), nn.LeakyReLU(0.2, inplace=True))\n        self.e3 = nn.Sequential(nn.Conv2d(128, 256, 4, 2, 1), nn.BatchNorm2d(256), nn.LeakyReLU(0.2, inplace=True))\n        self.d1 = nn.Sequential(nn.ConvTranspose2d(256, 128, 4, 2, 1), nn.BatchNorm2d(128), nn.ReLU(inplace=True))\n        self.d2 = nn.Sequential(nn.ConvTranspose2d(256, 64, 4, 2, 1), nn.BatchNorm2d(64), nn.ReLU(inplace=True))\n        self.final = nn.ConvTranspose2d(128, out_ch, 4, 2, 1)\n\n    def forward(self, x):\n        en1 = self.e1(x)\n        en2 = self.e2(en1)\n        en3 = self.e3(en2)\n        de1 = self.d1(en3)\n        de2 = self.d2(torch.cat([de1, en2], 1))\n        return self.final(torch.cat([de2, en1], 1))\n\n\ndef load_2d_ensemble(paths, device):\n    models = []\n    for path in paths:\n        if not path.exists():\n            print(f\"  ⚠  Checkpoint {path.name} not found. Skipping.\")\n            continue\n        model = Pix2PixGenerator(in_ch=2 * Z_CONTEXT + 1).to(device)\n        ckpt  = torch.load(path, map_location=device)\n        sd    = ckpt.get(\"state_dict\", ckpt)\n        clean = OrderedDict(\n            (k[7:] if k.startswith(\"module.\") else k, v) for k, v in sd.items()\n        )\n        model.load_state_dict(clean)\n        model.eval()\n        models.append(model)\n        print(f\"  ✓ Loaded 2D: {path.name}\")\n    return models\n\n\ndef read_volume_pil(path: Path) -> np.ndarray:\n    \"\"\"Read a multi-frame TIF using PIL (avoids tifffile codec issues).\n    Returns array of shape (Z, H, W).\"\"\"\n    with Image.open(str(path)) as img:\n        frames = []\n        for frame in ImageSequence.Iterator(img):\n            frames.append(np.array(frame.copy()))\n    vol = np.stack(frames, axis=0)\n    return vol[None, ...] if vol.ndim == 2 else vol\n\n\ndef apply_cc_filter(mask: np.ndarray, min_size: int = MIN_CC_SIZE):\n    labeled, n = ndi.label(mask)\n    if n == 0:\n        return mask\n    sizes = np.bincount(labeled.ravel())\n    mask[sizes[labeled] < min_size] = 0\n    return mask\n\n\ndef predict_with_tta_2d(model, inp, device):\n    logits = [torch.sigmoid(model(inp)).cpu().numpy()[0, 0]]\n    for dims in [(2,), (3,), (2, 3)]:\n        img_f = torch.flip(inp, dims=dims)\n        with torch.no_grad():\n            p = torch.sigmoid(model(img_f)).cpu().numpy()[0, 0]\n        p = np.flip(p, axis=tuple(d - 2 for d in dims))\n        logits.append(p)\n    for k in [1, 2, 3]:\n        img_r = torch.rot90(inp, k=k, dims=(2, 3))\n        with torch.no_grad():\n            p = torch.sigmoid(model(img_r)).cpu().numpy()[0, 0]\n        p = np.rot90(p, k=-k, axes=(0, 1))\n        logits.append(p)\n    return np.mean(logits, axis=0)\n\n\ndef run_2d_inference(models_2d, tif_path: Path) -> np.ndarray:\n    \"\"\"Run full 2.5D ensemble on a single TIF. Returns binary uint8 volume.\"\"\"\n    print(f\"  Running 2.5D inference on {tif_path.name} …\")\n    vol        = read_volume_pil(tif_path)\n    z_max, h, w = vol.shape\n\n    vol_f      = vol.astype(np.float32)\n    vol_padded = np.pad(vol_f, ((Z_CONTEXT, Z_CONTEXT), (0, 0), (0, 0)), mode=\"reflect\")\n    output     = np.zeros((z_max, h, w), dtype=np.uint8)\n\n    for z in tqdm(range(z_max), desc=\"    slices\"):\n        slice_stack = vol_padded[z: z + 2 * Z_CONTEXT + 1].astype(np.float32)\n        s_mean, s_std = slice_stack.mean(), slice_stack.std() + 1e-6\n        slice_stack = (slice_stack - s_mean) / s_std\n\n        inp          = torch.from_numpy(slice_stack).unsqueeze(0).to(DEVICE)\n        ens_pred     = np.zeros((h, w), dtype=np.float32)\n\n        with torch.no_grad():\n            for model in models_2d:\n                if TTA:\n                    ens_pred += predict_with_tta_2d(model, inp, DEVICE)\n                else:\n                    ens_pred += torch.sigmoid(model(inp)).cpu().numpy()[0, 0]\n        ens_pred /= len(models_2d)\n\n        mask = (ens_pred > THRESHOLD).astype(np.uint8)\n        mask = apply_cc_filter(mask, min_size=MIN_CC_SIZE)\n        output[z] = mask\n\n    return output\n\n\n# ============================================================================\n# ██████████████████████  STEP 2 – VISUALIZE FIRST FILE  ██████████████████████\n# ============================================================================\n\ndef visualize_binary_predictions(binary_vol: np.ndarray,\n                                  slice_indices: list,\n                                  min_area: int,\n                                  title_prefix: str = \"\"):\n    \"\"\"\n    For each requested slice index, show the raw binary prediction and\n    the cleaned mask after connected-component area filtering.\n    \"\"\"\n    print(f\"\\n  Plotting {len(slice_indices)} sample slices for {title_prefix} …\")\n    n = len(slice_indices)\n    fig, axes = plt.subplots(n, 2, figsize=(10, 4 * n))\n    if n == 1:\n        axes = [axes]\n\n    for row, page_idx in enumerate(slice_indices):\n        page_idx = min(page_idx, binary_vol.shape[0] - 1)\n        raw_frame = binary_vol[page_idx]\n\n        # Scale to 0-255\n        if raw_frame.max() == 1:\n            frame_8bit = (raw_frame * 255).astype(np.uint8)\n        else:\n            frame_8bit = cv2.normalize(raw_frame, None, 0, 255,\n                                       cv2.NORM_MINMAX, dtype=cv2.CV_8U)\n\n        # CC area filter\n        num_labels, labels, stats, _ = cv2.connectedComponentsWithStats(\n            frame_8bit, connectivity=8\n        )\n        keep = np.where(stats[:, cv2.CC_STAT_AREA] >= min_area)[0]\n        keep = keep[keep != 0]\n        mask = np.isin(labels, keep).astype(np.uint8) * 255\n\n        axes[row][0].imshow(frame_8bit, cmap=\"gray\")\n        axes[row][0].set_title(f\"{title_prefix} — Raw slice {page_idx}\")\n        axes[row][0].axis(\"off\")\n\n        axes[row][1].imshow(mask, cmap=\"gray\")\n        axes[row][1].set_title(f\"Cleaned mask (area ≥ {min_area}) — slice {page_idx}\")\n        axes[row][1].axis(\"off\")\n\n    plt.tight_layout()\n    plt.savefig(f\"/kaggle/working/sample_slices_{title_prefix}.png\", dpi=100)\n    plt.show()\n    print(f\"  Figure saved to /kaggle/working/sample_slices_{title_prefix}.png\")\n\n\n# ============================================================================\n# ██████████████████████  STEP 3 – APPLY MASK TO ORIGINAL TIF  ██████████████████████\n# ============================================================================\n\ndef apply_binary_mask_to_volume(original_vol: np.ndarray,\n                                 binary_mask: np.ndarray,\n                                 dilation_xy: int = 8,\n                                 dilation_z: int  = 2) -> np.ndarray:\n    \"\"\"\n    Zero out voxels in the original volume that are CONFIDENTLY background\n    according to the 2.5D binary mask.\n\n    Strategy\n    --------\n    1. Take the raw 2.5D binary mask  (1 = surface detected, 0 = background).\n    2. Dilate it generously in XY and Z so that:\n         - thin / slightly mis-aligned surface detections are not clipped\n         - only regions that are solidly empty in ALL nearby slices get zeroed\n    3. Apply: original voxel -> 0  only where dilated_mask == 0.\n\n    This means the 3D model sees the full original signal everywhere the 2.5D\n    model detected *anything nearby*, and only truly empty space is blacked out.\n\n    Parameters\n    ----------\n    dilation_xy : XY dilation radius in pixels (default 8)\n    dilation_z  : Z dilation radius in slices  (default 2)\n    \"\"\"\n    # Build 3-D dilation structuring element (cylinder: circular XY, full Z)\n    dz, dr = dilation_z, dilation_xy\n    struct = np.zeros((2*dz+1, 2*dr+1, 2*dr+1), dtype=bool)\n    for iz in range(-dz, dz+1):\n        for iy in range(-dr, dr+1):\n            for ix in range(-dr, dr+1):\n                if iy*iy + ix*ix <= dr*dr:\n                    struct[dz+iz, dr+iy, dr+ix] = True\n\n    # Dilate: any foreground voxel \"protects\" its entire neighbourhood\n    fg_bool      = binary_mask.astype(bool)\n    dilated_mask = ndi.binary_dilation(fg_bool, structure=struct)   # (Z, H, W) bool\n\n    # Zero out only voxels that remain background even after dilation\n    masked = original_vol.astype(np.float32).copy()\n    masked[~dilated_mask] = 0.0\n\n    n_zeroed = int((~dilated_mask).sum())\n    n_total  = int(dilated_mask.size)\n    print(f\"  Masking: {n_zeroed:,} / {n_total:,} voxels zeroed \"\n          f\"({100*n_zeroed/n_total:.1f}% blacked out, \"\n          f\"{100*(n_total-n_zeroed)/n_total:.1f}% kept for 3D model)\")\n\n    return masked\n\n\n# ============================================================================\n# ██████████████████████  STEP 4 – 3D TRANSUNET  ██████████████████████\n# ============================================================================\n\ndef get_3d_ensemble_models():\n    base_path   = \"/kaggle/input/vsd-model/keras/transunet\"\n    path_combo  = f\"{base_path}/3/transunet.seresnext50.160px.comboloss.weights.h5\"\n    path_tpu    = \"/kaggle/input/train-vesuvius-surface-3d-detection-on-tpu/model.weights.h5\"\n\n    models = []\n    for label, path in [(\"ComboLoss\", path_combo), (\"TPU\", path_tpu)]:\n        if os.path.exists(path):\n            print(f\"  ✓ Loading 3D Model ({label}): {path}\")\n            m = TransUNet(\n                input_shape=(160, 160, 160, 1),\n                encoder_name=\"seresnext50\",\n                classifier_activation=None,\n                num_classes=NUM_CLASSES,\n            )\n            m.load_weights(path)\n            models.append(m)\n\n    if not models:\n        print(\"  ⚠  No 3D weights found. Using random-init model.\")\n        m = TransUNet(\n            input_shape=(160, 160, 160, 1),\n            encoder_name=\"seresnext50\",\n            classifier_activation=None,\n            num_classes=NUM_CLASSES,\n        )\n        models = [m]\n\n    return models\n\n\ndef val_transformation(image):\n    pipeline = Compose([NormalizeIntensity(keys=[\"image\"], nonzero=True, channel_wise=False)])\n    return pipeline({\"image\": image})[\"image\"]\n\n\ndef predict_with_tta_3d(inputs, swi):\n    logits = [swi(inputs)]\n    for k in [1, 2, 3]:\n        img_r = np.rot90(inputs, k=k, axes=(2, 3))\n        p     = swi(img_r)\n        p     = np.rot90(p, k=-k, axes=(2, 3))\n        logits.append(p)\n    return np.mean(logits, axis=0)\n\n\ndef ensemble_predict_3d(inputs, models_3d):\n    weights = [0.75, 0.25] if len(models_3d) == 2 else [1.0] * len(models_3d)\n    ensemble_probs = []\n\n    for i, model in enumerate(models_3d):\n        swi = SlidingWindowInference(\n            model,\n            num_classes=NUM_CLASSES,\n            roi_size=PATCH_SIZE,\n            sw_batch_size=BATCH_SIZE,\n            mode=\"gaussian\",\n            overlap=OVERLAP,\n        )\n        logits = predict_with_tta_3d(inputs, swi)\n        probs  = ops.softmax(logits, axis=-1)\n        fg     = probs[..., 1]\n        ensemble_probs.append(fg * weights[i])\n\n    return np.squeeze(np.sum(ensemble_probs, axis=0))\n\n\ndef build_anisotropic_struct(z_radius, xy_radius):\n    if z_radius == 0 and xy_radius == 0:\n        return None\n    z, r = z_radius, xy_radius\n    struct = np.zeros((2 * z + 1, 2 * r + 1, 2 * r + 1), dtype=bool)\n    for dz in range(-z, z + 1):\n        for dy in range(-r, r + 1):\n            for dx in range(-r, r + 1):\n                if dy * dy + dx * dx <= r * r:\n                    struct[z + dz, r + dy, r + dx] = True\n    return struct\n\n\ndef topo_postprocess(probs):\n    strong = probs >= T_HIGH\n    weak   = probs >= T_LOW\n\n    if not strong.any():\n        return np.zeros_like(probs, dtype=np.uint8)\n\n    struct_hyst = ndi.generate_binary_structure(3, 3)\n    mask = ndi.binary_propagation(strong, mask=weak, structure=struct_hyst)\n\n    if not mask.any():\n        return np.zeros_like(probs, dtype=np.uint8)\n\n    struct = build_anisotropic_struct(Z_RADIUS, XY_RADIUS)\n    if struct is not None:\n        mask = ndi.binary_closing(mask, structure=struct)\n\n    if DUST_MIN_SIZE > 0:\n        mask = remove_small_objects(mask.astype(bool), min_size=DUST_MIN_SIZE)\n\n    return mask.astype(np.uint8)\n\n\n# ============================================================================\n# ██████████████████████  MAIN PIPELINE  ██████████████████████\n# ============================================================================\n\ndef main():\n    start_time = time.time()\n\n    print(\"=\" * 70)\n    print(\"VESUVIUS COMBINED PIPELINE\")\n    print(\"  1 → 2.5D ensemble inference\")\n    print(\"  2 → Visualise first test file (sample slices)\")\n    print(\"  3 → Apply 2.5D binary mask to original TIFs\")\n    print(\"  4 → 3D TransUNet ensemble inference\")\n    print(\"=\" * 70)\n\n    # ------------------------------------------------------------------ #\n    #  Load 2.5D models                                                    #\n    # ------------------------------------------------------------------ #\n    print(\"\\n[STEP 1] Loading 2.5D Pix2Pix ensemble …\")\n    models_2d = load_2d_ensemble(CHECKPOINT_PATHS_2D, DEVICE)\n    if not models_2d:\n        raise RuntimeError(\n            \"No 2.5D models loaded. Check that CHECKPOINT_PATHS_2D \"\n            \"point to existing files.\"\n        )\n\n    # Gather test files\n    test_files = sorted(TEST_IMG_DIR.glob(\"*.tif\"))\n    print(f\"  Found {len(test_files)} test TIF(s).\")\n\n    # ------------------------------------------------------------------ #\n    #  Per-file: 2.5D inference → optional visualisation → save mask      #\n    # ------------------------------------------------------------------ #\n    binary_masks = {}   # image_id -> np.ndarray (z, h, w) uint8\n\n    for file_idx, tif_path in enumerate(test_files):\n        image_id = tif_path.stem\n        print(f\"\\n── File {file_idx+1}/{len(test_files)}: {tif_path.name} ──\")\n\n        binary_vol = run_2d_inference(models_2d, tif_path)\n        binary_masks[image_id] = binary_vol\n\n        # Save intermediate 2.5D prediction\n        out_mask_path = OUTPUT_DIR_2D / tif_path.name\n        tifffile.imwrite(str(out_mask_path), binary_vol, compression=\"deflate\")\n        print(f\"  2.5D mask saved → {out_mask_path}\")\n\n        # ---- Visualise first file only -------------------------------- #\n        if file_idx == 0:\n            print(f\"\\n[STEP 2] Visualising sample slices for first file: {tif_path.name}\")\n            visualize_binary_predictions(\n                binary_vol,\n                slice_indices=PLOT_SLICE_INDICES,\n                min_area=PLOT_MIN_AREA,\n                title_prefix=image_id,\n            )\n\n    # Free 2.5D models from GPU before loading 3D models\n    del models_2d\n    torch.cuda.empty_cache()\n\n    # ------------------------------------------------------------------ #\n    #  Load 3D models                                                      #\n    # ------------------------------------------------------------------ #\n    print(\"\\n[STEP 3 & 4] Loading 3D TransUNet ensemble …\")\n    models_3d = get_3d_ensemble_models()\n\n    # ------------------------------------------------------------------ #\n    #  Per-file: apply 2.5D mask to original TIF → 3D inference           #\n    # ------------------------------------------------------------------ #\n    test_df = pd.read_csv(f\"{DATA_PATH}/test.csv\")\n\n    with zipfile.ZipFile(str(ZIP_PATH_3D), \"w\", compression=zipfile.ZIP_DEFLATED) as zf:\n        for file_idx, (_, row) in enumerate(test_df.iterrows()):\n            image_id = str(row[\"id\"])\n            tif_path = TEST_IMG_DIR / f\"{image_id}.tif\"\n\n            print(f\"\\n── [3D] File {file_idx+1}/{len(test_df)}: {image_id} ──\")\n\n            # Load original volume via PIL (avoids tifffile codec errors)\n            original_vol = read_volume_pil(tif_path).astype(np.float32)\n\n            # Apply 2.5D binary mask\n            if image_id in binary_masks:\n                print(\"  Applying 2.5D mask to zero-out background voxels …\")\n                masked_vol = apply_binary_mask_to_volume(\n                    original_vol,\n                    binary_masks[image_id],\n                    dilation_xy=MASK_DILATION_XY,\n                    dilation_z=MASK_DILATION_Z,\n                )\n            else:\n                print(f\"  ⚠  No 2.5D mask found for {image_id}. Using raw volume.\")\n                masked_vol = original_vol\n\n            # Save masked volume (optional, for debugging) using PIL\n            masked_path = MASKED_VOL_DIR / f\"{image_id}.tif\"\n            masked_frames = [Image.fromarray(masked_vol[z]) for z in range(masked_vol.shape[0])]\n            masked_frames[0].save(\n                str(masked_path),\n                save_all=True,\n                append_images=masked_frames[1:],\n                compression=\"tiff_deflate\",\n            )\n\n            # Prepare for 3D model: shape must be (1, Z, H, W, 1)\n            # masked_vol is (Z, H, W) float32\n            vol_input = masked_vol[np.newaxis, ..., np.newaxis]   # (1, Z, H, W, 1)\n            vol_input = val_transformation(vol_input)\n\n            # Ensemble 3D inference\n            probs      = ensemble_predict_3d(vol_input, models_3d)\n            final_mask = topo_postprocess(probs)\n\n            print(f\"  Foreground voxels: {final_mask.sum():,}\")\n\n            out_path = OUTPUT_DIR_3D / f\"{image_id}.tif\"\n            tifffile.imwrite(str(out_path), final_mask)\n            zf.write(str(out_path), arcname=f\"{image_id}.tif\")\n            os.remove(str(out_path))\n\n            del original_vol, masked_vol, vol_input, probs, final_mask\n            gc.collect()\n\n    elapsed = time.time() - start_time\n    print(\"\\n\" + \"=\" * 70)\n    print(\"PIPELINE COMPLETE\")\n    print(f\"  Submission zip : {ZIP_PATH_3D}\")\n    print(f\"  2.5D masks     : {OUTPUT_DIR_2D}\")\n    print(f\"  Masked volumes : {MASKED_VOL_DIR}\")\n    print(f\"  Total time     : {elapsed:.1f}s  ({elapsed/60:.1f} min)\")\n    print(\"=\" * 70)\n\n\n# Works both as a script and as a Kaggle notebook cell\nmain()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-02T21:03:21.707388Z","iopub.execute_input":"2026-03-02T21:03:21.708120Z","iopub.status.idle":"2026-03-02T21:07:06.706855Z","shell.execute_reply.started":"2026-03-02T21:03:21.708092Z","shell.execute_reply":"2026-03-02T21:07:06.706116Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import zipfile\nimport io\nimport numpy as np\nimport tifffile\nimport matplotlib.pyplot as plt\nimport cv2\n\ndef visualize_results(zip_path=None,\n                      folder_path=None,\n                      num_files=1,\n                      num_slices=5,\n                      min_area=0,\n                      title_prefix=\"result\"):\n    \"\"\"\n    Vizualizează rezultate 2.5D sau 3D dintr-un ZIP sau folder.\n    \n    Parametri:\n        zip_path    : calea către arhiva ZIP (ex: submission.zip)\n        folder_path : folder cu .tif (ex: masks_2d)\n        num_files   : câte fișiere să vizualizeze\n        num_slices  : câte slice-uri per fișier\n        min_area    : filtrare CC (0 = dezactivat)\n    \"\"\"\n\n    # -------------------------------------------------------------\n    # 1) Colectează fișierele TIFF\n    # -------------------------------------------------------------\n    tif_list = []\n\n    if zip_path:\n        print(f\"Deschid ZIP: {zip_path}\")\n        with zipfile.ZipFile(zip_path, \"r\") as z:\n            tif_list = [(\"zip\", f) for f in z.namelist() if f.endswith(\".tif\")]\n    elif folder_path:\n        print(f\"Deschid folder: {folder_path}\")\n        folder = Path(folder_path)\n        tif_list = [(\"folder\", f) for f in sorted(folder.glob(\"*.tif\"))]\n    else:\n        print(\"⚠ Trebuie să specifici zip_path sau folder_path.\")\n        return\n\n    if not tif_list:\n        print(\"⚠ Nu există fișiere .tif de vizualizat.\")\n        return\n\n    tif_list = tif_list[:num_files]\n\n    # -------------------------------------------------------------\n    # 2) Vizualizare\n    # -------------------------------------------------------------\n    for source_type, tif_name in tif_list:\n        print(f\"\\n--- Vizualizare {tif_name} ---\")\n\n        # Citire TIFF\n        if source_type == \"zip\":\n            with zipfile.ZipFile(zip_path, \"r\") as z:\n                with z.open(tif_name) as f:\n                    vol = tifffile.imread(io.BytesIO(f.read()))\n        else:\n            vol = tifffile.imread(str(tif_name))\n\n        if vol.ndim == 2:\n            vol = vol[None, ...]\n\n        z_max = vol.shape[0]\n\n        # Selectăm slice-uri\n        if z_max <= num_slices:\n            indices = list(range(z_max))\n        else:\n            indices = np.linspace(0, z_max - 1, num_slices, dtype=int)\n\n        # ---------------------------------------------------------\n        # Plot\n        # ---------------------------------------------------------\n        for idx in indices:\n            frame = vol[idx]\n\n            # Normalizare\n            if frame.max() <= 1:\n                frame_8bit = (frame * 255).astype(np.uint8)\n            else:\n                frame_8bit = cv2.normalize(frame, None, 0, 255,\n                                           cv2.NORM_MINMAX, dtype=cv2.CV_8U)\n\n            # Filtrare CC\n            if min_area > 0:\n                num_labels, labels, stats, _ = cv2.connectedComponentsWithStats(\n                    frame_8bit, connectivity=8\n                )\n                keep = np.where(stats[:, cv2.CC_STAT_AREA] >= min_area)[0]\n                keep = keep[keep != 0]\n                mask = np.isin(labels, keep).astype(np.uint8) * 255\n            else:\n                mask = frame_8bit\n\n            # Plot\n            plt.figure(figsize=(12, 5))\n            plt.subplot(1, 2, 1)\n            plt.imshow(frame_8bit, cmap=\"gray\")\n            plt.title(f\"{title_prefix} – slice {idx} (raw)\")\n            plt.axis(\"off\")\n\n            plt.subplot(1, 2, 2)\n            plt.imshow(mask, cmap=\"gray\")\n            plt.title(f\"{title_prefix} – slice {idx} (filtered)\")\n            plt.axis(\"off\")\n\n            plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-02T21:11:03.683885Z","iopub.execute_input":"2026-03-02T21:11:03.684608Z","iopub.status.idle":"2026-03-02T21:11:03.696028Z","shell.execute_reply.started":"2026-03-02T21:11:03.684578Z","shell.execute_reply":"2026-03-02T21:11:03.695319Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# visualize_results(\n#     folder_path=\"/kaggle/working/masks_2d\",\n#     num_files=1,\n#     num_slices=5,\n#     min_area=15,\n#     title_prefix=\"2.5D\"\n# )\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-02T21:11:16.471698Z","iopub.execute_input":"2026-03-02T21:11:16.472432Z","iopub.status.idle":"2026-03-02T21:11:17.560025Z","shell.execute_reply.started":"2026-03-02T21:11:16.472405Z","shell.execute_reply":"2026-03-02T21:11:17.559489Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# visualize_results(\n#     zip_path=\"/kaggle/working/submission.zip\",\n#     num_files=1,\n#     num_slices=5,\n#     min_area=0,\n#     title_prefix=\"3D\"\n# )","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-02T21:11:26.871738Z","iopub.execute_input":"2026-03-02T21:11:26.872258Z","iopub.status.idle":"2026-03-02T21:11:27.761691Z","shell.execute_reply.started":"2026-03-02T21:11:26.872235Z","shell.execute_reply":"2026-03-02T21:11:27.761133Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}