{"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":[{"sourceId":99552,"databundleVersionId":13747926,"sourceType":"competition"},{"sourceId":12780021,"sourceType":"datasetVersion","datasetId":8079690},{"sourceId":12885547,"sourceType":"datasetVersion","datasetId":8030314},{"sourceId":13050805,"sourceType":"datasetVersion","datasetId":8009945},{"sourceId":227202990,"sourceType":"kernelVersion"},{"sourceId":549965,"sourceType":"modelInstanceVersion","isSourceIdPinned":true,"modelInstanceId":419249,"modelId":436908}],"dockerImageVersionId":31090,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"from pathlib import Path\nfrom typing import List, Tuple, Dict, Optional\nimport cv2\n\nimport sys\nimport matplotlib.pyplot as plt \nsys.path.append('/kaggle/input/rsna-iad-vesselfm-codebase')\nimport torch\nimport torch.nn.functional as F\nimport numpy as np\nfrom tqdm import tqdm\nfrom monai.inferers import SlidingWindowInfererAdapt\nfrom skimage.morphology import remove_small_objects\nfrom skimage.exposure import equalize_hist\nfrom utils.data import generate_transforms\nfrom utils.io import determine_reader_writer\nimport os\nfrom monai.transforms import LoadImaged, Spacingd, LoadImage\nfrom monai.networks.nets import DynUNet\nimport SimpleITK as sitk\nimport yaml\nimport torch.nn as nn\nfrom scipy.ndimage import label\nfrom tqdm import  tqdm\nimport pydicom\nfrom concurrent.futures import ThreadPoolExecutor\nfrom collections import Counter\nimport pandas as pd\nfrom matplotlib.widgets import Slider\nimport ipywidgets as widgets\nfrom IPython.display import HTML\nfrom matplotlib.animation import FuncAnimation\nfrom monai.transforms import (\n    Compose,\n    EnsureChannelFirstd,\n    EnsureTyped,\n    Resized,\n    RandCropByLabelClassesd,\n    SpatialPadd,\n    ConcatItemsd,\n    ToTensord,\n    Lambdad,\n)\n\nfrom monai.transforms import MapTransform\nfrom skimage import morphology, filters","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:50:41.316051Z","iopub.execute_input":"2025-09-16T12:50:41.316806Z","iopub.status.idle":"2025-09-16T12:51:10.874449Z","shell.execute_reply.started":"2025-09-16T12:50:41.316775Z","shell.execute_reply":"2025-09-16T12:51:10.873857Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"yaml_path = '/kaggle/input/rsna-iad-vesselfm-codebase/configs/inference.yaml'\nMAX_WORKERS = 4\nwith open(yaml_path, 'r') as f:\n    config = yaml.safe_load(f)\n\n\nclass CFG:\n    ckpt_path = '/kaggle/input/rsna-iad-vesselfm-13-classes/pytorch/default/2/vesselfm_13_classes_dynunet-val_dice0.5047.ckpt'\n    model_structure = {\n        'in_channels': 1,\n        'out_channels': 14,\n        'spatial_dims': 3,\n        'strides': [[1, 1, 1], [2, 2, 2], [2, 2, 2], [2, 2, 2], [2, 2, 2]] , # 5 levels\n        'kernel_size': [[3, 3, 3], [3, 3, 3], [3, 3, 3], [3, 3, 3], [3, 3, 3]],\n        'upsample_kernel_size': [[2, 2, 2], [2, 2, 2], [2, 2, 2], [2, 2, 2]],\n        'filters': [32, 64, 128, 256, 320],\n        'res_block': True}\n    device = 'cuda:0'\n    thrd = 0.1\n\n    #sliding window\n    batch_size= 1\n    patch_size= [128, 128, 128]\n    overlap= 0.5\n    mode= \"constant\"\n    sigma_scale= 0.125\n    padding_mode= \"constant\"\n\n    #volume transform\n    transforms_config = config['transforms_config']\n    transforms_config.insert(1, {'Resize': {\n        'spatial_size': (128, 384, 384),\n        'mode': \"trilinear\"\n    }},)\n    \n    tta = config['tta']\n    post = config['post']\n    merging = config['merging']","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:51:10.875614Z","iopub.execute_input":"2025-09-16T12:51:10.876255Z","iopub.status.idle":"2025-09-16T12:51:10.89423Z","shell.execute_reply.started":"2025-09-16T12:51:10.876233Z","shell.execute_reply":"2025-09-16T12:51:10.893504Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def load_dicom_series(series_dir: Path) -> Tuple[np.ndarray, np.ndarray]:\n    \"\"\"\n    Load a DICOM series into a 3D HU volume with affine matrix.\n\n    Args:\n        series_dir (Path): directory containing DICOM files\n\n    Returns:\n        volume (np.ndarray): 3D array (Z, Y, X) in HU\n        affine (np.ndarray): 4x4 voxel-to-world affine\n    \"\"\"\n    # Collect DICOMs\n    dcm_paths = [Path(series_dir) / f for f in os.listdir(series_dir) if f.lower().endswith(\".dcm\")]\n    if not dcm_paths:\n        raise FileNotFoundError(f\"No DICOM files found in {series_dir}\")\n    \n    slices = [pydicom.dcmread(str(p), force=True) for p in dcm_paths]\n\n    # Orientation & sorting\n    orientation = np.array(slices[0].ImageOrientationPatient).reshape(2, 3)\n    row_cos, col_cos = orientation\n    normal = np.cross(row_cos, col_cos)\n    slices.sort(key=lambda ds: np.dot(ds.ImagePositionPatient, normal))\n\n    # HU scaling\n    slope = float(getattr(slices[0], \"RescaleSlope\", 1.0))\n    intercept = float(getattr(slices[0], \"RescaleIntercept\", 0.0))\n    volume = np.stack([ds.pixel_array for ds in slices]).astype(np.float32)\n    volume = volume * slope + intercept\n\n    # Spacing\n    Δx, Δy = [float(x) for x in slices[0].PixelSpacing]\n    positions = [np.array(ds.ImagePositionPatient) for ds in slices]\n    slice_positions = [np.dot(p, normal) for p in positions]\n    Δz = float(np.mean(np.diff(slice_positions)))\n\n    # Affine\n    slice_cos = normal\n    origin = np.array(slices[0].ImagePositionPatient)\n    affine = np.eye(4)\n    affine[0:3, 0] = row_cos * Δx\n    affine[0:3, 1] = col_cos * Δy\n    affine[0:3, 2] = slice_cos * Δz\n    affine[0:3, 3] = origin\n\n    return volume\n\n\nclass ModalityIntensityScalingd(MapTransform):\n    \"\"\"\n    Modality-agnostic normalization:\n      - Robust z-score within 2–98 percentiles\n      - Then rescale to [0, 1]\n    \"\"\"\n    def __init__(self, keys=(\"Image\",)):\n        super().__init__(keys)\n\n    def __call__(self, data):\n        d = dict(data)\n        img = d[self.keys[0]]\n\n        # robust stats\n        p2, p98 = np.percentile(img, (2, 98))\n        mask = (img >= p2) & (img <= p98)\n        mean = np.mean(img[mask])\n        std = np.std(img[mask]) + 1e-6\n\n        img = (img - mean) / std\n        img = np.clip((img - img.min()) / (img.max() - img.min() + 1e-6), 0, 1)\n\n        d[self.keys[0]] = img.astype(np.float32)\n        return d\n\ndef _generate_transforms(vol_size, input_size, mode):\n    if mode == \"train\":\n        return Compose([\n            EnsureChannelFirstd(keys=[\"Image\", \"Mask\"], channel_dim=\"no_channel\"),\n            EnsureTyped(keys=[\"Image\", \"Mask\"]),\n            Resized(keys=[\"Image\", \"Mask\"], spatial_size=vol_size, mode=[\"trilinear\", \"nearest\"]),\n            RandCropByLabelClassesd(\n                keys=[\"Image\", \"Mask\"],\n                label_key=\"Mask\",\n                spatial_size=input_size,\n                num_classes=13,\n                ratios=[1] * 13,\n                num_samples=4,\n                image_key=\"Image\",\n                allow_smaller=True,\n            ),\n            ModalityIntensityScalingd(keys=[\"Image\"]),\n            SpatialPadd(keys=[\"Image\", \"Mask\"], spatial_size=input_size, mode=\"constant\", method=\"symmetric\"),\n            ConcatItemsd(keys=[\"Image\", \"Mask\"], name=[\"Image\", \"Mask\"], dim=0),\n            ToTensord(keys=[\"Image\", \"Mask\"]),\n        ])\n    elif mode == 'val':\n        return Compose([\n            EnsureChannelFirstd(keys=[\"Image\", \"Mask\"], channel_dim=\"no_channel\"),\n            EnsureTyped(keys=[\"Image\", \"Mask\"]),\n            Resized(keys=[\"Image\", \"Mask\"], spatial_size=vol_size, mode=[\"trilinear\", \"nearest\"]),\n            ModalityIntensityScalingd(keys=[\"Image\"]),\n            ToTensord(keys=[\"Image\", \"Mask\"]),\n        ])\n    else:\n        return Compose([\n            EnsureChannelFirstd(keys=[\"Image\"], channel_dim=\"no_channel\"),\n            EnsureTyped(keys=[\"Image\"]),\n            Resized(keys=[\"Image\"], spatial_size=vol_size, mode=[\"trilinear\"]),\n            ModalityIntensityScalingd(keys=[\"Image\"]),\n            ToTensord(keys=[\"Image\"]),\n        ])\n\n\n\ndef cca_remove(mask_preds, min_size = 1000, sigma = 0.5):\n    '''\n    mask_preds: (D, H, W)\n    '''\n    #to 2d mip\n    mip_2d = mask_preds.max(axis=0)\n    \n    # Remove small bright spots (choose an area threshold)\n    smoothed_mip = filters.gaussian(mip_2d, sigma=sigma)\n    cleaned_mip = morphology.remove_small_objects(smoothed_mip > 0, min_size=min_size)  \n    \n    # Morphological opening for noise removal\n    selem = morphology.disk(2)  \n    cleaned_mip = morphology.opening(cleaned_mip, selem)\n    \n    # Convert boolean mask back to vessel intensity (optional thresholding)\n    cleaned_mip = cleaned_mip.astype(np.uint8) * mip_2d.max()\n    \n    # Median filter smoothing\n    cleaned_mip_mask = filters.median(cleaned_mip, morphology.disk(1))\n\n    masked_vol = (cleaned_mip_mask[None]!=0) * mask_preds\n\n    new_mip_mask = masked_vol.max(axis=0)\n    \n    return masked_vol, new_mip_mask","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:51:10.895149Z","iopub.execute_input":"2025-09-16T12:51:10.895489Z","iopub.status.idle":"2025-09-16T12:51:10.922303Z","shell.execute_reply.started":"2025-09-16T12:51:10.895466Z","shell.execute_reply":"2025-09-16T12:51:10.921673Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def load_model(cfg):\n    ckpt = torch.load(cfg.ckpt_path, map_location=cfg.device, weights_only=False)['state_dict']\n    ckpt = {k.replace(\"model.\", \"\"): v for k, v in ckpt.items()}\n    model = DynUNet(**CFG.model_structure)\n    model.load_state_dict(ckpt)\n    model.eval()\n    return model.to(cfg.device)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:51:10.924062Z","iopub.execute_input":"2025-09-16T12:51:10.924268Z","iopub.status.idle":"2025-09-16T12:51:10.939652Z","shell.execute_reply.started":"2025-09-16T12:51:10.924251Z","shell.execute_reply":"2025-09-16T12:51:10.938765Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"sliding_window_inferer = SlidingWindowInfererAdapt(\n            roi_size=CFG.patch_size, sw_batch_size=1, overlap=0.5,\n        )","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:51:10.94059Z","iopub.execute_input":"2025-09-16T12:51:10.940843Z","iopub.status.idle":"2025-09-16T12:51:10.954041Z","shell.execute_reply.started":"2025-09-16T12:51:10.940826Z","shell.execute_reply":"2025-09-16T12:51:10.953312Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"transforms = generate_transforms(CFG.transforms_config)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:51:10.954805Z","iopub.execute_input":"2025-09-16T12:51:10.955058Z","iopub.status.idle":"2025-09-16T12:51:10.970117Z","shell.execute_reply.started":"2025-09-16T12:51:10.955036Z","shell.execute_reply":"2025-09-16T12:51:10.969395Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"model = load_model(CFG)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:51:10.970737Z","iopub.execute_input":"2025-09-16T12:51:10.970944Z","iopub.status.idle":"2025-09-16T12:51:14.065861Z","shell.execute_reply.started":"2025-09-16T12:51:10.970927Z","shell.execute_reply":"2025-09-16T12:51:14.065083Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"df = pd.read_csv('/kaggle/input/rsna-intracranial-aneurysm-detection/train.csv')\ndf_loc = pd.read_csv('/kaggle/input/rsna-intracranial-aneurysm-detection/train_localizers.csv')\nseg_path = Path('/kaggle/input/rsna-intracranial-aneurysm-detection/segmentations')\nseries_path = Path('/kaggle/input/rsna-intracranial-aneurysm-detection/series')\nseg_uids = [name.split('.nii')[0] for name in os.listdir(seg_path) if 'cowseg' not in name]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:51:14.06677Z","iopub.execute_input":"2025-09-16T12:51:14.067297Z","iopub.status.idle":"2025-09-16T12:51:14.220157Z","shell.execute_reply.started":"2025-09-16T12:51:14.067273Z","shell.execute_reply":"2025-09-16T12:51:14.219574Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#uid = '1.2.826.0.1.3680043.8.498.10004044428023505108375152878107656647'\n#uid = '1.2.826.0.1.3680043.8.498.10004044428023505108375152878107656647'\nuid = '1.2.826.0.1.3680043.8.498.10005158603912009425635473100344077317'\nprint(uid in seg_uids)\npath = series_path/uid\nmask_path = seg_path/f'{uid}_cowseg.nii'\nvol_path = seg_path/f'{uid}.nii'\ndf[df.SeriesInstanceUID\t== uid]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:51:14.220902Z","iopub.execute_input":"2025-09-16T12:51:14.221115Z","iopub.status.idle":"2025-09-16T12:51:14.249699Z","shell.execute_reply.started":"2025-09-16T12:51:14.221098Z","shell.execute_reply":"2025-09-16T12:51:14.249058Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"vol_ = load_dicom_series(path)\nvol = np.flip(vol_, axis=1).copy()\nvol = np.flip(vol, axis=2).copy()\n\ntransform_test = _generate_transforms((128, 384, 384), (128, 128, 128), 'test')\n_input = {'Image': vol}\n\ntest_data = transform_test(_input)\ntest_vol = test_data['Image'].to(CFG.device)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:51:14.25227Z","iopub.execute_input":"2025-09-16T12:51:14.252689Z","iopub.status.idle":"2025-09-16T12:51:17.2081Z","shell.execute_reply.started":"2025-09-16T12:51:14.252671Z","shell.execute_reply":"2025-09-16T12:51:17.207322Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"with torch.no_grad():\n    pred_logits = sliding_window_inferer(test_vol[None], model)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:51:17.208983Z","iopub.execute_input":"2025-09-16T12:51:17.209219Z","iopub.status.idle":"2025-09-16T12:51:21.757326Z","shell.execute_reply.started":"2025-09-16T12:51:17.2092Z","shell.execute_reply":"2025-09-16T12:51:21.756593Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"test_mip, _ = test_vol[0].max(dim=0)\ntest_mip = test_mip.cpu()\n\npred_probs = pred_logits.softmax(dim=1)\npred_mask = pred_logits.argmax(dim=1)\n\nmask_mip, _ = pred_mask[0].max(dim=0)\nmask_mip = mask_mip.cpu()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:51:21.758511Z","iopub.execute_input":"2025-09-16T12:51:21.758814Z","iopub.status.idle":"2025-09-16T12:51:21.843855Z","shell.execute_reply.started":"2025-09-16T12:51:21.758789Z","shell.execute_reply":"2025-09-16T12:51:21.843294Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"new_mask, new_mask_mip = cca_remove(pred_mask[0].cpu().numpy())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:51:21.844526Z","iopub.execute_input":"2025-09-16T12:51:21.844746Z","iopub.status.idle":"2025-09-16T12:51:22.148876Z","shell.execute_reply.started":"2025-09-16T12:51:21.84473Z","shell.execute_reply":"2025-09-16T12:51:22.148078Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"fig, ax = plt.subplots(ncols=3, figsize = (16, 16))\nax[0].imshow(test_mip)\nax[1].imshow(mask_mip)\nax[1].set_title('before process')\nax[2].imshow(new_mask_mip)\nax[2].set_title('after process')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:51:22.149824Z","iopub.execute_input":"2025-09-16T12:51:22.150828Z","execution_failed":"2025-09-16T12:53:01.971Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Filtering Invalid YOLO Location Prediction ","metadata":{}},{"cell_type":"code","source":"!tar xfvz /kaggle/input/ultralytics-offlineinstall-yolo12-weights/archive.tar.gz\n!pip install --no-index --find-links=./packages ultralytics\n!rm -rf ./packages","metadata":{"trusted":true,"_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2025-09-16T12:51:22.758078Z","iopub.execute_input":"2025-09-16T12:51:22.758298Z","iopub.status.idle":"2025-09-16T12:52:09.868219Z","shell.execute_reply.started":"2025-09-16T12:51:22.758281Z","shell.execute_reply":"2025-09-16T12:52:09.866955Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ML/DL\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nfrom torch.cuda.amp import autocast\nimport timm\n\n# Transformations\nimport albumentations as A\nfrom albumentations.pytorch import ToTensorV2\n\n# YOLO\nfrom ultralytics import YOLO\nimport matplotlib.colors as mcolors\nimport matplotlib.patches as mpatches\nimport matplotlib.cm as cm\n\n# Set device\ndevice = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n\n# Optimization settings\ntorch.set_float32_matmul_precision('medium')\ntorch.backends.cudnn.benchmark = True\ntorch.backends.cudnn.deterministic = False\ntorch.backends.cuda.matmul.allow_tf32 = True\ntorch.backends.cudnn.allow_tf32 = True","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:52:09.869405Z","iopub.execute_input":"2025-09-16T12:52:09.870065Z","iopub.status.idle":"2025-09-16T12:52:14.278341Z","shell.execute_reply.started":"2025-09-16T12:52:09.870038Z","shell.execute_reply":"2025-09-16T12:52:14.27782Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def load_yolo_models():\n    \"\"\"Load all YOLO models\"\"\"\n    models = []\n    for config in YOLO_MODEL_CONFIGS:\n        model = YOLO(config[\"path\"])\n        model.to(device)\n        \n        model_dict = {\n            \"model\": model,\n            \"weight\": config[\"weight\"],\n            \"name\": config[\"name\"],\n            \"fold\": config[\"fold\"]\n        }\n        models.append(model_dict)\n    return models","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:52:14.27916Z","iopub.execute_input":"2025-09-16T12:52:14.279824Z","iopub.status.idle":"2025-09-16T12:52:14.284058Z","shell.execute_reply.started":"2025-09-16T12:52:14.279796Z","shell.execute_reply":"2025-09-16T12:52:14.283228Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# YOLO label mappings\nYOLO_LABELS_TO_IDX = {\n    'Anterior Communicating Artery': 0,\n    'Basilar Tip': 1,\n    'Left Anterior Cerebral Artery': 2,\n    'Left Infraclinoid Internal Carotid Artery': 3,\n    'Left Middle Cerebral Artery': 4,\n    'Left Posterior Communicating Artery': 5,\n    'Left Supraclinoid Internal Carotid Artery': 6,\n    'Other Posterior Circulation': 7,\n    'Right Anterior Cerebral Artery': 8,\n    'Right Infraclinoid Internal Carotid Artery': 9,\n    'Right Middle Cerebral Artery': 10,\n    'Right Posterior Communicating Artery': 11,\n    'Right Supraclinoid Internal Carotid Artery': 12\n}\nYOLO_LABELS = sorted(list(YOLO_LABELS_TO_IDX.keys()))\n\n\n# ====================================================\nIMG_SIZE = 512\nBATCH_SIZE = int(os.getenv(\"YOLO_BATCH_SIZE\", \"32\"))\nMAX_WORKERS = 4\n\nYOLO_MODEL_CONFIGS = [\n    {\n        \"path\": \"/kaggle/input/rsna-sergio-models/cv_y11m_with_mix_up_mosaic_fold0/weights/best.pt\",\n        \"fold\": \"0\",\n        \"weight\": 1.0,\n        \"name\": \"YOLOv11n_fold0\"\n    },\n    {\n        \"path\": \"/kaggle/input/rsna-sergio-models/cv_y11m_with_mix_up_mosaic_fold1/weights/best.pt\",\n        \"fold\": \"1\",\n        \"weight\": 1.0,\n        \"name\": \"YOLOv11n_fold1\"\n    },  \n    #{\n    #    \"path\": \"/kaggle/input/rsna-sergio-models/cv_y11m_with_mix_up_mosaic_fold2/weights/best.pt\",\n    #    \"fold\": \"2\",\n    #    \"weight\": 1.0,\n    #    \"name\": \"YOLOv11n_fold2\"\n    #}\n]\n\nYOLO_MODELS = load_yolo_models()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:52:14.284921Z","iopub.execute_input":"2025-09-16T12:52:14.285188Z","iopub.status.idle":"2025-09-16T12:52:15.997662Z","shell.execute_reply.started":"2025-09-16T12:52:14.285157Z","shell.execute_reply":"2025-09-16T12:52:15.996871Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def read_dicom_frames_hu(path: Path) -> List[Tuple[float, np.ndarray]]:\n    \"\"\"Read DICOM file and return list of (slice_position, HU frame)\"\"\"\n    ds = pydicom.dcmread(str(path), force=True)\n    pix = ds.pixel_array\n    slope = float(getattr(ds, 'RescaleSlope', 1.0))\n    intercept = float(getattr(ds, 'RescaleIntercept', 0.0))\n\n    # Compute slice location using orientation + position\n    try:\n        orientation = np.array(ds.ImageOrientationPatient).reshape(2, 3)\n        row_cos, col_cos = orientation\n        normal = np.cross(row_cos, col_cos)  # slice normal vector\n        position = np.array(ds.ImagePositionPatient)\n        slice_loc = float(np.dot(position, normal))  # projection along normal\n    except Exception:\n        # Fallback: SliceLocation / InstanceNumber\n        slice_loc = float(getattr(ds, \"SliceLocation\", getattr(ds, \"InstanceNumber\", 0.0)))\n\n    frames: List[Tuple[float, np.ndarray]] = []\n\n    if pix.ndim == 2:\n        img = pix.astype(np.float32)\n        frames.append((slice_loc, img * slope + intercept))\n    elif pix.ndim == 3:\n        # RGB or multi-frame\n        if pix.shape[-1] == 3 and pix.shape[0] != 3:\n            try:\n                gray = cv2.cvtColor(pix.astype(np.uint8), cv2.COLOR_BGR2GRAY).astype(np.float32)\n            except Exception:\n                gray = pix[..., 0].astype(np.float32)\n            frames.append((slice_loc, gray * slope + intercept))\n        else:\n            for i in range(pix.shape[0]):\n                frm = pix[i].astype(np.float32)\n                # tiny offset ensures consistent ordering for multi-frame\n                frames.append((slice_loc + i * 1e-3, frm * slope + intercept))\n    return frames\n\n\ndef min_max_normalize(img: np.ndarray) -> np.ndarray:\n    \"\"\"Min-max normalization to 0-255 with optional flipping\"\"\"\n    mn, mx = float(img.min()), float(img.max())\n    if mx - mn < 1e-6:\n        norm = np.zeros_like(img, dtype=np.uint8)\n    else:\n        norm = (img - mn) / (mx - mn)\n        norm = (norm * 255.0).clip(0, 255).astype(np.uint8)\n    return norm\n\n\ndef process_dicom_file(dcm_path: Path) -> List[Tuple[float, np.ndarray]]:\n    \"\"\"Process single DICOM file -> list of (slice_loc, image) tuples\"\"\"\n    try:\n        frames = read_dicom_frames_hu(dcm_path)\n        processed_slices = []\n        for loc, f in frames:\n            img_u8 = min_max_normalize(f)\n            if img_u8.ndim == 2:\n                img_u8 = cv2.cvtColor(img_u8, cv2.COLOR_GRAY2BGR)\n            processed_slices.append((loc, img_u8))\n        return processed_slices\n    except Exception as e:\n        print(f\"Failed processing {dcm_path.name}: {e}\")\n        return []\n\n\ndef collect_series_slices(series_dir: Path) -> List[Path]:\n    \"\"\"Collect all DICOM files in a series directory (recursively).\"\"\"\n    dcm_paths: List[Path] = []\n    try:\n        for root, _, files in os.walk(series_dir):\n            for f in files:\n                if f.lower().endswith('.dcm'):\n                    dcm_paths.append(Path(root) / f)\n    except Exception as e:\n        print(f\"Failed to walk series dir {series_dir}: {e}\")\n    return dcm_paths\n\n\ndef slice_sort_key(path: Path) -> float:\n    \"\"\"Compute a robust slice sort key (orientation + position) for a single DICOM file\"\"\"\n    try:\n        ds = pydicom.dcmread(str(path), stop_before_pixels=True, force=True)\n        orientation = np.array(ds.ImageOrientationPatient).reshape(2, 3)\n        row_cos, col_cos = orientation\n        normal = np.cross(row_cos, col_cos)\n        position = np.array(ds.ImagePositionPatient)\n        return float(np.dot(position, normal))\n    except Exception:\n        # fallback\n        try:\n            return float(getattr(ds, \"SliceLocation\", getattr(ds, \"InstanceNumber\", 0.0)))\n        except:\n            return 0.0\n\ndef process_dicom_for_yolo(series_path):\n    series_path = Path(series_path)\n    dicom_files = collect_series_slices(series_path)\n    \n    # Sort DICOM files by orientation+position before processing\n    dicom_files.sort(key=slice_sort_key)\n    \n    with ThreadPoolExecutor(max_workers=MAX_WORKERS) as executor:\n        results = list(executor.map(process_dicom_file, dicom_files))\n    \n    # Flatten into (loc, img)\n    all_slices_with_loc = [item for sublist in results for item in sublist]\n    \n    # Already sorted by dicom_files order, but double-check (safe)\n    all_slices_with_loc.sort(key=lambda x: x[0])\n    \n    # Extract just the images\n    all_slices = [img for _, img in all_slices_with_loc]\n    \n    # Now dicom_files matches the sorted slices\n    dcm_list = [f.stem for f in dicom_files]\n    return all_slices","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:52:15.998547Z","iopub.execute_input":"2025-09-16T12:52:15.999385Z","iopub.status.idle":"2025-09-16T12:52:16.013784Z","shell.execute_reply.started":"2025-09-16T12:52:15.999357Z","shell.execute_reply":"2025-09-16T12:52:16.012992Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"@torch.no_grad()\ndef predict_yolo_ensemble(slices):\n    if not slices:\n        return 0.1, np.ones(len(YOLO_LABELS)) * 0.1\n    ensemble_cls_preds = []\n    ensemble_loc_preds = []\n    total_weight = 0.0\n    vol_size = (len(slices), slices[0].shape[0], slices[0].shape[1])\n    location_preds = {f'MODEL{i}': [] for i in range(len(YOLO_MODELS))}\n    \n    for model_idx, model_dict in enumerate(YOLO_MODELS):\n        model = model_dict[\"model\"]\n        weight = model_dict[\"weight\"]\n        all_detections = []\n        try:\n            max_conf_all = 0.0\n            per_class_max = np.zeros(len(YOLO_LABELS), dtype=np.float32)\n            \n            # Process in batches\n            for i in range(0, len(slices), BATCH_SIZE):\n                batch_slices = slices[i:i+BATCH_SIZE]\n                z_idxes = [i + batch_idx  for batch_idx in range(len(batch_slices))]\n                results = model.predict(\n                    batch_slices, \n                    verbose=False, \n                    batch=len(batch_slices), \n                    device=\"cuda:0\", \n                    conf=0.01\n                )\n                \n                for z_idx, r in enumerate(results):\n                    if r is None or r.boxes is None or r.boxes.conf is None or len(r.boxes) == 0:\n                        continue\n                    try:\n                        confs = r.boxes.conf\n                        clses = r.boxes.cls\n                        for j in range(len(confs)):\n                            c = float(confs[j].item())\n                            k = int(clses[j].item())\n                            if c > max_conf_all:\n                                max_conf_all = c\n                            if 0 <= k < len(YOLO_LABELS) and c > per_class_max[k]:\n                                per_class_max[k] = c\n                            x1, y1, x2, y2 = r.boxes.xyxy[j].cpu().numpy()\n                            x_center = (x1 + x2) / 2\n                            y_center = (y1 + y2) / 2\n                            point = np.array([round(z_idxes[z_idx]), round(y_center), round(x_center)])\n                            location_preds[f'MODEL{model_idx}'].append([*point, float(c), k, model_idx]) #z, y, x, prob, class, model index\n                            \n                    except Exception as e:\n                        print(e)\n                        try:\n                            batch_max = float(r.boxes.conf.max().item())\n                            if batch_max > max_conf_all:\n                                max_conf_all = batch_max\n                        except Exception:\n                            pass\n            ensemble_cls_preds.append(max_conf_all * weight)\n            ensemble_loc_preds.append(per_class_max * weight)\n            total_weight += weight\n            \n        except Exception as e:\n            print(f\"Error in model {model_dict['name']}: {e}\")\n            ensemble_cls_preds.append(0.1 * weight)\n            ensemble_loc_preds.append(np.ones(len(YOLO_LABELS)) * 0.1 * weight)\n            total_weight += weight\n\n    location_preds = {k: np.array(v) for k, v in location_preds.items()}\n    \n    if total_weight > 0:\n        final_cls_pred = sum(ensemble_cls_preds) / len(ensemble_cls_preds)\n        final_loc_preds = sum(ensemble_loc_preds) / total_weight\n    else:\n        final_cls_pred = 0.1\n        final_loc_preds = np.ones(len(YOLO_LABELS)) * 0.1\n    return final_cls_pred, final_loc_preds, location_preds","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:52:16.014507Z","iopub.execute_input":"2025-09-16T12:52:16.014811Z","iopub.status.idle":"2025-09-16T12:52:16.034483Z","shell.execute_reply.started":"2025-09-16T12:52:16.01479Z","shell.execute_reply":"2025-09-16T12:52:16.033801Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"all_slices = process_dicom_for_yolo(path)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:52:16.035315Z","iopub.execute_input":"2025-09-16T12:52:16.036101Z","iopub.status.idle":"2025-09-16T12:52:17.192413Z","shell.execute_reply.started":"2025-09-16T12:52:16.036078Z","shell.execute_reply":"2025-09-16T12:52:17.191867Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"final_cls_pred, final_loc_preds, location_preds = predict_yolo_ensemble(all_slices)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:52:17.193126Z","iopub.execute_input":"2025-09-16T12:52:17.193327Z","iopub.status.idle":"2025-09-16T12:52:26.89885Z","shell.execute_reply.started":"2025-09-16T12:52:17.19331Z","shell.execute_reply":"2025-09-16T12:52:26.898255Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"SEG_LABELS_TO_IDX = {\n    \"BG\": 0,\n    \"Other Posterior Circulation\": 1,\n    \"Basilar Tip\": 2,\n    \"Right Posterior Communicating Artery\": 3,\n    \"Left Posterior Communicating Artery\": 4,\n    \"Right Infraclinoid Internal Carotid Artery\": 5,\n    \"Left Infraclinoid Internal Carotid Artery\": 6,\n    \"Right Supraclinoid Internal Carotid Artery\": 7,\n    \"Left Supraclinoid Internal Carotid Artery\": 8,\n    \"Right Middle Cerebral Artery\": 9,\n    \"Left Middle Cerebral Artery\": 10,\n    \"Right Anterior Cerebral Artery\": 11,\n    \"Left Anterior Cerebral Artery\": 12,\n    \"Anterior Communicating Artery\": 13,\n}\n\n# Invert SEG map\nSEG_IDX_TO_LABELS = {v: k for k, v in SEG_LABELS_TO_IDX.items()}\n\n# Build YOLO→SEG mapping\nYOLO_TO_SEG = {\n    yolo_idx: SEG_LABELS_TO_IDX[name]\n    for name, yolo_idx in YOLO_LABELS_TO_IDX.items()\n    if name in SEG_LABELS_TO_IDX\n}\n\nmax_yolo_idx = max(YOLO_TO_SEG.keys())\nyolo2seg_lookup = np.zeros(max_yolo_idx + 1, dtype=int)\nfor yolo_idx, seg_idx in YOLO_TO_SEG.items():\n    yolo2seg_lookup[yolo_idx] = seg_idx","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:52:26.899639Z","iopub.execute_input":"2025-09-16T12:52:26.899891Z","iopub.status.idle":"2025-09-16T12:52:26.905442Z","shell.execute_reply.started":"2025-09-16T12:52:26.899867Z","shell.execute_reply":"2025-09-16T12:52:26.904587Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def resize_mask_3d(volume, target_size):\n    \"\"\"\n    Resize a 3D volume [D, H, W] to [d, h, w] using nearest neighbor interpolation (SimpleITK).\n    \"\"\"\n    d, h, w = target_size\n    image = sitk.GetImageFromArray(volume)\n    old_size = image.GetSize()\n    old_spacing = image.GetSpacing()\n    \n    new_size = [w, h, d]  # SITK uses (x, y, z)\n    new_spacing = [\n        old_spacing[i] * (old_size[i] / new_size[i]) for i in range(3)\n    ]\n    \n    resampler = sitk.ResampleImageFilter()\n    resampler.SetSize(new_size)\n    resampler.SetOutputSpacing(new_spacing)\n    resampler.SetInterpolator(sitk.sitkNearestNeighbor)\n    \n    resized_img = resampler.Execute(image)\n    reszied_mask = sitk.GetArrayFromImage(resized_img)\n    return np.flip(np.flip(reszied_mask, axis=2), axis=1)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:52:26.90606Z","iopub.execute_input":"2025-09-16T12:52:26.906252Z","iopub.status.idle":"2025-09-16T12:52:26.92134Z","shell.execute_reply.started":"2025-09-16T12:52:26.906237Z","shell.execute_reply":"2025-09-16T12:52:26.92072Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"all_slices_vol = np.array(all_slices)[..., 0]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:52:26.922038Z","iopub.execute_input":"2025-09-16T12:52:26.92223Z","iopub.status.idle":"2025-09-16T12:52:26.946092Z","shell.execute_reply.started":"2025-09-16T12:52:26.922215Z","shell.execute_reply":"2025-09-16T12:52:26.945493Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"resized_mask = resize_mask_3d(new_mask, all_slices_vol.shape)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:52:26.946873Z","iopub.execute_input":"2025-09-16T12:52:26.947143Z","iopub.status.idle":"2025-09-16T12:52:27.32183Z","shell.execute_reply.started":"2025-09-16T12:52:26.947119Z","shell.execute_reply":"2025-09-16T12:52:27.321161Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"loc_preds_2d_1 = location_preds['MODEL0'][:, [1, 2]] #y, x\nloc_preds_2d_2 = location_preds['MODEL1'][:, [1, 2]] #y, x\nconf1 = location_preds['MODEL0'][:, 3]\nconf2 = location_preds['MODEL1'][:, 3]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:52:27.32264Z","iopub.execute_input":"2025-09-16T12:52:27.323136Z","iopub.status.idle":"2025-09-16T12:52:27.327491Z","shell.execute_reply.started":"2025-09-16T12:52:27.32311Z","shell.execute_reply":"2025-09-16T12:52:27.326639Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(12, 12))\n\n\n#label map\ncmap_base = plt.get_cmap(\"jet\", len(SEG_LABELS_TO_IDX))\ncolors = [cmap_base(i) for i in range(len(SEG_LABELS_TO_IDX))]\ncolors[0] = (0.25, 0.25, 0.25, 1.0)\ncmap_custom = mcolors.ListedColormap(colors)\n\n# Image in grayscale\nax.imshow(vol_.max(0), cmap=\"gray\", alpha=0.5)\n\n# Mask with modified cmap\nax.imshow(resized_mask.max(0), alpha=0.25, cmap=cmap_custom)\n\n\n# Normalize confidence [0,1] for colormap\nnorm = plt.Normalize(0, 1)\ncmap_conf = cm.get_cmap(\"plasma\")  # high contrast confidence heatmap\n\n# Overlay detections\nsc1 = ax.scatter(\n    loc_preds_2d_1[:, 1], loc_preds_2d_1[:, 0],\n    c=conf1, cmap=cmap_conf, norm=norm,\n    s=conf1 * 400, alpha=0.9, marker=\"o\",\n    edgecolors=\"k\", linewidths=0.3,\n    label=\"Fold0 predicted points\"\n)\n\nsc2 = ax.scatter(\n    loc_preds_2d_2[:, 1], loc_preds_2d_2[:, 0],\n    c=conf2, cmap=cmap_conf, norm=norm,\n    s=conf2 * 400, alpha=0.9, marker=\"^\",\n    linewidths=0.8,\n    label=\"Fold1 predicted points\"\n)\n\n# Colorbar axis (to the right of the plot)\ncbar = fig.colorbar(sc1, ax=ax, orientation=\"vertical\", fraction=0.046, pad=0.04)\ncbar.set_label(\"Confidence (0 = low, 1 = high)\", fontsize=12)\n\n# Legends\nunique_classes = np.unique(resized_mask)\nhandles = [\n    mpatches.Patch(color=cmap_custom(idx), label=label)\n    for label, idx in SEG_LABELS_TO_IDX.items()\n    if idx in unique_classes\n]\n\n# Predicted points legend\nscatter_legend = ax.legend(\n    markerscale=2, fontsize=12, handlelength=2, handletextpad=0.8,\n    loc=\"lower right\", title=\"Predicted Points\"\n)\nax.add_artist(scatter_legend)\n\n# Segmentation legend\nax.legend(\n    handles=handles,\n    loc='upper center', bbox_to_anchor=(0.5, 1.05),\n    ncol=3, fancybox=True, shadow=True\n)\n\nax.axis(\"off\")\nax.set_title(\"MIP view on segmentation prediction and predicted points\", fontsize=20, pad=40)\n\nplt.tight_layout(rect=[0, 0, 0.95, 1])  # leave space for colorbar on right\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:52:27.330892Z","iopub.execute_input":"2025-09-16T12:52:27.331154Z","iopub.status.idle":"2025-09-16T12:52:29.184381Z","shell.execute_reply.started":"2025-09-16T12:52:27.331129Z","shell.execute_reply":"2025-09-16T12:52:29.183721Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# for i in range(0, 100, 5):\n#     depth_1 = location_preds['MODEL0'][:, 0]\n#     depth_2 = location_preds['MODEL1'][:, 0]\n#     loc_slice_preds_1 = location_preds['MODEL0'][depth_1==i]\n#     loc_slice_preds_2 = location_preds['MODEL1'][depth_2==i]\n#     fig, ax = plt.subplots(ncols=3, figsize = (16, 16))\n#     ax[0].imshow(resized_mask[i])\n#     ax[0].set_title('predicted mask')\n#     ax[1].imshow(vol_[i])\n#     ax[1].scatter(loc_slice_preds_1[:, 2], loc_slice_preds_1[:, 1], color = 'red', label ='yolo-fold0 preds')\n#     ax[1].scatter(loc_slice_preds_2[:, 2], loc_slice_preds_2[:, 1], color = 'white', label ='yolo-fold1 preds')\n#     ax[1].legend()\n#     ax[1].set_title('vol from nii preprocess')\n#     ax[2].imshow(all_slices_vol[i])\n#     ax[2].set_title('vol from yolo preprocess')\n#     plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:52:41.139203Z","iopub.execute_input":"2025-09-16T12:52:41.139483Z","iopub.status.idle":"2025-09-16T12:52:43.494336Z","shell.execute_reply.started":"2025-09-16T12:52:41.139462Z","shell.execute_reply":"2025-09-16T12:52:43.493608Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Invalid Point Checking","metadata":{}},{"cell_type":"code","source":"def check_valid(yolo_loc_cls_preds, seg_loc_cls_preds):\n    return yolo2seg_lookup[yolo_loc_cls_preds]==seg_loc_cls_preds","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:52:45.456486Z","iopub.execute_input":"2025-09-16T12:52:45.457087Z","iopub.status.idle":"2025-09-16T12:52:45.460456Z","shell.execute_reply.started":"2025-09-16T12:52:45.457061Z","shell.execute_reply":"2025-09-16T12:52:45.459708Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"loc_cls_preds_1 = location_preds['MODEL0'][:, -2].astype('int32')\nloc_cls_preds_2 = location_preds['MODEL1'][:, -2].astype('int32')\n\nloc_preds_1  = location_preds['MODEL0'][:, [0, 1, 2]].astype('int32')\nloc_preds_2  = location_preds['MODEL1'][:, [0, 1, 2]].astype('int32')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:52:45.593138Z","iopub.execute_input":"2025-09-16T12:52:45.593364Z","iopub.status.idle":"2025-09-16T12:52:45.597935Z","shell.execute_reply.started":"2025-09-16T12:52:45.593348Z","shell.execute_reply":"2025-09-16T12:52:45.597162Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"seg_cls_preds_1 = resized_mask[loc_preds_1[:, 0], loc_preds_1[:, 1], loc_preds_1[:, 2]]\nseg_cls_preds_2 = resized_mask[loc_preds_2[:, 0], loc_preds_2[:, 1], loc_preds_2[:, 2]]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:52:45.746913Z","iopub.execute_input":"2025-09-16T12:52:45.747177Z","iopub.status.idle":"2025-09-16T12:52:45.751758Z","shell.execute_reply.started":"2025-09-16T12:52:45.747155Z","shell.execute_reply":"2025-09-16T12:52:45.750978Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"valid_mask_1 = check_valid(loc_cls_preds_1, seg_cls_preds_1)\nvalid_mask_2 = check_valid(loc_cls_preds_2, seg_cls_preds_2)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:52:46.68097Z","iopub.execute_input":"2025-09-16T12:52:46.681484Z","iopub.status.idle":"2025-09-16T12:52:46.685555Z","shell.execute_reply.started":"2025-09-16T12:52:46.681459Z","shell.execute_reply":"2025-09-16T12:52:46.684735Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print('yolo-fold0 invalid point percentage:', round(100 * np.sum(np.logical_not(valid_mask_1))/len(valid_mask_1), 2), '%')\nprint('yolo-fold1 invalid point percentage:', round(100 * np.sum(np.logical_not(valid_mask_2))/len(valid_mask_2), 2), '%')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-16T12:52:46.834275Z","iopub.execute_input":"2025-09-16T12:52:46.835048Z","iopub.status.idle":"2025-09-16T12:52:46.840115Z","shell.execute_reply.started":"2025-09-16T12:52:46.835013Z","shell.execute_reply":"2025-09-16T12:52:46.839389Z"}},"outputs":[],"execution_count":null}]}