{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":117682,"databundleVersionId":14443416,"sourceType":"competition"}],"dockerImageVersionId":31192,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# ============================================================\n#  Vesuvius Challenge – Enhanced Surface Detection\n#  Competition Submission: Fixed & Optimized\n# ============================================================\n\nimport numpy as np\nimport pandas as pd\nfrom pathlib import Path\nfrom PIL import Image, ImageSequence\nimport zipfile\nfrom io import BytesIO\nfrom scipy import ndimage\nfrom tqdm import tqdm\nfrom skimage import measure, filters\n\n# ------------------------------------------------------------\n#  OPTIMIZED PARAMETERS\n# ------------------------------------------------------------\nBASE_PERCENTILE = 25.0           # Base brightness threshold\nMIN_COMPONENT_VOXELS = 200       # Minimum component size\nCLOSING_ITERS = 1                # Gap filling iterations\nADAPTIVE_BLOCK_SIZE = 51         # Local thresholding area\nMIN_SOLIDITY = 0.5               # Shape compactness requirement\nMIN_EXTENT = 0.3                 # Bounding box fill requirement\n\n# ------------------------------------------------------------\n#  COMPETITION PATHS\n# ------------------------------------------------------------\nTEST_DIR = Path(\"/kaggle/input/vesuvius-challenge-surface-detection/test_images/\")\nTEST_CSV = Path(\"/kaggle/input/vesuvius-challenge-surface-detection/test.csv\")\nOUT_ZIP = Path(\"/kaggle/working/submission.zip\")  # Competition requirement\n\n# ------------------------------------------------------------\n#  Load test metadata\n# ------------------------------------------------------------\ntest_df = pd.read_csv(TEST_CSV)\nprint(f\"Found {len(test_df)} test volumes to process.\")\n\n# ------------------------------------------------------------\n#  Load 3D TIFF stack\n# ------------------------------------------------------------\ndef load_stack(path: Path) -> np.ndarray:\n    \"\"\"\n    Read multi-page TIFF into 3D numpy array (Slices, Height, Width)\n    \"\"\"\n    try:\n        with Image.open(path) as tif:\n            arr = [np.array(frame) for frame in ImageSequence.Iterator(tif)]\n        return np.stack(arr)\n    except Exception as e:\n        raise RuntimeError(f\"Failed to load {path}: {e}\")\n\n# ------------------------------------------------------------\n#  Enhanced preprocessing with background normalization\n# ------------------------------------------------------------\ndef enhanced_preprocessing(volume: np.ndarray) -> np.ndarray:\n    \"\"\"\n    Multi-stage preprocessing:\n    1. 3D Gaussian smoothing for noise reduction\n    2. Background subtraction to normalize illumination\n    3. Contrast enhancement\n    \"\"\"\n    # Step 1: Mild 3D Gaussian smoothing\n    smoothed = ndimage.gaussian_filter(volume.astype(np.float32), sigma=(0.3, 1.2, 1.2))\n    \n    # Step 2: Background subtraction (rolling ball approximation)\n    background = ndimage.uniform_filter(smoothed, size=(3, 35, 35))\n    normalized = smoothed - (background * 0.85)\n    \n    # Step 3: Clip extreme values and enhance contrast\n    normalized = np.clip(normalized, np.percentile(normalized, 1), np.percentile(normalized, 99))\n    \n    return normalized.astype(np.uint16)\n\n# ------------------------------------------------------------\n#  Adaptive multi-threshold segmentation\n# ------------------------------------------------------------\ndef adaptive_segmentation(volume: np.ndarray) -> np.ndarray:\n    \"\"\"\n    Hybrid thresholding approach:\n    - Global percentile threshold (base)\n    - Local adaptive threshold (detail preservation)\n    - Combined via logical OR\n    \"\"\"\n    masks = []\n    \n    # Calculate volume-wide statistics for adaptive parameters\n    volume_std = np.std(volume)\n    slice_stds = [np.std(sl) for sl in volume]\n    mean_slice_std = np.mean(slice_stds)\n    \n    for i, sl in enumerate(volume):\n        # Dynamic percentile based on slice characteristics\n        slice_complexity = (slice_stds[i] - mean_slice_std) / (volume_std + 1e-6)\n        dynamic_percentile = BASE_PERCENTILE + slice_complexity * 8\n        dynamic_percentile = max(18.0, min(35.0, dynamic_percentile))\n        \n        # Global threshold\n        global_thresh = np.percentile(sl, 100 - dynamic_percentile)\n        global_mask = sl > global_thresh\n        \n        # Local adaptive threshold (for textured regions)\n        try:\n            local_thresh = filters.threshold_local(sl, block_size=ADAPTIVE_BLOCK_SIZE, method='gaussian')\n            local_mask = sl > (local_thresh * 0.95)\n        except:\n            local_mask = np.zeros_like(global_mask)\n        \n        # Combine masks: global OR local\n        combined_mask = global_mask | local_mask\n        masks.append(combined_mask.astype(np.uint8))\n    \n    return np.stack(masks)\n\n# ------------------------------------------------------------\n#  Surface-aware component cleaning - FIXED VERSION\n# ------------------------------------------------------------\ndef surface_aware_cleaning(mask: np.ndarray, original_volume: np.ndarray) -> np.ndarray:\n    \"\"\"\n    Intelligent component filtering:\n    - Size-based filtering\n    - Intensity-weighted importance  \n    - Shape characteristics (solidity, extent)\n    \"\"\"\n    # Use 6-connectivity for 3D labeling\n    structure = ndimage.generate_binary_structure(rank=3, connectivity=1)\n    labeled, num_components = ndimage.label(mask, structure=structure)\n    \n    if num_components == 0:\n        return mask.astype(np.uint8)\n\n    # CRITICAL FIX: Define intensity thresholds before using them\n    intensity_70 = np.percentile(original_volume, 70)\n    intensity_85 = np.percentile(original_volume, 85)\n    \n    final_mask = np.zeros_like(mask, dtype=bool)\n    \n    for label in range(1, num_components + 1):\n        component_mask = (labeled == label)\n        component_voxels = np.sum(component_mask)\n        \n        # Basic size filter - remove very small components\n        if component_voxels < 50:\n            continue\n            \n        # Calculate component intensity statistics\n        component_intensities = original_volume[component_mask]\n        mean_intensity = np.mean(component_intensities)\n        max_intensity = np.max(component_intensities)\n        \n        # Size-based rules\n        if component_voxels >= 1000:  # Large components - keep most\n            if mean_intensity > intensity_70 * 0.7:\n                final_mask |= component_mask\n            continue\n        \n        # For medium components, check shape and intensity\n        if component_voxels >= 200:\n            # Calculate 2D shape properties from the most representative slice\n            slice_sums = np.sum(component_mask, axis=(1, 2))\n            main_slice = np.argmax(slice_sums)\n            \n            slice_mask = component_mask[main_slice]\n            if np.any(slice_mask):\n                try:\n                    labeled_slice = measure.label(slice_mask)\n                    regions = measure.regionprops(labeled_slice)\n                    if regions:\n                        region = regions[0]\n                        solidity = region.solidity\n                        extent = region.extent\n                        \n                        # Keep components with good shape OR high intensity\n                        if (solidity > MIN_SOLIDITY and extent > MIN_EXTENT) or mean_intensity > intensity_85:\n                            final_mask |= component_mask\n                except:\n                    # Fallback: keep based on intensity alone\n                    if mean_intensity > intensity_85:\n                        final_mask |= component_mask\n        else:\n            # Small components: only keep if very bright\n            if mean_intensity > intensity_85 or max_intensity > np.percentile(original_volume, 90):\n                final_mask |= component_mask\n    \n    return final_mask.astype(np.uint8)\n\n# ------------------------------------------------------------\n#  Direction-aware morphological processing\n# ------------------------------------------------------------\ndef advanced_morphology(mask: np.ndarray) -> np.ndarray:\n    \"\"\"\n    Enhanced morphological operations:\n    - Directional opening to preserve stroke-like structures\n    - Multi-scale closing for gap filling\n    \"\"\"\n    # Create directional structuring elements\n    horizontal_se = np.zeros((1, 1, 5), dtype=bool)  # Emphasize horizontal strokes\n    horizontal_se[0, 0, :] = True\n    \n    vertical_se = np.zeros((1, 5, 1), dtype=bool)    # Emphasize vertical strokes  \n    vertical_se[0, :, 0] = True\n    \n    # Apply directional opening to clean noise while preserving strokes\n    horizontal_cleaned = ndimage.binary_opening(mask, structure=horizontal_se)\n    vertical_cleaned = ndimage.binary_opening(mask, structure=vertical_se)\n    \n    # Combine directional results\n    directional_combined = horizontal_cleaned | vertical_cleaned\n    \n    # Preserve original mask in areas where directional cleaning removed too much\n    preserved_original = mask & (ndimage.binary_dilation(directional_combined, structure=horizontal_se))\n    \n    # Final gentle closing to connect nearby components\n    small_structure = ndimage.generate_binary_structure(rank=3, connectivity=1)\n    if CLOSING_ITERS > 0:\n        final_mask = ndimage.binary_closing(\n            preserved_original, structure=small_structure, iterations=CLOSING_ITERS\n        )\n    else:\n        final_mask = preserved_original\n    \n    # Final size-based cleanup\n    labeled, num = ndimage.label(final_mask, structure=small_structure)\n    if num > 0:\n        component_sizes = ndimage.sum(final_mask, labeled, index=np.arange(1, num + 1))\n        keep_labels = np.where(component_sizes >= MIN_COMPONENT_VOXELS)[0] + 1\n        final_mask = np.isin(labeled, keep_labels)\n    \n    return final_mask.astype(np.uint8)\n\n# ------------------------------------------------------------\n#  Main processing pipeline\n# ------------------------------------------------------------\ndef process_volume(volume: np.ndarray) -> np.ndarray:\n    \"\"\"\n    Enhanced processing pipeline:\n    1. Advanced preprocessing\n    2. Adaptive segmentation  \n    3. Surface-aware cleaning\n    4. Directional morphology\n    \"\"\"\n    print(\"  → Preprocessing...\")\n    processed_vol = enhanced_preprocessing(volume)\n    \n    print(\"  → Segmenting...\")\n    raw_mask = adaptive_segmentation(processed_vol)\n    \n    print(\"  → Cleaning components...\")\n    cleaned_mask = surface_aware_cleaning(raw_mask, processed_vol)\n    \n    print(\"  → Morphological refinement...\")\n    final_mask = advanced_morphology(cleaned_mask)\n    \n    return final_mask.astype(np.uint8)\n\n# ------------------------------------------------------------\n#  Write 3D mask to multi-page TIFF inside ZIP\n# ------------------------------------------------------------\ndef write_stack_to_zip(array3d: np.ndarray, zip_handle, name: str):\n    \"\"\"\n    Save 3D uint8 array as multi-page TIFF directly into ZIP archive\n    \"\"\"\n    try:\n        pages = [Image.fromarray(s.astype(np.uint8)) for s in array3d]\n\n        buffer = BytesIO()\n        pages[0].save(\n            buffer,\n            format=\"TIFF\",\n            save_all=True,\n            append_images=pages[1:],\n            compression=\"tiff_lzw\"  # Reduce file size\n        )\n        zip_handle.writestr(name, buffer.getvalue())\n\n    except Exception as e:\n        raise RuntimeError(f\"Could not write {name} to ZIP: {e}\")\n\n# ------------------------------------------------------------\n#  MAIN PROCESSING - Create submission.zip\n# ------------------------------------------------------------\nprint(\"🚀 Starting Vesuvius Challenge Surface Detection Pipeline...\")\nprint(\"=\" * 60)\n\nwith zipfile.ZipFile(OUT_ZIP, \"w\", compression=zipfile.ZIP_DEFLATED) as z:\n    for _, row in tqdm(test_df.iterrows(), total=len(test_df), desc=\"Processing Volumes\"):\n        \n        vol_id = row[\"id\"]\n        fname = f\"{vol_id}.tif\"\n        fpath = TEST_DIR / fname\n\n        # Check if file exists\n        if not fpath.exists():\n            print(f\"⚠ Missing file: {fname}\")\n            continue\n\n        try:\n            print(f\"\\n📖 Processing {fname}...\")\n            \n            # Load volume\n            vol = load_stack(fpath)\n\n            # Validate dimensions\n            if vol.ndim != 3:\n                print(f\"⚠ Bad shape {vol.shape}, skipping {fname}\")\n                continue\n\n            # Process through pipeline\n            mask = process_volume(vol)\n            \n            # Write to submission zip\n            write_stack_to_zip(mask, z, fname)\n\n            print(f\"✅ Completed {fname} | Input: {vol.shape} | Output: {mask.shape}\")\n\n        except Exception as e:\n            print(f\"❌ Error on {fname}: {e}\")\n            continue\n\nprint(\"=\" * 60)\nprint(f\"🎉 PROCESSING COMPLETE!\")\nprint(f\"📦 Submission file created: {OUT_ZIP}\")\nprint(f\"📊 Processed {len(test_df)} test volumes\")\nprint(\"Ready for competition submission! 🚀\")","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-11-20T09:59:04.688519Z","iopub.execute_input":"2025-11-20T09:59:04.688856Z"}},"outputs":[],"execution_count":null}]}