{"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":"gpu","dataSources":[{"sourceType":"competition","sourceId":99552,"databundleVersionId":13441085}],"dockerImageVersionId":31090,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"!pip install -q monai timm pydicom nibabel torchio","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-09T07:09:07.937671Z","iopub.execute_input":"2025-09-09T07:09:07.937987Z","iopub.status.idle":"2025-09-09T07:09:11.379774Z","shell.execute_reply.started":"2025-09-09T07:09:07.937953Z","shell.execute_reply":"2025-09-09T07:09:11.378771Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nimport multiprocessing as mp\nimport numpy as np\nimport pandas as pd\nimport torch\nimport torch.nn as nn\nfrom torch.utils.data import Dataset, DataLoader\nimport monai\nfrom monai.data import CacheDataset, DataLoader as MonaiDataLoader\nfrom monai.transforms import (\n    EnsureChannelFirstd, Resized, NormalizeIntensityd, RandFlipd,\n    RandRotated, Compose, ToTensord\n)\nimport pydicom\nfrom pathlib import Path\nfrom tqdm import tqdm\nfrom functools import partial\nfrom sklearn.metrics import roc_auc_score\nimport SimpleITK as sitk\nimport h5py\nimport torchio as tio\nfrom collections import OrderedDict\nfrom typing import Tuple, List\nfrom scipy import ndimage\nimport ast\nfrom sklearn.model_selection import train_test_split\nimport nibabel as nb\nimport shutil\n\n\n# Config\nDATA_DIR = \"/kaggle/input/rsna-intracranial-aneurysm-detection/series/\"\nTRAIN_CSV = \"/kaggle/input/rsna-intracranial-aneurysm-detection/train.csv\"\nLOCALIZERS_CSV = \"/kaggle/input/rsna-intracranial-aneurysm-detection/train_localizers.csv\"\nSEGMENTATIONS_DIR = \"/kaggle/input/rsna-intracranial-aneurysm-detection/segmentations/\"\nOUTPUT_DIR = \"/kaggle/working/preprocessed/\"\nBATCH_SIZE = 8\nNUM_EPOCHS = 50\nLR = 1e-4\nDEVICE = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\nTARGET_SIZE = (128, 128, 64)  # Downsample resolution\nTARGET_SPACING_MM = 1.0  # Isotropic spacing\nCTA_WINDOW = (-1000, 1000)  # HU window for CTA\nMRI_Z_CLIP = 3.0  # Z-score clip for MRI\nLRU_CAPACITY = 10  # Memory cache capacity\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# Utility Functions\ndef _safe_zoom(volume: np.ndarray, zoom_factors: Tuple[float, ...], order: int = 1) -> np.ndarray:\n    \"\"\"Robust wrapper around ndimage.zoom to avoid rank mismatch and invalid factors.\"\"\"\n    volume = np.nan_to_num(volume, copy=False)\n    zf = tuple(float(max(1e-6, f)) for f in zoom_factors)  # avoid zeros/negatives\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 (out_h, out_w) using safe zoom.\"\"\"\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# DICOM Processor Class\nclass DICOMProcessor:\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    def _cache_put(self, key: str, vol: np.ndarray):\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        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 _sort_slices(self, ds_list: List[pydicom.dataset.FileDataset]) -> List[pydicom.dataset.FileDataset]:\n        try:\n            orient = np.array(ds_list[0].ImageOrientationPatient, dtype=np.float32)\n            row, col = orient[:3], orient[3:]\n            normal = np.cross(row, col)\n            return sorted(ds_list, key=lambda ds: float(np.dot(np.array(ds.ImagePositionPatient, dtype=np.float32), normal)))\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        try:\n            dy, dx = map(float, ds_sorted[0].PixelSpacing)\n        except Exception:\n            dy, dx = 1.0, 1.0\n        if has_multiframe:\n            dz = float(getattr(ds_sorted[0], \"SpacingBetweenSlices\", getattr(ds_sorted[0], \"SliceThickness\", 1.0)))\n        else:\n            zs = [np.linalg.norm(np.array(ds_sorted[i].ImagePositionPatient, dtype=np.float32) - np.array(ds_sorted[i-1].ImagePositionPatient, dtype=np.float32)) for i in range(1, len(ds_sorted))]\n            dz = np.median([z for z in zs if z > 0]) if zs else float(getattr(ds_sorted[0], \"SliceThickness\", 1.0))\n        return (dz if dz > 0 else 1.0, dy if dy > 0 else 1.0, dx if dx > 0 else 1.0)\n\n    def _choose_base_shape(self, ds_list: List[pydicom.dataset.FileDataset]) -> Tuple[int, int]:\n        shapes = [(int(ds.Rows), int(ds.Columns)) for ds in ds_list if hasattr(ds, \"Rows\") and hasattr(ds, \"Columns\")]\n        if not shapes:\n            shapes = [(ds.pixel_array.shape[-2], ds.pixel_array.shape[-1]) for ds in ds_list if hasattr(ds, \"pixel_array\")]\n        vals, counts = np.unique(shapes, return_counts=True, axis=0)\n        return tuple(int(x) for x in vals[counts.argmax()]) if vals.size else (512, 512)\n\n    def _normalize_by_modality(self, volume: np.ndarray, modality_tag: str) -> np.ndarray:\n        volume = np.nan_to_num(volume, copy=False)\n        if modality_tag == \"CT\":\n            c, w = self.cta_window\n            lo, hi = c - w / 2.0, c + w / 2.0\n            return np.clip((volume - lo) / (hi - lo + 1e-6), 0, 1).astype(np.float32)\n        else:\n            mean = float(volume.mean())\n            std = float(volume.std() + 1e-6)\n            v = np.clip((volume - mean) / std, -self.mri_z_clip, self.mri_z_clip)\n            return ((v + self.mri_z_clip) / (2.0 * self.mri_z_clip)).astype(np.float32)\n\n    def load_dicom_series(self, series_path: str) -> np.ndarray:\n        series_id = os.path.basename(series_path)\n        m = self._cache_get(series_id)\n        if m is not None and m.shape == self.target_size:\n            return m\n\n        try:\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 as e:\n                            print(f\"[DICOM read] {series_id}: {e}\")\n                            continue\n            if not dicoms:\n                raise ValueError(f\"No valid DICOM files in {series_path}\")\n\n            dicoms = self._sort_slices(dicoms)\n            has_multiframe = any(getattr(ds, \"NumberOfFrames\", 1) > 1 for ds in dicoms)\n            spacing = self._get_spacing(dicoms, has_multiframe)\n            base_h, base_w = self._choose_base_shape(dicoms)\n            modality_tag = getattr(dicoms[0], \"Modality\", \"\").upper()\n\n            vol_slices = []\n            for ds in dicoms:\n                arr = ds.pixel_array\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                else:\n                    arr = arr[np.newaxis, ...]\n                for sl in arr:\n                    sl = sl.astype(np.float32)\n                    if getattr(ds, \"PhotometricInterpretation\", \"MONOCHROME2\") == \"MONOCHROME1\":\n                        sl = sl.max() - sl\n                    slope = float(getattr(ds, \"RescaleSlope\", 1.0))\n                    intercept = float(getattr(ds, \"RescaleIntercept\", 0.0))\n                    sl = sl * slope + intercept\n                    sl = _resize_slice(sl, base_h, base_w)\n                    vol_slices.append(sl)\n\n            if not vol_slices:\n                raise ValueError(f\"No valid slices in {series_id}\")\n            volume = np.stack(vol_slices, axis=0)\n\n            # Resample to target spacing\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 grid\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            volume = self._normalize_by_modality(volume, modality_tag)\n            self._cache_put(series_id, volume)\n            return volume\n\n        except Exception as e:\n            print(f\"[Processor] Error processing {series_id}: {e}\")\n            vol = np.zeros(self.target_size, dtype=np.float32)\n            self._cache_put(series_id, vol)\n            return vol\n            \ndef process_row(args):\n    row_dict, data_dir, label_cols, processor = args\n    series_id = row_dict[\"SeriesInstanceUID\"]\n    series_path = os.path.join(data_dir, series_id)\n    labels = [row_dict.get(col, 0) for col in label_cols]\n    volume = processor.load_dicom_series(series_path)\n    return series_path, volume, np.array(labels, dtype=np.float32)\n\ndef preprocess_dataset():\n    if os.path.exists(OUTPUT_DIR):\n        shutil.rmtree(OUTPUT_DIR)\n    os.makedirs(OUTPUT_DIR, exist_ok=True)\n\n    train_df = pd.read_csv(TRAIN_CSV)\n    # train_df = train_df.sample(n=50, random_state=42)  \n    print(f\"Preprocessing {len(train_df)} series (with both positive & negative cases)\")\n\n    localizers_df = pd.read_csv(LOCALIZERS_CSV)\n    def parse_coordinates(coord_str):\n        try:\n            d = ast.literal_eval(coord_str)\n            return pd.Series({\"center_x\": d.get(\"x\", np.nan), \"center_y\": d.get(\"y\", np.nan), \"center_z\": d.get(\"z\", np.nan) if \"z\" in d else np.nan})\n        except:\n            return pd.Series({\"center_x\": np.nan, \"center_y\": np.nan, \"center_z\": np.nan})\n    coord_df = localizers_df[\"coordinates\"].apply(parse_coordinates)\n    localizers_df = pd.concat([localizers_df.drop(columns=[\"coordinates\"]), coord_df], axis=1)\n    localizers_df = localizers_df.groupby(\"SeriesInstanceUID\")[[\"center_x\", \"center_y\", \"center_z\"]].mean().reset_index()\n\n    train_df = train_df.merge(localizers_df, on=\"SeriesInstanceUID\", how=\"left\")\n    for col in [\"center_x\", \"center_y\", \"center_z\"]:\n        train_df[col] = train_df[col].fillna(0)\n\n    print(\"Label distribution (Aneurysm Present):\")\n    print(train_df[\"Aneurysm Present\"].value_counts())\n\n    processor = DICOMProcessor(target_size=TARGET_SIZE, target_spacing_mm=TARGET_SPACING_MM, cta_window=CTA_WINDOW, mri_z_clip=MRI_Z_CLIP)\n    args = [(row.to_dict(), DATA_DIR, LABEL_COLS, processor) for _, row in train_df.iterrows()]\n\n    series_data = []\n    results = []\n    for arg in tqdm(args, desc=\"Preprocessing\"):\n        results.append(process_row(arg))\n\n    for series_path, volume, labels in results:\n        if volume is not None and not np.all(volume == 0):\n            series_id = os.path.basename(series_path)\n            np.save(os.path.join(OUTPUT_DIR, f\"{series_id}.npy\"), volume)\n            # Removed segmentation saving\n            series_data.append({\n                \"image\": os.path.join(OUTPUT_DIR, f\"{series_id}.npy\"),\n                \"labels\": labels.tolist()\n            })\n\n    pd.DataFrame(series_data).to_csv(\"series_metadata.csv\", index=False)\n    print(\"Preprocessing done!\")\n\nclass AneurysmDataset(Dataset):\n    def __init__(self, data: pd.DataFrame, transform=None):\n        self.data = data\n        self.transform = transform\n        self.segmentations_dir = SEGMENTATIONS_DIR\n\n    def __len__(self):\n        return len(self.data)\n\n    def __getitem__(self, idx):\n        row = self.data.iloc[idx]\n        image = np.load(row[\"image\"]).astype(np.float32)\n        if image.ndim == 3:\n            image = np.expand_dims(image, axis=0)\n        \n        # Load and preprocess segmentation mask on-the-fly\n        series_id = os.path.basename(row[\"image\"]).replace(\".npy\", \"\")\n        seg_path = os.path.join(self.segmentations_dir, f\"{series_id}.nii.gz\")\n        if os.path.exists(seg_path):\n            try:\n                seg_img = nb.load(seg_path)\n                seg_data = seg_img.get_fdata()\n                tz, ty, tx = TARGET_SIZE\n                z, y, x = seg_data.shape\n                seg_data = _safe_zoom(seg_data, (tz / z, ty / y, tx / x), order=0).astype(np.float32)\n                seg_data = (seg_data > 0).astype(np.float32)\n            except Exception as e:\n                print(f\"Error loading segmentation for {series_id}: {e}\")\n                seg_data = np.zeros(TARGET_SIZE, dtype=np.float32)\n        else:\n            seg_data = np.zeros(TARGET_SIZE, dtype=np.float32)\n        \n        if seg_data.ndim == 3:\n            seg_data = np.expand_dims(seg_data, axis=0)\n        \n        labels = np.array(ast.literal_eval(row[\"labels\"]), dtype=np.float32) if isinstance(row[\"labels\"], str) else np.array(row[\"labels\"], dtype=np.float32)\n        \n        if self.transform:\n            subject = tio.Subject(image=tio.ScalarImage(tensor=image), segmentation=tio.LabelMap(tensor=seg_data))\n            transformed = self.transform(subject)\n            image = transformed[\"image\"].numpy()\n            seg_data = transformed[\"segmentation\"].numpy()\n\n        image_tensor = torch.from_numpy(image).float()\n        segmentation_tensor = torch.from_numpy(seg_data).float()\n        labels_tensor = torch.from_numpy(labels)\n        return image_tensor, segmentation_tensor, labels_tensor\n\n\ntrain_transforms = tio.Compose([tio.RandomFlip(axes=(0,), p=0.5), tio.RandomAffine(degrees=10, p=0.5), tio.ToCanonical()])\nval_transforms = tio.Compose([tio.ToCanonical()])\n\n# # Import additional MONAI modules for U-Net and DiceLoss\nfrom monai.networks.nets import UNet\nfrom monai.losses import DiceLoss\n\nclass MultiTaskUNet(nn.Module):\n    def __init__(self, num_classes=14):\n        super().__init__()\n        self.unet = UNet(\n            spatial_dims=3,\n            in_channels=1,\n            out_channels=1,  # For binary segmentation\n            channels=(16, 32, 64, 128, 256),\n            strides=(2, 2, 2, 2),\n            num_res_units=2,\n            norm=\"batch\",\n            dropout=0.1\n        )\n        # Classification head: Adjust input size to match seg_logits (1 channel after pooling)\n        self.pool = nn.AdaptiveAvgPool3d(1)\n        self.fc = nn.Linear(1, num_classes)  # Changed from 256 to 1 to match seg_logits channels\n\n    def forward(self, x):\n        # U-Net forward pass\n        seg_logits = self.unet(x)  # Shape: (B, 1, D, H, W)\n        \n        # Classification head: Pool segmentation output\n        features = self.pool(seg_logits)  # Shape: (B, 1, 1, 1, 1)\n        features = features.view(features.size(0), -1)  # Shape: (B, 1)\n        class_logits = self.fc(features)  # Shape: (B, num_classes)\n        \n        return class_logits, seg_logits\n\n\n# Training Loop\ndef train_one_epoch(model, loader, optimizer, scaler, criterion_class, criterion_seg):\n    model.train()\n    losses = []\n    for batch in tqdm(loader, desc=\"Training\"):\n        images, seg_masks, labels = batch[0].to(DEVICE), batch[1].to(DEVICE), batch[2].to(DEVICE)\n        optimizer.zero_grad()\n        with torch.cuda.amp.autocast():\n            class_outputs, seg_outputs = model(images)\n            loss_class = criterion_class(class_outputs, labels)\n            loss_seg = criterion_seg(seg_outputs, seg_masks)\n            loss = loss_class + loss_seg  # Combined loss (equal weight; adjust if needed)\n        scaler.scale(loss).backward()\n        scaler.step(optimizer)\n        scaler.update()\n        losses.append(loss.item())\n    return np.mean(losses)\n\ndef validate(model, loader, criterion_class, criterion_seg):\n    model.eval()\n    preds, truths = [], []\n    losses_class, losses_seg = [], []\n    with torch.no_grad():\n        for batch in tqdm(loader, desc=\"Validating\"):\n            images, seg_masks, labels = batch[0].to(DEVICE), batch[1].to(DEVICE), batch[2].to(DEVICE)\n            class_outputs, seg_outputs = model(images)\n            loss_class = criterion_class(class_outputs, labels)\n            loss_seg = criterion_seg(seg_outputs, seg_masks)\n            losses_class.append(loss_class.item())\n            losses_seg.append(loss_seg.item())\n            preds.append(class_outputs.cpu().numpy())\n            truths.append(labels.cpu().numpy())\n    preds = np.concatenate(preds)\n    truths = np.concatenate(truths)\n    \n    aucs = []\n    for i in range(len(LABEL_COLS)):\n        y_true = truths[:, i]\n        y_pred = preds[:, i]\n        if len(np.unique(y_true)) == 1:\n            print(f\"Warning: Only one class present in column {LABEL_COLS[i]}. AUC set to 0.5.\")\n            aucs.append(0.5)\n        else:\n            aucs.append(roc_auc_score(y_true, y_pred))\n    \n    final_score = (aucs[-1] + np.mean(aucs[:-1])) / 2\n    return np.mean(losses_class), np.mean(losses_seg), final_score, aucs\n\ndef main():\n    # Preprocess\n    preprocess_dataset()\n\n    # Data\n    metadata_file = \"/kaggle/working/series_metadata.csv\"\n    metadata = pd.read_csv(metadata_file)\n    metadata = metadata.sample(frac=1, random_state=42).reset_index(drop=True)\n    train_idx = int(0.8 * len(metadata))\n    train_data = metadata[:train_idx]\n    val_data = metadata[train_idx:]\n\n    train_dataset = AneurysmDataset(train_data, train_transforms)\n    val_dataset = AneurysmDataset(val_data, val_transforms)\n    train_loader = DataLoader(train_dataset, batch_size=BATCH_SIZE, shuffle=True, num_workers=0)\n    val_loader = DataLoader(val_dataset, batch_size=BATCH_SIZE, shuffle=False, num_workers=0)\n\n    # Model\n    model = MultiTaskUNet(num_classes=len(LABEL_COLS)).to(DEVICE)\n    criterion_class = nn.BCEWithLogitsLoss()\n    criterion_seg = DiceLoss(sigmoid=True)  # DiceLoss for segmentation\n    optimizer = torch.optim.AdamW(model.parameters(), lr=LR)\n    scaler = torch.cuda.amp.GradScaler()\n\n    # Early stopping parameters\n    patience = 5\n    counter = 0\n    best_score = 0\n    best_epoch = 0\n\n    for epoch in range(NUM_EPOCHS):\n        train_loss = train_one_epoch(model, train_loader, optimizer, scaler, criterion_class, criterion_seg)\n        val_loss_class, val_loss_seg, val_score, aucs = validate(model, val_loader, criterion_class, criterion_seg)\n\n        print(f\"Epoch {epoch+1}/{NUM_EPOCHS}, \"\n              f\"Train Loss: {train_loss:.4f}, \"\n              f\"Val Loss Class: {val_loss_class:.4f}, \"\n              f\"Val Loss Seg: {val_loss_seg:.4f}, \"\n              f\"Val AUC: {val_score:.4f}\")\n\n        # Nếu tốt hơn -> lưu lại\n        if val_score > best_score:\n            best_score = val_score\n            best_epoch = epoch\n            torch.save(model.state_dict(), \"best_model.pth\")\n            counter = 0  # reset bộ đếm\n        else:\n            counter += 1\n            print(f\"No improvement. EarlyStopping counter: {counter}/{patience}\")\n\n        # Nếu chờ đủ patience mà không cải thiện -> dừng\n        if counter >= patience:\n            print(f\"Early stopping at epoch {epoch+1}. Best epoch was {best_epoch+1} with score {best_score:.4f}\")\n            break","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-09T07:21:30.622918Z","iopub.execute_input":"2025-09-09T07:21:30.623266Z","iopub.status.idle":"2025-09-09T07:21:30.674291Z","shell.execute_reply.started":"2025-09-09T07:21:30.623244Z","shell.execute_reply":"2025-09-09T07:21:30.673528Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"if __name__ == \"__main__\":\n    main()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-09T07:21:33.442665Z","iopub.execute_input":"2025-09-09T07:21:33.443312Z","iopub.status.idle":"2025-09-09T07:25:32.117674Z","shell.execute_reply.started":"2025-09-09T07:21:33.443287Z","shell.execute_reply":"2025-09-09T07:25:32.11691Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# metadata_file = \"/kaggle/working/series_metadata.csv\"\n# metadata = pd.read_csv(metadata_file)\n# metadata.head(5)","metadata":{"trusted":true,"execution":{"execution_failed":"2025-09-09T07:08:56.991Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# print(metadata['labels'].iloc[0])","metadata":{"trusted":true,"execution":{"execution_failed":"2025-09-09T07:08:56.991Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# print(metadata['image'].iloc[0])","metadata":{"trusted":true,"execution":{"execution_failed":"2025-09-09T07:08:56.991Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Second unet\n","metadata":{}},{"cell_type":"code","source":"# class MultiTaskUNet(nn.Module):\n#     def __init__(self, num_classes=14):\n#         super().__init__()\n#         self.unet = UNet(\n#             spatial_dims=3,\n#             in_channels=1,\n#             out_channels=1,\n#             channels=(16, 32, 64, 128, 256),\n#             strides=(2, 2, 2, 2),\n#             num_res_units=2,\n#             norm=\"batch\",\n#             dropout=0.1\n#         )\n#         # Separate feature extractor for classification\n#         self.feature_extractor = nn.Sequential(\n#             nn.Conv3d(1, 16, kernel_size=3, padding=1),\n#             nn.ReLU(),\n#             nn.MaxPool3d(2),\n#             nn.Conv3d(16, 32, kernel_size=3, padding=1),\n#             nn.ReLU(),\n#             nn.MaxPool3d(2),\n#             nn.Conv3d(32, 64, kernel_size=3, padding=1),\n#             nn.ReLU(),\n#             nn.MaxPool3d(2),\n#             nn.Conv3d(64, 128, kernel_size=3, padding=1),\n#             nn.ReLU(),\n#             nn.MaxPool3d(2),\n#             nn.Conv3d(128, 256, kernel_size=3, padding=1),\n#             nn.ReLU()\n#         )\n#         self.pool = nn.AdaptiveAvgPool3d(1)\n#         self.fc = nn.Linear(256, num_classes)\n\n#     def forward(self, x):\n#         # Segmentation path\n#         seg_logits = self.unet(x)\n        \n#         # Classification path\n#         features = self.feature_extractor(x)  # Shape: (B, 256, D', H', W')\n#         features = self.pool(features)  # Shape: (B, 256, 1, 1, 1)\n#         features = features.view(features.size(0), -1)  # Shape: (B, 256)\n#         class_logits = self.fc(features)  # Shape: (B, num_classes)\n        \n#         return class_logits, seg_logits","metadata":{"trusted":true,"execution":{"execution_failed":"2025-09-09T07:08:56.991Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Submit","metadata":{}},{"cell_type":"code","source":"# import os\n# import shutil\n# from collections import defaultdict\n# import numpy as np\n# import torch\n# import torch.nn as nn\n# import pandas as pd\n# import polars as pl\n# import pydicom\n# from torch.cuda.amp import autocast\n# from tqdm import tqdm\n\n# import monai\n# import kaggle_evaluation.rsna_inference_server\n\n# # Configuration from your previous code\n# ID_COL = 'SeriesInstanceUID'\n# LABEL_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# DEVICE = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n# TARGET_SIZE = (128, 128, 64)\n# CTA_WINDOW = (-1000, 1000)\n\n# # DICOM Tag Allowlist\n# DICOM_TAG_ALLOWLIST = [\n#     'BitsAllocated', 'BitsStored', 'Columns', 'FrameOfReferenceUID', 'HighBit',\n#     'ImageOrientationPatient', 'ImagePositionPatient', 'InstanceNumber', 'Modality',\n#     'PatientID', 'PhotometricInterpretation', 'PixelRepresentation', 'PixelSpacing',\n#     'PlanarConfiguration', 'RescaleIntercept', 'RescaleSlope', 'RescaleType', 'Rows',\n#     'SOPClassUID', 'SOPInstanceUID', 'SamplesPerPixel', 'SliceThickness',\n#     'SpacingBetweenSlices', 'StudyInstanceUID', 'TransferSyntaxUID',\n# ]\n\n# # Utility Functions (simplified from your code)\n# def _safe_zoom(volume: np.ndarray, zoom_factors: tuple, order: int = 1) -> np.ndarray:\n#     \"\"\"Robust wrapper around ndimage.zoom.\"\"\"\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\n# def _resize_slice(arr: np.ndarray, out_h: int, out_w: int) -> np.ndarray:\n#     \"\"\"Resize a 2D slice to (out_h, out_w).\"\"\"\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# # Model Definition\n# class EfficientNet3D(nn.Module):\n#     def __init__(self, num_classes=14):\n#         super().__init__()\n#         self.backbone = monai.networks.nets.EfficientNetBN(\"efficientnet-b0\", spatial_dims=3, in_channels=1, pretrained=False)\n#         self.pool = nn.AdaptiveAvgPool3d(1)\n#         self.fc = nn.Linear(self.backbone._fc.in_features, num_classes)\n#         self.backbone._fc = nn.Identity()\n\n#     def forward(self, x):\n#         if x.dim() != 5:\n#             raise ValueError(f\"Expected 5D input tensor, got {x.shape}\")\n#         x = self.backbone(x)\n#         if x.dim() < 3:\n#             x = x.unsqueeze(-1).unsqueeze(-1).unsqueeze(-1)\n#         x = self.pool(x)\n#         x = x.view(x.size(0), -1)\n#         x = self.fc(x)\n#         return x\n\n# # Prediction Function\n# def predict(series_path: str) -> pl.DataFrame:\n#     \"\"\"Make a prediction for the given series.\"\"\"\n#     series_id = os.path.basename(series_path)\n    \n#     # Load all DICOM files\n#     all_filepaths = []\n#     for root, _, files in os.walk(series_path):\n#         for file in files:\n#             if file.endswith('.dcm'):\n#                 all_filepaths.append(os.path.join(root, file))\n#     all_filepaths.sort()\n\n#     # Process DICOM series into a volume\n#     vol_slices = []\n#     for filepath in tqdm(all_filepaths, desc=f\"Processing {series_id}\"):\n#         try:\n#             ds = pydicom.dcmread(filepath, force=True)\n#             arr = ds.pixel_array\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#             else:\n#                 arr = arr[np.newaxis, ...]\n#             for sl in arr:\n#                 sl = sl.astype(np.float32)\n#                 if getattr(ds, \"PhotometricInterpretation\", \"MONOCHROME2\") == \"MONOCHROME1\":\n#                     sl = sl.max() - sl\n#                 slope = float(getattr(ds, \"RescaleSlope\", 1.0))\n#                 intercept = float(getattr(ds, \"RescaleIntercept\", 0.0))\n#                 sl = sl * slope + intercept\n#                 sl = _resize_slice(sl, TARGET_SIZE[1], TARGET_SIZE[2])  # Resize 2D slices\n#                 vol_slices.append(sl)\n#         except Exception as e:\n#             print(f\"[DICOM read] {series_id}: {e}\")\n#             continue\n\n#     if not vol_slices:\n#         print(f\"No valid slices in {series_id}\")\n#         volume = np.zeros(TARGET_SIZE, dtype=np.float32)\n#     else:\n#         volume = np.stack(vol_slices, axis=0)\n#         # Resample and resize to target size\n#         z, y, x = volume.shape\n#         volume = _safe_zoom(volume, (TARGET_SIZE[0] / z, TARGET_SIZE[1] / y, TARGET_SIZE[2] / x), order=1).astype(np.float32)\n\n#     # Normalize based on modality (assuming CTA for test set)\n#     c, w = CTA_WINDOW\n#     lo, hi = c - w / 2.0, c + w / 2.0\n#     volume = np.clip((volume - lo) / (hi - lo + 1e-6), 0, 1).astype(np.float32)\n\n#     # Prepare input tensor\n#     volume = np.expand_dims(volume, axis=0)  # Add channel dimension: (1, D, H, W)\n#     volume = torch.from_numpy(volume).float().unsqueeze(0).to(DEVICE)  # Add batch dimension: (1, 1, D, H, W)\n\n#     # Load trained model\n#     model = EfficientNet3D(num_classes=len(LABEL_COLS)).to(DEVICE)\n#     model.load_state_dict(torch.load(\"best_model.pth\"))  # Load your trained model\n#     model.eval()\n\n#     # Inference\n#     with torch.no_grad(), autocast():\n#         outputs = model(volume)\n#         predictions = torch.sigmoid(outputs).cpu().numpy()[0]  # Convert logits to probabilities\n\n#     # Create prediction DataFrame\n#     predictions = pl.DataFrame(\n#         data=[[series_id] + predictions.tolist()],\n#         schema=[ID_COL, *LABEL_COLS],\n#         orient='row',\n#     )\n\n#     # Clean up temporary files\n#     shutil.rmtree('/kaggle/shared', ignore_errors=True)\n\n#     return predictions.drop(ID_COL)\n\n# # Inference Server Setup\n# inference_server = kaggle_evaluation.rsna_inference_server.RSNAInferenceServer(predict)\n\n# if os.getenv('KAGGLE_IS_COMPETITION_RERUN'):\n#     inference_server.serve()\n# else:\n#     inference_server.run_local_gateway()\n#     display(pl.read_parquet('/kaggle/working/submission.parquet'))","metadata":{"trusted":true,"execution":{"execution_failed":"2025-09-09T07:08:56.992Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## For segmentation","metadata":{}},{"cell_type":"code","source":"\n# # Preprocessing Function\n# def preprocess_dataset():\n#     os.makedirs(OUTPUT_DIR, exist_ok=True)\n#     train_df = pd.read_csv(TRAIN_CSV)\n\n#     # Keep only aneurysm cases (for now, can expand to include negatives later)\n#     train_df = train_df[train_df[\"Aneurysm Present\"] == 1].copy()\n#     print(f\"Preprocessing {len(train_df)} series with aneurysm\")\n\n#     # --- Load and parse localizers ---\n#     localizers_df = pd.read_csv(LOCALIZERS_CSV)\n\n#     def parse_coordinates(coord_str):\n#         try:\n#             d = ast.literal_eval(coord_str)\n#             return pd.Series({\n#                 \"center_x\": d.get(\"x\", np.nan),\n#                 \"center_y\": d.get(\"y\", np.nan),\n#                 \"center_z\": d.get(\"z\", np.nan) if \"z\" in d else np.nan\n#             })\n#         except Exception:\n#             return pd.Series({\"center_x\": np.nan, \"center_y\": np.nan, \"center_z\": np.nan})\n\n#     coord_df = localizers_df[\"coordinates\"].apply(parse_coordinates)\n#     localizers_df = pd.concat([localizers_df.drop(columns=[\"coordinates\"]), coord_df], axis=1)\n\n#     # Average coordinates per series\n#     localizers_df = localizers_df.groupby(\"SeriesInstanceUID\")[[\"center_x\", \"center_y\", \"center_z\"]].mean().reset_index()\n\n#     # Merge with training data\n#     train_df = train_df.merge(localizers_df, on=\"SeriesInstanceUID\", how=\"left\")\n#     train_df = train_df[:100]\n\n#     # Fill NaNs with global mean\n#     for col in [\"center_x\", \"center_y\", \"center_z\"]:\n#         train_df[col] = train_df[col].fillna(train_df[col].mean())\n\n#     # --- Prepare multiprocessing args ---\n#     processor = 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#     )\n\n#     args = [(row._asdict(), DATA_DIR, LABEL_COLS, processor) \n#             for row in train_df.itertuples(index=False)]\n\n#     series_data = []\n#     with mp.Pool(processes=min(mp.cpu_count(), 4)) as pool:\n#         results = list(tqdm(pool.imap(process_row, args),\n#                             total=len(train_df), desc=\"Preprocessing\"))\n\n#     for series_path, volume, labels in results:\n#         if volume is not None and not np.all(volume == 0):\n#             series_id = os.path.basename(series_path)\n#             np.save(os.path.join(OUTPUT_DIR, f\"{series_id}.npy\"), volume)\n            \n#             # Load and process segmentation mask\n#             seg_path = os.path.join(SEGMENTATIONS_DIR, f\"{series_id}.nii.gz\")  # Adjust extension if needed\n#             try:\n#                 seg_mask = sitk.ReadImage(seg_path)\n#                 seg_array = sitk.GetArrayFromImage(seg_mask)\n#                 # Resize segmentation mask to match volume shape\n#                 seg_array = _safe_zoom(seg_array, (TARGET_SIZE[0] / seg_array.shape[0],\n#                                                  TARGET_SIZE[1] / seg_array.shape[1],\n#                                                  TARGET_SIZE[2] / seg_array.shape[2]), order=0)\n#                 seg_array = (seg_array > 0).astype(np.float32)  # Binary mask\n#             except Exception as e:\n#                 print(f\"[Segmentation] Error loading {series_id}: {e}\")\n#                 seg_array = np.zeros(TARGET_SIZE, dtype=np.float32)\n\n#             # Crop volume using segmentation mask (simple bounding box approach)\n#             if np.any(seg_array):\n#                 coords = np.where(seg_array)\n#                 z_min, z_max = max(0, coords[0].min() - 10), min(TARGET_SIZE[0], coords[0].max() + 10)\n#                 y_min, y_max = max(0, coords[1].min() - 10), min(TARGET_SIZE[1], coords[1].max() + 10)\n#                 x_min, x_max = max(0, coords[2].min() - 10), min(TARGET_SIZE[2], coords[2].max() + 10)\n#                 volume_cropped = volume[z_min:z_max, y_min:y_max, x_min:x_max]\n#                 # Resize back to TARGET_SIZE if cropped region is smaller\n#                 if volume_cropped.shape != TARGET_SIZE:\n#                     volume_cropped = _safe_zoom(volume_cropped, (TARGET_SIZE[0] / volume_cropped.shape[0],\n#                                                                TARGET_SIZE[1] / volume_cropped.shape[1],\n#                                                                TARGET_SIZE[2] / volume_cropped.shape[2]))\n#             else:\n#                 volume_cropped = volume\n\n#             np.save(os.path.join(OUTPUT_DIR, f\"{series_id}_cropped.npy\"), volume_cropped)\n#             series_data.append({\n#                 \"image\": os.path.join(OUTPUT_DIR, f\"{series_id}_cropped.npy\"),\n#                 \"labels\": labels.tolist()\n#             })\n\n#     pd.DataFrame(series_data).to_csv(os.path.join(OUTPUT_DIR, \"series_metadata.csv\"), index=False)\n#     print(\"Preprocessing done!\")\n\n# # Dataset\n# class AneurysmDataset(Dataset):\n#     def __init__(self, data_file: str, transform=None):\n#         self.data = pd.read_csv(data_file)\n#         self.transform = transform\n\n#     def __len__(self):\n#         return len(self.data)\n\n#     def __getitem__(self, idx):\n#         row = self.data.iloc[idx]\n#         image = np.load(row[\"image\"]).astype(np.float32)\n#         # Ensure image is 3D (D, H, W) and add channel dimension\n#         if image.ndim == 3:\n#             image = np.expand_dims(image, axis=0)  # Shape: (1, D, H, W)\n#         elif image.ndim == 4 and image.shape[0] == 1:\n#             image = image.squeeze(0)  # Ensure no extra batch dimension from transform\n        \n#         # Parse labels string using ast.literal_eval\n#         labels = np.array(ast.literal_eval(row[\"labels\"]), dtype=np.float32) if isinstance(row[\"labels\"], str) else np.array(row[\"labels\"], dtype=np.float32)\n        \n#         # Light augmentation with torchio\n#         if self.transform:\n#             subject = tio.Subject(image=tio.ScalarImage(tensor=image))\n#             transformed = self.transform(subject)\n#             image = transformed[\"image\"].numpy()\n\n#         # Ensure 5D tensor for model (batch is added by DataLoader)\n#         image_tensor = torch.from_numpy(image).float()  # Shape: (1, D, H, W)\n#         labels_tensor = torch.from_numpy(labels)\n#         return image_tensor, labels_tensor","metadata":{"trusted":true,"execution":{"execution_failed":"2025-09-09T07:08:56.992Z"}},"outputs":[],"execution_count":null}]}