{"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":"\"\"\"\nVESUVIUS CHALLENGE - TOPOLOGY-AWARE SOLUTION\nOptimized for: SurfaceDice@τ + VOI + TopoScore\nTarget: Score > 0.4 (currently 0.32)\n\nScoring breakdown:\n- 30% TopoScore (preserve topology: no spurious holes/handles)\n- 35% SurfaceDice@τ (surface proximity within τ=2.0)\n- 35% VOI_score (avoid splits/merges of components)\n\"\"\"\n\nimport os\nimport warnings\nwarnings.filterwarnings('ignore')\n\nos.environ[\"KERAS_BACKEND\"] = \"jax\"\nos.environ['XLA_PYTHON_CLIENT_PREALLOCATE'] = 'false'\n\nimport keras\nimport numpy as np\nimport pandas as pd\nimport zipfile\nimport tifffile\nfrom pathlib import Path\nfrom tqdm import tqdm\n\nfrom skimage.morphology import (\n    remove_small_objects, binary_closing, binary_opening,\n    binary_dilation, binary_erosion, ball\n)\nfrom skimage.filters import gaussian, sobel, threshold_otsu\nfrom scipy import ndimage\nfrom scipy.ndimage import uniform_filter, distance_transform_edt\n\nprint(\"=\"*60)\nprint(\"VESUVIUS - TOPOLOGY-AWARE SOLUTION\")\nprint(\"=\"*60)\nprint(f\"Keras: {keras.version()}\")\nprint(\"=\"*60)\n\n# ============================================================================\n# TOPOLOGY-AWARE CONFIGURATION\n# ============================================================================\n\nclass TopoConfig:\n    # Paths\n    ROOT_DIR = Path(\"/kaggle/input/vesuvius-challenge-surface-detection\")\n    TEST_DIR = ROOT_DIR / \"test_images\"\n    OUTPUT_DIR = Path(\"/kaggle/working/submission_masks\")\n    ZIP_PATH = Path(\"/kaggle/working/submission.zip\")\n    \n    # Detection parameters - OPTIMIZED for SurfaceDice@τ=2.0\n    GRADIENT_PERCENTILE = 50        # Balanced\n    INTENSITY_PERCENTILE = 53       # Moderate\n    VARIANCE_PERCENTILE = 48        # Sensitive to texture\n    MORPHOLOGY_PERCENTILE = 60      # Conservative\n    \n    # Voting - Focus on CONSENSUS\n    MIN_METHODS_AGREE = 2           # 2 out of 4\n    GRADIENT_WEIGHT = 0.40          # Edges most important for surface\n    INTENSITY_WEIGHT = 0.30         # Important\n    VARIANCE_WEIGHT = 0.20          # Texture\n    MORPHOLOGY_WEIGHT = 0.10        # Supplementary\n    SOFT_VOTING_THRESHOLD = 0.40    # 40% weighted confidence\n    \n    # TOPOLOGY PRESERVATION (Critical for TopoScore)\n    PRESERVE_TOPOLOGY = True\n    MIN_COMPONENT_SIZE = 1000       # Larger to avoid spurious components (VOI)\n    MAX_COMPONENTS = 10             # Limit components to avoid over-segmentation (VOI)\n    FILL_HOLES_2D = True            # Fill holes slice-by-slice (TopoScore k=1,k=2)\n    FILL_HOLES_3D = False           # DISABLED - causes explosion to 100% (was True)\n    REMOVE_HANDLES = False          # Don't create handles (TopoScore k=1)\n    \n    # SURFACE SMOOTHING (Critical for SurfaceDice@τ=2.0)\n    SMOOTH_SURFACE = True\n    GAUSSIAN_SIGMA = 1.0            # Smooth for τ=2.0 tolerance\n    MORPHOLOGY_ITERATIONS = 2       # Smooth boundaries\n    \n    # Post-processing - TOPOLOGY-AWARE\n    CLOSING_RADIUS = 3              # Connect nearby surfaces (sheet continuity)\n    DILATION_RADIUS = 1             # Slight expansion for τ=2.0 tolerance\n    OPENING_RADIUS = 1              # Remove thin protrusions (topology noise)\n    \n    # Component analysis - AVOID SPLITS/MERGES (VOI)\n    KEEP_LARGEST_COMPONENT = False  # Keep multiple if significant\n    COMPONENT_SIZE_RATIO = 0.05     # Keep components > 5% of largest\n    \n    VERBOSE = True\n\nconfig = TopoConfig()\nconfig.OUTPUT_DIR.mkdir(exist_ok=True, parents=True)\n\nprint(\"\\nTopology-Aware Configuration:\")\nprint(f\"  Preserve topology: {config.PRESERVE_TOPOLOGY}\")\nprint(f\"  Surface smoothing: {config.SMOOTH_SURFACE}\")\nprint(f\"  Min component size: {config.MIN_COMPONENT_SIZE}\")\nprint(f\"  Soft voting threshold: {config.SOFT_VOTING_THRESHOLD}\")\nprint(\"=\"*60)\n\n# ============================================================================\n# TOPOLOGY-PRESERVING SURFACE DETECTOR\n# ============================================================================\n\nclass TopologyAwareSurfaceDetector:\n    \"\"\"Surface detector optimized for topology metrics\"\"\"\n    \n    def __init__(self, config):\n        self.config = config\n    \n    def normalize_volume(self, volume):\n        \"\"\"Robust normalization\"\"\"\n        volume = volume.astype(np.float32)\n        p1, p99 = np.percentile(volume, [1, 99])\n        volume = np.clip(volume, p1, p99)\n        \n        v_min, v_max = volume.min(), volume.max()\n        if v_max > v_min:\n            volume = (volume - v_min) / (v_max - v_min)\n        \n        return volume\n    \n    def detect_edges_multiscale(self, volume):\n        \"\"\"Multi-scale edge detection for surface boundaries\"\"\"\n        edges_list = []\n        \n        # Multiple scales for better surface detection\n        for sigma in [0.8, 1.2, 1.6]:\n            volume_smooth = gaussian(volume, sigma=sigma)\n            \n            # Sobel gradients\n            grad_z = sobel(volume_smooth, axis=0)\n            grad_y = sobel(volume_smooth, axis=1)\n            grad_x = sobel(volume_smooth, axis=2)\n            grad_mag = np.sqrt(grad_z**2 + grad_y**2 + grad_x**2)\n            \n            edges_list.append(grad_mag)\n        \n        # Take maximum across scales (catch edges at different levels)\n        edges = np.max(edges_list, axis=0)\n        \n        # Threshold\n        try:\n            thresh = threshold_otsu(edges[edges > 0])\n        except:\n            thresh = np.percentile(edges, config.GRADIENT_PERCENTILE)\n        \n        mask = edges > thresh * 0.65  # Moderate threshold\n        return mask\n    \n    def detect_intensity(self, volume):\n        \"\"\"Intensity-based detection\"\"\"\n        try:\n            thresh = threshold_otsu(volume)\n        except:\n            thresh = np.percentile(volume, config.INTENSITY_PERCENTILE)\n        \n        mask = volume > thresh * 0.88\n        return mask\n    \n    def detect_variance(self, volume):\n        \"\"\"Texture-based detection\"\"\"\n        local_mean = uniform_filter(volume, size=5)\n        local_mean_sq = uniform_filter(volume**2, size=5)\n        local_var = local_mean_sq - local_mean**2\n        \n        try:\n            thresh = threshold_otsu(local_var[local_var > 0])\n        except:\n            thresh = np.percentile(local_var, config.VARIANCE_PERCENTILE)\n        \n        mask = local_var > thresh * 0.75\n        return mask\n    \n    def detect_morphology(self, volume):\n        \"\"\"Morphological detection\"\"\"\n        tophat = ndimage.white_tophat(volume, structure=ball(2))\n        \n        try:\n            thresh = threshold_otsu(tophat[tophat > 0])\n        except:\n            thresh = np.percentile(tophat, config.MORPHOLOGY_PERCENTILE)\n        \n        mask = tophat > thresh * 0.85\n        return mask\n    \n    def ensemble_detection(self, volume):\n        \"\"\"Weighted ensemble detection\"\"\"\n        volume_norm = self.normalize_volume(volume)\n        self.current_volume = volume_norm  # Store for safety check\n        \n        print(\"  Multi-method detection...\")\n        \n        # Method 1: Edges (highest weight for SurfaceDice)\n        mask1 = self.detect_edges_multiscale(volume_norm)\n        print(f\"    Edges: {mask1.sum() / mask1.size:.2%}\")\n        \n        # Method 2: Intensity\n        mask2 = self.detect_intensity(volume_norm)\n        print(f\"    Intensity: {mask2.sum() / mask2.size:.2%}\")\n        \n        # Method 3: Variance\n        mask3 = self.detect_variance(volume_norm)\n        print(f\"    Variance: {mask3.sum() / mask3.size:.2%}\")\n        \n        # Method 4: Morphology\n        mask4 = self.detect_morphology(volume_norm)\n        print(f\"    Morphology: {mask4.sum() / mask4.size:.2%}\")\n        \n        # Weighted soft voting\n        combined = (mask1.astype(float) * config.GRADIENT_WEIGHT +\n                   mask2.astype(float) * config.INTENSITY_WEIGHT +\n                   mask3.astype(float) * config.VARIANCE_WEIGHT +\n                   mask4.astype(float) * config.MORPHOLOGY_WEIGHT)\n        \n        mask = combined > config.SOFT_VOTING_THRESHOLD\n        \n        ensemble_ratio = mask.sum() / mask.size\n        print(f\"  Ensemble: {ensemble_ratio:.2%}\")\n        \n        # Safety check on ensemble result\n        if ensemble_ratio > 0.80:\n            print(f\"    WARNING: Ensemble too high, using stricter threshold...\")\n            mask = combined > (config.SOFT_VOTING_THRESHOLD + 0.15)  # Raise threshold\n            print(f\"    Adjusted: {mask.sum() / mask.size:.2%}\")\n        \n        return mask\n    \n    def preserve_topology(self, mask):\n        \"\"\"Topology-preserving operations for TopoScore - FIXED\"\"\"\n        \n        initial_ratio = mask.sum() / mask.size\n        print(f\"    Before topology ops: {initial_ratio:.2%}\")\n        \n        # 1. Remove small components FIRST (before filling holes)\n        mask = remove_small_objects(mask, min_size=config.MIN_COMPONENT_SIZE, connectivity=3)\n        \n        after_small = mask.sum() / mask.size\n        print(f\"    After small removal: {after_small:.2%}\")\n        \n        # 2. Fill holes slice-by-slice (more conservative than 3D)\n        if config.FILL_HOLES_2D and after_small > 0.01 and after_small < 0.70:\n            filled_count = 0\n            for i in range(mask.shape[0]):\n                if mask[i].sum() > 0:  # Only process non-empty slices\n                    original_sum = mask[i].sum()\n                    mask[i] = ndimage.binary_fill_holes(mask[i])\n                    filled_count += (mask[i].sum() - original_sum)\n            \n            print(f\"    Filled {filled_count} voxels in 2D\")\n        \n        after_2d = mask.sum() / mask.size\n        print(f\"    After 2D fill: {after_2d:.2%}\")\n        \n        # 3. SKIP 3D hole filling if already high foreground\n        if config.FILL_HOLES_3D and after_2d > 0.01 and after_2d < 0.50:\n            before_3d = mask.sum()\n            mask = ndimage.binary_fill_holes(mask)\n            after_3d = mask.sum()\n            print(f\"    Filled {after_3d - before_3d} voxels in 3D\")\n        else:\n            print(f\"    Skipped 3D fill (fg={after_2d:.2%})\")\n        \n        # 4. Safety check - if exploded, revert to before filling\n        final_ratio = mask.sum() / mask.size\n        if final_ratio > 0.75:\n            print(f\"    WARNING: Exploded to {final_ratio:.2%}, reverting...\")\n            # Start over with stricter approach\n            mask_orig = self.ensemble_detection(self.current_volume)\n            mask = remove_small_objects(mask_orig, min_size=config.MIN_COMPONENT_SIZE * 2, connectivity=3)\n            mask = binary_opening(mask, footprint=ball(2))  # Remove weak connections\n        \n        # 5. Light smoothing only\n        if config.SMOOTH_SURFACE and final_ratio < 0.70:\n            mask = binary_closing(mask, footprint=ball(1))  # Reduced from ball(2)\n        \n        final_ratio = mask.sum() / mask.size\n        print(f\"    Final after topology: {final_ratio:.2%}\")\n        \n        return mask\n    \n    def component_filtering(self, mask):\n        \"\"\"Filter components to avoid VOI merge/split errors\"\"\"\n        \n        # Label connected components\n        labeled_mask, num_features = ndimage.label(mask, structure=np.ones((3, 3, 3)))\n        \n        if num_features == 0:\n            return mask\n        \n        print(f\"  Found {num_features} components\")\n        \n        # Compute sizes\n        sizes = ndimage.sum(mask, labeled_mask, range(1, num_features + 1))\n        sizes = np.array(sizes)\n        max_size = sizes.max()\n        \n        # Keep significant components (avoid over-segmentation)\n        keep_labels = []\n        for i, size in enumerate(sizes):\n            # Keep if: large enough absolute size OR significant relative to largest\n            if size > config.MIN_COMPONENT_SIZE or size > max_size * config.COMPONENT_SIZE_RATIO:\n                keep_labels.append(i + 1)\n        \n        # Limit total components to avoid VOI issues\n        if len(keep_labels) > config.MAX_COMPONENTS:\n            # Keep largest N components\n            top_indices = np.argsort(sizes)[-config.MAX_COMPONENTS:]\n            keep_labels = [i + 1 for i in top_indices]\n        \n        print(f\"  Keeping {len(keep_labels)} components\")\n        \n        mask_filtered = np.isin(labeled_mask, keep_labels)\n        return mask_filtered\n    \n    def smooth_surface(self, mask):\n        \"\"\"Smooth surface for better SurfaceDice@τ=2.0\"\"\"\n        \n        if not config.SMOOTH_SURFACE:\n            return mask\n        \n        # Morphological smoothing\n        # Closing to fill small gaps\n        mask = binary_closing(mask, footprint=ball(config.CLOSING_RADIUS))\n        \n        # Opening to remove small protrusions\n        mask = binary_opening(mask, footprint=ball(config.OPENING_RADIUS))\n        \n        # Light dilation to expand within τ=2.0 tolerance\n        mask = binary_dilation(mask, footprint=ball(config.DILATION_RADIUS))\n        \n        # Final closing for smooth surface\n        mask = binary_closing(mask, footprint=ball(2))\n        \n        return mask\n    \n    def adaptive_refinement(self, mask):\n        \"\"\"Adaptive refinement based on statistics\"\"\"\n        fg_ratio = mask.sum() / mask.size\n        \n        # EMERGENCY: If exploded, start over with very strict settings\n        if fg_ratio > 0.85:\n            print(f\"  EMERGENCY: {fg_ratio:.2%} - Rebuilding with strict settings...\")\n            # Use only gradient method with high threshold\n            volume_norm = self.current_volume\n            edges = self.detect_edges_multiscale(volume_norm)\n            thresh = np.percentile(edges, 65)  # Much higher threshold\n            mask = edges > thresh\n            mask = remove_small_objects(mask, min_size=2000, connectivity=3)\n            mask = binary_opening(mask, footprint=ball(2))\n            fg_ratio = mask.sum() / mask.size\n            print(f\"  Rebuilt: {fg_ratio:.2%}\")\n        \n        elif fg_ratio > 0.55:\n            print(f\"  High foreground ({fg_ratio:.2%}), applying strong erosion...\")\n            mask = binary_erosion(mask, footprint=ball(3))\n            mask = remove_small_objects(mask, min_size=2000, connectivity=3)\n            mask = binary_closing(mask, footprint=ball(1))\n        \n        elif fg_ratio < 0.12:\n            print(f\"  Low foreground ({fg_ratio:.2%}), expanding...\")\n            mask = binary_dilation(mask, footprint=ball(2))\n            mask = binary_closing(mask, footprint=ball(3))\n        \n        return mask\n    \n    def process_volume(self, volume):\n        \"\"\"Complete topology-aware pipeline\"\"\"\n        \n        # 1. Ensemble detection\n        mask = self.ensemble_detection(volume)\n        \n        # 2. Preserve topology\n        print(\"  Topology preservation...\")\n        mask = self.preserve_topology(mask)\n        \n        # 3. Component filtering (VOI optimization)\n        print(\"  Component filtering...\")\n        mask = self.component_filtering(mask)\n        \n        # 4. Surface smoothing (SurfaceDice optimization)\n        print(\"  Surface smoothing...\")\n        mask = self.smooth_surface(mask)\n        \n        # 5. Adaptive refinement\n        print(\"  Adaptive refinement...\")\n        mask = self.adaptive_refinement(mask)\n        \n        return mask\n\n# ============================================================================\n# VOLUME LOADER\n# ============================================================================\n\ndef load_volume(filepath):\n    \"\"\"Load 3D TIFF volume\"\"\"\n    try:\n        return tifffile.imread(str(filepath))\n    except:\n        from PIL import Image\n        with Image.open(str(filepath)) as img:\n            frames = []\n            try:\n                while True:\n                    frames.append(np.array(img))\n                    img.seek(img.tell() + 1)\n            except EOFError:\n                pass\n        return np.stack(frames, axis=0)\n\n# ============================================================================\n# MAIN PROCESSING\n# ============================================================================\n\ndef create_submission():\n    \"\"\"Create submission with topology-aware processing\"\"\"\n    \n    print(\"\\n\" + \"=\"*60)\n    print(\"STARTING TOPOLOGY-AWARE PROCESSING\")\n    print(\"=\"*60)\n    \n    detector = TopologyAwareSurfaceDetector(config)\n    \n    # Find test volumes\n    try:\n        test_df = pd.read_csv(config.ROOT_DIR / \"test.csv\")\n        test_files = [config.TEST_DIR / f\"{row['id']}.tif\" for _, row in test_df.iterrows()]\n        print(f\"✓ Loaded {len(test_files)} files\")\n    except:\n        test_files = list(config.TEST_DIR.glob(\"*.tif\"))\n        print(f\"✓ Found {len(test_files)} files\")\n    \n    if len(test_files) == 0:\n        print(\"✗ No files found!\")\n        with zipfile.ZipFile(config.ZIP_PATH, 'w') as zf:\n            pass\n        return\n    \n    # Process volumes\n    with zipfile.ZipFile(config.ZIP_PATH, 'w', zipfile.ZIP_DEFLATED, compresslevel=9) as zf:\n        for idx, filepath in enumerate(test_files, 1):\n            volume_id = filepath.stem\n            \n            print(f\"\\n{'='*60}\")\n            print(f\"Volume {idx}/{len(test_files)}: {volume_id}\")\n            print(f\"{'='*60}\")\n            \n            try:\n                # Load\n                volume = load_volume(filepath)\n                print(f\"  Shape: {volume.shape}, Range: [{volume.min()}, {volume.max()}]\")\n                \n                # Process\n                mask = detector.process_volume(volume)\n                \n                # Stats\n                fg_ratio = mask.sum() / mask.size\n                print(f\"\\n  Final foreground: {fg_ratio:.2%}\")\n                \n                # Expected score contributions\n                print(f\"  Expected contributions:\")\n                print(f\"    - SurfaceDice@2.0: Smooth surface ✓\")\n                print(f\"    - VOI: {len(np.unique(ndimage.label(mask)[0])) - 1} components\")\n                print(f\"    - TopoScore: Holes filled, topology preserved ✓\")\n                \n                # Save\n                mask_uint8 = mask.astype(np.uint8) * 255\n                temp_path = config.OUTPUT_DIR / f\"{volume_id}.tif\"\n                tifffile.imwrite(str(temp_path), mask_uint8, compression='zlib')\n                \n                with open(temp_path, 'rb') as f:\n                    zf.writestr(f\"{volume_id}.tif\", f.read())\n                \n                temp_path.unlink()\n                print(f\"  ✓ Success\")\n                \n            except Exception as e:\n                print(f\"  ✗ Error: {e}\")\n                import traceback\n                traceback.print_exc()\n    \n    print(f\"\\n{'='*60}\")\n    print(f\"✓✓✓ SUBMISSION CREATED ✓✓✓\")\n    print(f\"File: {config.ZIP_PATH}\")\n    print(f\"\\nOptimized for:\")\n    print(f\"  • 30% TopoScore (topology preservation)\")\n    print(f\"  • 35% SurfaceDice@τ=2.0 (surface smoothing)\")\n    print(f\"  • 35% VOI (component consistency)\")\n    print(f\"{'='*60}\")\n\n# ============================================================================\n# EXECUTE\n# ============================================================================\n\nif __name__ == \"__main__\":\n    try:\n        create_submission()\n    except Exception as e:\n        print(f\"Fatal error: {e}\")\n        import traceback\n        traceback.print_exc()\n        with zipfile.ZipFile(config.ZIP_PATH, 'w') as zf:\n            pass","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-07T20:43:59.065853Z","iopub.execute_input":"2025-12-07T20:43:59.067075Z","iopub.status.idle":"2025-12-07T20:45:30.974999Z","shell.execute_reply.started":"2025-12-07T20:43:59.067040Z","shell.execute_reply":"2025-12-07T20:45:30.974080Z"}},"outputs":[],"execution_count":null}]}