{"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":"nvidiaTeslaT4","dataSources":[{"sourceType":"competition","sourceId":99552,"databundleVersionId":13441085},{"sourceType":"kernelVersion","sourceId":255517047},{"sourceType":"kernelVersion","sourceId":257435185}],"dockerImageVersionId":31090,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## References\n- https://www.kaggle.com/code/ahsuna123/voxel-by-voxel-3d-cnn-intracranial-aneurysms  \n\n## My Notebooks\n- Train: https://www.kaggle.com/code/ichigoe/voxel-by-voxel-3d-cnn-train/notebook\n- Inference: here","metadata":{}},{"cell_type":"code","source":"# =========================\n# Imports & Device Settings\n# =========================\nimport os\nimport gc\nimport shutil\nfrom collections import OrderedDict\nfrom typing import Tuple, List\n\nimport numpy as np\nimport polars as pl\nimport pydicom\nfrom scipy import ndimage\n\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nfrom torch import amp\n\nimport kaggle_evaluation.rsna_inference_server\n\n# Device configuration\nDEVICE = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\nUSE_AMP = torch.cuda.is_available()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-30T03:05:40.241718Z","iopub.execute_input":"2025-08-30T03:05:40.242355Z","iopub.status.idle":"2025-08-30T03:05:47.447437Z","shell.execute_reply.started":"2025-08-30T03:05:40.242327Z","shell.execute_reply":"2025-08-30T03:05:47.446827Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =====================\n# Competition Constants\n# =====================\nID_COL = 'SeriesInstanceUID'\nLABEL_COLS = [\n    'Left Infraclinoid Internal Carotid Artery',\n    'Right Infraclinoid Internal Carotid Artery',\n    'Left Supraclinoid Internal Carotid Artery',\n    'Right Supraclinoid Internal Carotid Artery',\n    'Left Middle Cerebral Artery',\n    'Right Middle Cerebral Artery',\n    'Anterior Communicating Artery',\n    'Left Anterior Cerebral Artery',\n    'Right Anterior Cerebral Artery',\n    'Left Posterior Communicating Artery',\n    'Right Posterior Communicating Artery',\n    'Basilar Tip',\n    'Other Posterior Circulation',\n    'Aneurysm Present',\n]\n\n# =========\n# File paths\n# =========\nMODEL_WEIGHTS_PATH = \"/kaggle/input/voxel-by-voxel-3d-cnn-train/model_weights.pth\"\n\n# ==========================\n# Processing configurations\n# ==========================\nTARGET_SIZE = (64, 64, 64)      # final volume dimensions (D,H,W)\nTARGET_SPACING_MM = 1.0         # isotropic spacing in millimeters\nCTA_WINDOW = (300.0, 700.0)     # CT windowing (center, width)\nMRI_Z_CLIP = 3.0                # z-score clipping for MRI\nLRU_CAPACITY = 8                # memory cache capacity","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-30T03:05:47.448271Z","iopub.execute_input":"2025-08-30T03:05:47.448693Z","iopub.status.idle":"2025-08-30T03:05:47.453389Z","shell.execute_reply.started":"2025-08-30T03:05:47.448664Z","shell.execute_reply":"2025-08-30T03:05:47.45283Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ==========================\n# Volume resizing utilities\n# ==========================\ndef _safe_zoom(volume: np.ndarray, zoom_factors: Tuple[float, ...], order: int = 1) -> np.ndarray:\n    \"\"\"Apply ndimage zoom with protection against invalid factors and shapes.\"\"\"\n    volume = np.nan_to_num(volume, copy=False)\n    zf = tuple(float(max(1e-6, f)) for f in zoom_factors)\n    if len(zf) != volume.ndim:\n        if len(zf) > volume.ndim:\n            zf = zf[:volume.ndim]\n        else:\n            zf = (1.0,) * (volume.ndim - len(zf)) + zf\n    return ndimage.zoom(volume, zf, order=order)\n\ndef _resize_slice(arr: np.ndarray, out_h: int, out_w: int) -> np.ndarray:\n    \"\"\"Resize a 2D slice to specified dimensions using bilinear interpolation.\"\"\"\n    h, w = arr.shape\n    if h == out_h and w == out_w:\n        return arr.astype(np.float32, copy=False)\n    zy = out_h / max(h, 1)\n    zx = out_w / max(w, 1)\n    return _safe_zoom(arr, (zy, zx), order=1).astype(np.float32, copy=False)\n\n# ==========================\n# DICOM series processor\n# ==========================\nclass DICOMProcessor:\n    \"\"\"Process DICOM series into normalized 3D volumes with LRU caching.\"\"\"\n\n    def __init__(\n        self,\n        target_size: Tuple[int, int, int] = TARGET_SIZE,\n        target_spacing_mm: float = TARGET_SPACING_MM,\n        cta_window: Tuple[float, float] = CTA_WINDOW,\n        mri_z_clip: float = MRI_Z_CLIP,\n        lru_capacity: int = LRU_CAPACITY,\n    ):\n        self.target_size = target_size\n        self.target_spacing_mm = target_spacing_mm\n        self.cta_window = cta_window\n        self.mri_z_clip = mri_z_clip\n        self.memory_cache = OrderedDict()\n        self.lru_capacity = lru_capacity\n        \n        # Processing adjustment counters\n        self.slope_adjustments = 0\n        self.intercept_adjustments = 0\n        self.adaptive_windowing_count = 0\n\n    def _cache_put(self, key: str, vol: np.ndarray):\n        \"\"\"Store volume in LRU cache.\"\"\"\n        self.memory_cache[key] = vol\n        self.memory_cache.move_to_end(key)\n        if len(self.memory_cache) > self.lru_capacity:\n            self.memory_cache.popitem(last=False)\n\n    def _cache_get(self, key: str):\n        \"\"\"Retrieve volume from LRU cache.\"\"\"\n        if key in self.memory_cache:\n            vol = self.memory_cache[key]\n            self.memory_cache.move_to_end(key)\n            return vol\n        return None\n\n    def _validate_and_apply_rescale(self, sl: np.ndarray, ds) -> np.ndarray:\n        \"\"\"Validate slope/intercept values and apply rescaling with robust error handling.\"\"\"\n        slope = float(getattr(ds, \"RescaleSlope\", 1.0))\n        intercept = float(getattr(ds, \"RescaleIntercept\", 0.0))\n        \n        # Validate slope\n        if slope <= 0 or not np.isfinite(slope) or abs(slope) > 1000:\n            slope = 1.0\n            self.slope_adjustments += 1\n        \n        # Validate intercept with range clamping for extreme values\n        if not np.isfinite(intercept):\n            intercept = 0.0\n            self.intercept_adjustments += 1\n        elif abs(intercept) > 10000:\n            intercept = np.clip(intercept, -2000, 0)\n            self.intercept_adjustments += 1\n        \n        # Apply rescaling\n        rescaled = sl * slope + intercept\n        \n        # Validate result\n        if np.any(~np.isfinite(rescaled)):\n            rescaled = np.nan_to_num(rescaled, copy=False)\n        \n        # Post-rescale range check\n        min_val, max_val = rescaled.min(), rescaled.max()\n        if min_val < -5000 or max_val > 10000:\n            rescaled = np.clip(rescaled, -3000, 5000)\n        \n        return rescaled\n\n    def load_dicom_series(self, series_path: str) -> np.ndarray:\n        \"\"\"Load and process DICOM series into normalized 3D volume.\"\"\"\n        series_id = os.path.basename(series_path)\n\n        # Check memory cache first\n        cached_volume = self._cache_get(series_id)\n        if cached_volume is not None and isinstance(cached_volume, np.ndarray) and cached_volume.shape == self.target_size:\n            return cached_volume\n\n        try:\n            # Collect DICOM datasets with pixel data\n            dicoms = []\n            for root, _, files in os.walk(series_path):\n                for f in files:\n                    if f.endswith(\".dcm\"):\n                        try:\n                            ds = pydicom.dcmread(os.path.join(root, f), force=True)\n                            if hasattr(ds, \"PixelData\"):\n                                dicoms.append(ds)\n                        except Exception:\n                            continue\n            \n            if not dicoms:\n                raise ValueError(f\"No valid DICOM files with pixel data in {series_path}\")\n\n            # Sort slices spatially\n            dicoms = self._sort_slices(dicoms)\n\n            # Determine if multiframe series\n            has_multiframe = any(getattr(ds, \"NumberOfFrames\", 1) > 1 for ds in dicoms)\n\n            # Calculate spatial spacing\n            spacing = self._get_spacing(dicoms, has_multiframe=has_multiframe)\n\n            # Determine consistent base shape\n            base_h, base_w = self._choose_base_shape(dicoms)\n\n            # Extract modality\n            modality_tag = (getattr(dicoms[0], \"Modality\", \"\") or \"\").upper()\n\n            # Process pixel data from all datasets\n            vol_slices = []\n            for ds in dicoms:\n                arr = ds.pixel_array\n                # Standardize to (N,H,W) where N=number of frames (1 if 2D)\n                if arr.ndim >= 3:\n                    h, w = arr.shape[-2], arr.shape[-1]\n                    n = int(np.prod(arr.shape[:-2]))\n                    arr = arr.reshape(n, h, w)\n                    frames = arr\n                else:\n                    frames = arr[np.newaxis, ...]  # shape (1,H,W)\n\n                for sl in frames:\n                    sl = sl.astype(np.float32)\n\n                    # Handle inverted grayscale\n                    if getattr(ds, \"PhotometricInterpretation\", \"MONOCHROME2\") == \"MONOCHROME1\":\n                        sl = sl.max() - sl\n\n                    # Apply validated rescaling\n                    sl = self._validate_and_apply_rescale(sl, ds)\n\n                    # Standardize slice dimensions\n                    sl = _resize_slice(sl, base_h, base_w)\n                    vol_slices.append(sl)\n\n            if len(vol_slices) == 0:\n                raise ValueError(\"No valid slices extracted.\")\n\n            # Create 3D volume\n            volume = np.stack(vol_slices, axis=0).astype(np.float32)  # (D,H,W)\n\n            # Apply modality-specific normalization\n            volume = self._normalize_by_modality(volume, modality_tag)\n\n            # Resample to isotropic spacing\n            if self.target_spacing_mm is not None:\n                dz, dy, dx = spacing\n                z, y, x = volume.shape\n                newD = max(1, int(round(z * dz / self.target_spacing_mm)))\n                newH = max(1, int(round(y * dy / self.target_spacing_mm)))\n                newW = max(1, int(round(x * dx / self.target_spacing_mm)))\n                volume = _safe_zoom(volume, (newD / z, newH / y, newW / x), order=1)\n\n            # Resize to target dimensions\n            tz, ty, tx = self.target_size\n            z, y, x = volume.shape\n            volume = _safe_zoom(volume, (tz / z, ty / y, tx / x), order=1).astype(np.float32)\n\n            # Cache and return\n            self._cache_put(series_id, volume)\n            return volume\n\n        except Exception:\n            vol = np.zeros(self.target_size, dtype=np.float32)\n            self._cache_put(series_id, vol)\n            return vol\n\n    def _sort_slices(self, ds_list: List[pydicom.dataset.FileDataset]) -> List[pydicom.dataset.FileDataset]:\n        \"\"\"Sort DICOM slices by spatial position along normal vector.\"\"\"\n        try:\n            orient = np.array(ds_list[0].ImageOrientationPatient, dtype=np.float32)\n            row = orient[:3]\n            col = orient[3:]\n            normal = np.cross(row, col)\n            \n            def sort_key(ds):\n                ipp = np.array(getattr(ds, \"ImagePositionPatient\", [0, 0, 0]), dtype=np.float32)\n                return float(np.dot(ipp, normal))\n            \n            return sorted(ds_list, key=sort_key)\n        except Exception:\n            return sorted(ds_list, key=lambda ds: getattr(ds, \"InstanceNumber\", 0))\n\n    def _get_spacing(self, ds_sorted: List[pydicom.dataset.FileDataset], has_multiframe: bool = False) -> Tuple[float, float, float]:\n        \"\"\"Calculate spatial spacing in millimeters (dz, dy, dx).\"\"\"\n        try:\n            dy, dx = map(float, ds_sorted[0].PixelSpacing)\n        except Exception:\n            ps = getattr(ds_sorted[0], \"PixelSpacing\", [1.0, 1.0])\n            dy, dx = float(ps[0]), float(ps[1])\n\n        if has_multiframe:\n            dz = float(getattr(ds_sorted[0], \"SpacingBetweenSlices\", getattr(ds_sorted[0], \"SliceThickness\", 1.0)))\n        else:\n            zs = []\n            for i in range(1, len(ds_sorted)):\n                p0 = np.array(getattr(ds_sorted[i-1], \"ImagePositionPatient\", [0, 0, 0]), dtype=np.float32)\n                p1 = np.array(getattr(ds_sorted[i], \"ImagePositionPatient\", [0, 0, 0]), dtype=np.float32)\n                d = np.linalg.norm(p1 - p0)\n                if d > 0:\n                    zs.append(d)\n            if zs:\n                dz = float(np.median(zs))\n            else:\n                dz = float(getattr(ds_sorted[0], \"SliceThickness\", 1.0))\n\n        # Validate spacing values\n        dz = dz if (dz > 0 and np.isfinite(dz)) else 1.0\n        dy = dy if (dy > 0 and np.isfinite(dy)) else 1.0\n        dx = dx if (dx > 0 and np.isfinite(dx)) else 1.0\n        return (dz, dy, dx)\n\n    def _choose_base_shape(self, ds_list: List[pydicom.dataset.FileDataset]) -> Tuple[int, int]:\n        \"\"\"Select most common image dimensions as base shape.\"\"\"\n        shapes = []\n        for ds in ds_list:\n            try:\n                h, w = int(ds.Rows), int(ds.Columns)\n            except Exception:\n                arr = ds.pixel_array\n                h, w = arr.shape[-2], arr.shape[-1]\n            shapes.append((h, w))\n        vals, counts = np.unique(shapes, return_counts=True, axis=0)\n        base = tuple(vals[counts.argmax()])\n        return int(base[0]), int(base[1])\n\n    def _normalize_by_modality(self, volume: np.ndarray, modality_tag: str) -> np.ndarray:\n        \"\"\"Apply modality-specific normalization with adaptive windowing for CT.\"\"\"\n        volume = np.nan_to_num(volume, copy=False)\n        \n        if modality_tag == \"CT\":\n            min_val, max_val = volume.min(), volume.max()\n            \n            # Check if values are in normal CT range\n            if min_val >= -2000 and max_val <= 4000:\n                # Normal range: use standard windowing\n                c, w = self.cta_window\n                lo, hi = c - w / 2.0, c + w / 2.0\n            else:\n                # Extreme range: use adaptive windowing\n                self.adaptive_windowing_count += 1\n                \n                # Percentile-based adaptive window\n                p1, p99 = np.percentile(volume, [1, 99])\n                margin = (p99 - p1) * 0.1\n                lo = p1 - margin\n                hi = p99 + margin\n                \n                # Ensure minimum window width\n                if hi - lo < 100:\n                    center = (hi + lo) / 2\n                    lo = center - 50\n                    hi = center + 50\n            \n            v = np.clip(volume, lo, hi)\n            v = (v - lo) / (hi - lo + 1e-6)\n            return v.astype(np.float32, copy=False)\n        else:\n            # MRI processing with robust statistics validation\n            mean = float(volume.mean())\n            std = float(volume.std() + 1e-6)\n            \n            # Validate statistics\n            if std < 1e-6 or not np.isfinite(mean) or not np.isfinite(std):\n                return np.full_like(volume, 0.5, dtype=np.float32)\n            \n            # Check dynamic range\n            min_val, max_val = volume.min(), volume.max()\n            if max_val - min_val < 1e-6:\n                return np.full_like(volume, 0.5, dtype=np.float32)\n            \n            v = (volume - mean) / std\n            zc = float(self.mri_z_clip)\n            v = np.clip(v, -zc, zc)\n            v = (v + zc) / (2.0 * zc)\n            return v.astype(np.float32, copy=False)\n\n# =======================\n# 3D CNN model\n# =======================\nclass Simple3DCNN(nn.Module):\n    \"\"\"Lightweight 3D CNN for multi-label aneurysm classification.\"\"\"\n    \n    def __init__(self, num_classes: int = len(LABEL_COLS)):\n        super(Simple3DCNN, self).__init__()\n        self.conv1 = nn.Conv3d(1, 16, kernel_size=3, padding=1)\n        self.bn1 = nn.BatchNorm3d(16)\n        self.pool1 = nn.MaxPool3d(2)\n\n        self.conv2 = nn.Conv3d(16, 32, kernel_size=3, padding=1)\n        self.bn2 = nn.BatchNorm3d(32)\n        self.pool2 = nn.MaxPool3d(2)\n\n        self.conv3 = nn.Conv3d(32, 64, kernel_size=3, padding=1)\n        self.bn3 = nn.BatchNorm3d(64)\n        self.pool3 = nn.MaxPool3d(2)\n\n        self.conv4 = nn.Conv3d(64, 128, kernel_size=3, padding=1)\n        self.bn4 = nn.BatchNorm3d(128)\n        self.pool4 = nn.MaxPool3d(2)\n\n        self.adaptive_pool = nn.AdaptiveAvgPool3d((2, 2, 2))\n        self.fc1 = nn.Linear(128 * 2 * 2 * 2, 256)\n        self.dropout1 = nn.Dropout(0.5)\n        self.fc2 = nn.Linear(256, 128)\n        self.dropout2 = nn.Dropout(0.3)\n        self.fc3 = nn.Linear(128, num_classes)\n\n    def forward(self, x):\n        \"\"\"Forward pass through 3D CNN.\"\"\"\n        x = self.pool1(F.relu(self.bn1(self.conv1(x))))\n        x = self.pool2(F.relu(self.bn2(self.conv2(x))))\n        x = self.pool3(F.relu(self.bn3(self.conv3(x))))\n        x = self.pool4(F.relu(self.bn4(self.conv4(x))))\n        x = self.adaptive_pool(x)\n        x = x.view(x.size(0), -1)\n        x = F.relu(self.fc1(x))\n        x = self.dropout1(x)\n        x = F.relu(self.fc2(x))\n        x = self.dropout2(x)\n        x = self.fc3(x)  # return logits\n        return x","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-30T03:05:47.455092Z","iopub.execute_input":"2025-08-30T03:05:47.455305Z","iopub.status.idle":"2025-08-30T03:05:47.49204Z","shell.execute_reply.started":"2025-08-30T03:05:47.455287Z","shell.execute_reply":"2025-08-30T03:05:47.491247Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =====================================\n# Initialize processor and load model\n# =====================================\n# Create DICOM processor with LRU caching\nprocessor = DICOMProcessor(\n    target_size=TARGET_SIZE,\n    target_spacing_mm=TARGET_SPACING_MM,\n    cta_window=CTA_WINDOW,\n    mri_z_clip=MRI_Z_CLIP,\n    lru_capacity=LRU_CAPACITY,\n)\n\n# Initialize model and load trained weights\nmodel = Simple3DCNN(num_classes=len(LABEL_COLS)).to(DEVICE)\ntry:\n    state = torch.load(MODEL_WEIGHTS_PATH, map_location='cpu')\n    model.load_state_dict(state, strict=True)\n    print(f\"Loaded weights from {MODEL_WEIGHTS_PATH}\")\nexcept Exception as e:\n    print(f\"[Weight loading warning] {e}\\nProceeding with randomly initialized weights.\")\nmodel.eval()\n\n# ==============\n# Inference API\n# ==============\n@torch.no_grad()\ndef predict(series_path: str) -> pl.DataFrame:\n    \"\"\"Main prediction function called by evaluation server.\"\"\"\n    try:\n        # Process DICOM series to normalized volume\n        volume = processor.load_dicom_series(series_path)  # (D,H,W) in [0,1]\n        volume_tensor = torch.from_numpy(volume).unsqueeze(0).unsqueeze(0).to(DEVICE)  # (1,1,D,H,W)\n\n        # Run inference with mixed precision\n        with amp.autocast(device_type='cuda', enabled=USE_AMP):\n            logits = model(volume_tensor)               # (1,14)\n            probs = torch.sigmoid(logits).cpu().numpy().flatten().tolist()\n\n        # Format results as DataFrame\n        result_df = pl.DataFrame(data=[probs], schema=LABEL_COLS, orient='row')\n\n        # Clean up GPU memory\n        del volume_tensor\n        if torch.cuda.is_available():\n            torch.cuda.empty_cache()\n        gc.collect()\n\n    except Exception as e:\n        print(f\"[Prediction error] {e}\")\n        # Return neutral probabilities on failure\n        result_df = pl.DataFrame(data=[[0.5] * len(LABEL_COLS)], schema=LABEL_COLS, orient='row')\n\n    # Clean up temporary files\n    shutil.rmtree('/kaggle/shared', ignore_errors=True)\n    return result_df","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-30T03:05:47.492741Z","iopub.execute_input":"2025-08-30T03:05:47.49298Z","iopub.status.idle":"2025-08-30T03:05:47.877824Z","shell.execute_reply.started":"2025-08-30T03:05:47.492961Z","shell.execute_reply":"2025-08-30T03:05:47.877175Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ==========================\n# Run evaluation server\n# ==========================\ninference_server = kaggle_evaluation.rsna_inference_server.RSNAInferenceServer(predict)\n\nif os.getenv('KAGGLE_IS_COMPETITION_RERUN'):\n    inference_server.serve()\nelse:\n    inference_server.run_local_gateway()\n    try:\n        sub = pl.read_parquet('/kaggle/working/submission.parquet')\n        print(sub.head())\n    except Exception as e:\n        print(f\"Submission parquet not found: {e}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-30T03:05:47.878501Z","iopub.execute_input":"2025-08-30T03:05:47.878761Z","iopub.status.idle":"2025-08-30T03:06:14.11064Z","shell.execute_reply.started":"2025-08-30T03:05:47.878741Z","shell.execute_reply":"2025-08-30T03:06:14.109819Z"}},"outputs":[],"execution_count":null}]}