{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":99552,"databundleVersionId":13190393,"sourceType":"competition"}],"dockerImageVersionId":31089,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"from skimage.filters import sobel, frangi\nfrom skimage.morphology import white_tophat, disk\nfrom scipy.ndimage import gaussian_gradient_magnitude\nfrom __future__ import annotations\nimport os\nimport ast\nimport math\nimport multiprocessing as mp\nfrom pathlib import Path\nimport sys\nfrom typing import List, Dict, Any, Tuple\n\nimport numpy as np\nimport pandas as pd\nimport pydicom\nimport cv2\nfrom sklearn.model_selection import StratifiedKFold\nimport matplotlib.pyplot as plt\n\ntry:\n    from scipy.ndimage import zoom as nd_zoom\nexcept ImportError:\n    raise SystemExit(\"scipy is required. Install via: pip install scipy\")\n    \n#current_dir = Path(__file__).parent\n#parent_dir = current_dir.parent\n#sys.path.insert(0, str(parent_dir))\n\nTARGET_DEPTH = 32\nTARGET_SIZE = 384\nHU_MIN = -1200.0\nHU_MAX = 4000.0\nSTORE_NORMALIZED = False  # Set True to revert to [0,1] scaling\n\n# Globals for worker processes\n\ndata_path = '/kaggle/input/rsna-intracranial-aneurysm-detection'","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-08-19T15:20:22.564697Z","iopub.execute_input":"2025-08-19T15:20:22.565051Z","iopub.status.idle":"2025-08-19T15:20:24.766408Z","shell.execute_reply.started":"2025-08-19T15:20:22.56502Z","shell.execute_reply":"2025-08-19T15:20:24.76551Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import SimpleITK as sitk\nfrom skimage.filters import frangi\nfrom skimage.morphology import skeletonize, remove_small_objects\nfrom scipy.ndimage import binary_closing, binary_opening, gaussian_filter","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-19T15:20:24.767584Z","iopub.execute_input":"2025-08-19T15:20:24.768119Z","iopub.status.idle":"2025-08-19T15:20:24.941097Z","shell.execute_reply.started":"2025-08-19T15:20:24.76809Z","shell.execute_reply":"2025-08-19T15:20:24.940276Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"data_path = '/kaggle/input/rsna-intracranial-aneurysm-detection'\nwindows = {\n    'CT': (40, 80),\n    'CTA': (50, 350),\n    'MRA': (600, 1200),\n    'MR': (600, 1200),\n    'MRI': (40, 80),\n}\n\nLABELS_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}\n\nIMG_SIZE = 512\nFACTOR = 1\nSEED = 42\nN_FOLDS = 5\nCORES = 4","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-19T15:20:24.941991Z","iopub.execute_input":"2025-08-19T15:20:24.942217Z","iopub.status.idle":"2025-08-19T15:20:24.950236Z","shell.execute_reply.started":"2025-08-19T15:20:24.942199Z","shell.execute_reply":"2025-08-19T15:20:24.949251Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def _load_series_dicom_paths(series_uid: str, root: Path) -> List[Path]:\n    series_dir = root / 'series' / series_uid\n    paths = []\n    for r, _, files in os.walk(series_dir):\n        for f in files:\n            if f.endswith('.dcm'):\n                paths.append(Path(r) / f)\n    return paths\n\n\ndef _read_dicom(path: Path):\n    ds = pydicom.dcmread(str(path), force=True)\n    arr = ds.pixel_array.astype(np.float32)\n    if hasattr(ds, 'RescaleSlope') and hasattr(ds, 'RescaleIntercept'):\n        arr = arr * float(ds.RescaleSlope) + float(ds.RescaleIntercept)\n    return ds, arr\n\n\ndef _extract_slice_position(ds) -> float:\n    # Prefer ImagePositionPatient z, fallback to InstanceNumber\n    if hasattr(ds, 'ImagePositionPatient') and len(ds.ImagePositionPatient) == 3:\n        try:\n            return float(ds.ImagePositionPatient[2])\n        except Exception:\n            pass\n    if hasattr(ds, 'InstanceNumber'):\n        try:\n            return float(ds.InstanceNumber)\n        except Exception:\n            pass\n    return 0.0\n\n\ndef _resample_depth(volume: np.ndarray, target_depth: int) -> np.ndarray:\n    if volume.shape[0] == target_depth:\n        return volume\n    depth_zoom = target_depth / volume.shape[0]\n    # zoom along depth only, order=1 linear\n    return nd_zoom(volume, (depth_zoom, 1.0, 1.0), order=1)\n\n\ndef _resize_inplane(volume: np.ndarray, target_hw: int) -> np.ndarray:\n    d, h, w = volume.shape\n    if h == target_hw and w == target_hw:\n        return volume\n    resized = np.empty((d, target_hw, target_hw), dtype=volume.dtype)\n    for i in range(d):\n        resized[i] = cv2.resize(volume[i], (target_hw, target_hw), interpolation=cv2.INTER_LINEAR)\n    return resized\n\n\ndef _clip_or_normalize(volume: np.ndarray) -> np.ndarray:\n    \"\"\"Either clip-only (raw HU retained in range) or clip+normalize to [0,1].\"\"\"\n    vol = np.clip(volume, HU_MIN, HU_MAX).astype(np.float32)\n    if STORE_NORMALIZED:\n        vol = (vol - HU_MIN) / (HU_MAX - HU_MIN)\n    return vol\n\n\ndef extract_details(vol):\n    mip = vol.max(axis=0, keepdims=True)\n    std_proj = vol.std(axis=0, keepdims=True)\n    edges = np.stack([sobel(slice_) for slice_ in vol], axis=0)\n    edge_proj = edges.max(axis=0, keepdims=True)\n    vesselness = np.stack([frangi(slice_) for slice_ in vol], axis=0)\n    vessel_proj = vesselness.max(axis=0, keepdims=True)\n    gradmag = np.stack([gaussian_gradient_magnitude(slice_, sigma=1) for slice_ in vol], axis=0)\n    grad_proj = gradmag.max(axis=0, keepdims=True)\n    extracted_vol = np.concatenate([vol, mip, std_proj, edge_proj, vessel_proj, grad_proj], axis=0) #(32 + 5, 384, 384)\n    return extracted_vol\n\n\ndef _process_single_series(uid: str, root: Path) -> Dict[str, Any]:\n    try:\n        dcm_paths = _load_series_dicom_paths(uid, root)\n        if not dcm_paths:\n            return {\"series_uid\": uid, \"volume_filename\": None, \"num_slices_raw\": 0}\n        slices: List[Tuple[float, np.ndarray]] = []\n        for p in dcm_paths:\n            try:\n                ds, arr = _read_dicom(p)\n                # If multi-frame (arr.ndim==3) stack frames individually\n                if arr.ndim == 3 and arr.shape[-1] != 3:\n                    for fi in range(arr.shape[0]):\n                        slices.append((_extract_slice_position(ds) + fi * 0.001, arr[fi].astype(np.float32)))\n                else:\n                    if arr.ndim == 3 and arr.shape[-1] == 3:\n                        # Convert RGB to grayscale\n                        arr = cv2.cvtColor(arr.astype(np.float32), cv2.COLOR_BGR2GRAY)\n                    slices.append((_extract_slice_position(ds), arr.astype(np.float32)))\n            except Exception:\n                continue\n        if not slices:\n            return {\"series_uid\": uid, \"volume_filename\": None, \"num_slices_raw\": 0}\n        # Sort by z\n        slices.sort(key=lambda x: x[0])\n        vol = np.stack([s[1] for s in slices], axis=0)  # (D, H, W)\n        num_raw = vol.shape[0]\n        # Clip HU range (optionally normalize based on STORE_NORMALIZED)\n        vol = _clip_or_normalize(vol)\n        # Depth resample\n        vol = _resample_depth(vol, TARGET_DEPTH)\n        # In-plane resize\n        vol = _resize_inplane(vol, TARGET_SIZE)\n        # Save\n        vol_filename = f\"{uid}_d{TARGET_DEPTH}_sz{TARGET_SIZE}.npz\"\n        # Save meta: [HU_MIN, HU_MAX, normalized_flag]\n        meta = np.array([HU_MIN, HU_MAX, 1.0 if STORE_NORMALIZED else 0.0], dtype=np.float32)\n        return {\"series_uid\": uid, \"volume_filename\": vol_filename, \"num_slices_raw\": num_raw, 'volume': vol, 'meta': meta}\n    except Exception as e:\n        return {\"series_uid\": uid, \"volume_filename\": None, \"error\": str(e), \"num_slices_raw\": 0}\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-19T15:20:24.952062Z","iopub.execute_input":"2025-08-19T15:20:24.952311Z","iopub.status.idle":"2025-08-19T15:20:24.975754Z","shell.execute_reply.started":"2025-08-19T15:20:24.95229Z","shell.execute_reply":"2025-08-19T15:20:24.974778Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"root = Path(data_path)\nprocessed = root / 'processed'\nvol_dir = processed / 'volumes_3d'\ntrain_df = pd.read_csv(root / 'train.csv')\nlabel_df = pd.read_csv(root / 'train_localizers.csv')\nmf_dicom_uids = pd.read_csv(root / 'multiframe_dicoms.csv') if (root / 'multiframe_dicoms.csv').exists() else pd.DataFrame(columns=['SeriesInstanceUID'])\n\nignore_uids = set([\n    '1.2.826.0.1.3680043.8.498.11145695452143851764832708867797988068',\n    '1.2.826.0.1.3680043.8.498.35204126697881966597435252550544407444',\n    '1.2.826.0.1.3680043.8.498.87480891990277582946346790136781912242',\n]) | set(mf_dicom_uids['SeriesInstanceUID'].tolist())\n\ntrain_df = train_df[~train_df['SeriesInstanceUID'].isin(ignore_uids)].reset_index(drop=True)\ntrain_df['fold_id'] = 0\n\nskf = StratifiedKFold(n_splits=N_FOLDS, random_state=SEED, shuffle=True)\nfor fold, (_, val_idx) in enumerate(skf.split(train_df['SeriesInstanceUID'], train_df['Aneurysm Present'])):\n    train_df.loc[val_idx, 'fold_id'] = fold\n\nuids = train_df['SeriesInstanceUID'].unique().tolist()\nprint(f\"Preparing 3D volumes for {len(uids)} series -> target shape ({TARGET_DEPTH}, {TARGET_SIZE}, {TARGET_SIZE})\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-19T15:20:24.976832Z","iopub.execute_input":"2025-08-19T15:20:24.97713Z","iopub.status.idle":"2025-08-19T15:20:25.089009Z","shell.execute_reply.started":"2025-08-19T15:20:24.977099Z","shell.execute_reply":"2025-08-19T15:20:25.088013Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"uid = uids[-1]\nseries_dict  = _process_single_series(uid, root)\nvol = series_dict['volume']","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-19T15:20:25.090126Z","iopub.execute_input":"2025-08-19T15:20:25.090459Z","iopub.status.idle":"2025-08-19T15:20:30.522246Z","shell.execute_reply.started":"2025-08-19T15:20:25.090432Z","shell.execute_reply":"2025-08-19T15:20:30.521232Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import SimpleITK as sitk","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-19T15:20:30.523215Z","iopub.execute_input":"2025-08-19T15:20:30.523542Z","iopub.status.idle":"2025-08-19T15:20:30.527649Z","shell.execute_reply.started":"2025-08-19T15:20:30.523505Z","shell.execute_reply":"2025-08-19T15:20:30.526845Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def vessel_enhancement(vol, sigmas=(1, 3)):\n    \"\"\"\n    Apply 3D Frangi filter to enhance vessel-like structures.\n    vol: (D, H, W) numpy array\n    \"\"\"\n    # skimage.frangi supports 3D directly if input is 3D\n    ves = frangi(vol, sigmas=sigmas, black_ridges=False)\n    ves = (ves - ves.min()) / (ves.max() - ves.min() + 1e-8)  # normalize [0,1]\n    return ves\n\n# ---------------------------\n# 2. Vessel skeletonization\n# ---------------------------\ndef vessel_skeleton(ves, thresh=0.2):\n    \"\"\"\n    Threshold vesselness map and apply 3D skeletonization.\n    \"\"\"\n    mask = ves > thresh\n    skel = skeletonize(mask.astype(np.uint8))\n    return mask, skel\n\n# ---------------------------\n# 3. Registration to atlas\n# ---------------------------\ndef register_to_atlas(atlas_np, moving_np, atlas_spacing=(1.0,1.0,1.0), moving_spacing=(1.0,1.0,1.0), do_bspline=True):\n    \"\"\"\n    Register moving volume to atlas using Rigid -> Affine -> (optional B-spline).\n    Returns: final_transform, fixed_atlas_img\n    \"\"\"\n\n    # Convert numpy → SimpleITK images\n    atlas_img = sitk.GetImageFromArray(atlas_np.astype(np.float32))\n    moving_img = sitk.GetImageFromArray(moving_np.astype(np.float32))\n    atlas_img.SetSpacing(atlas_spacing)\n    moving_img.SetSpacing(moving_spacing)\n\n    # --- Rigid registration ---\n    rigid_reg = sitk.ImageRegistrationMethod()\n    rigid_reg.SetMetricAsMattesMutualInformation(numberOfHistogramBins=32)\n    rigid_reg.SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=200)\n    rigid_reg.SetInterpolator(sitk.sitkLinear)\n    rigid_tx = sitk.VersorRigid3DTransform()\n    rigid_reg.SetInitialTransform(rigid_tx)\n    rigid_tx = rigid_reg.Execute(atlas_img, moving_img)\n\n    # --- Affine registration (fixed) ---\n    affine_reg = sitk.ImageRegistrationMethod()\n    affine_reg.SetMetricAsMattesMutualInformation(numberOfHistogramBins=32)\n    affine_reg.SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=200)\n    affine_reg.SetInterpolator(sitk.sitkLinear)\n\n    # Initialize affine from rigid\n    affine_tx = sitk.AffineTransform(3)\n    affine_tx.SetMatrix(rigid_tx.GetMatrix())\n    affine_tx.SetTranslation(rigid_tx.GetTranslation())\n\n    affine_reg.SetInitialTransform(affine_tx, inPlace=False)\n    affine_tx = affine_reg.Execute(atlas_img, moving_img)\n\n    final_tx = affine_tx\n\n    # --- Optional BSpline refinement ---\n    if do_bspline:\n        # Setup BSpline registration\n        grid_physical_spacing = [50.0, 50.0, 50.0]  # adjust to your volume size\n        image_spacing = atlas_img.GetSpacing()\n        image_size = atlas_img.GetSize()\n        mesh_size = [int(sz*spc/gsp + 0.5)\n                     for sz, spc, gsp in zip(image_size, image_spacing, grid_physical_spacing)]\n        bspline_tx = sitk.BSplineTransformInitializer(atlas_img, mesh_size)\n    \n        bspline_reg = sitk.ImageRegistrationMethod()\n        bspline_reg.SetMetricAsMattesMutualInformation(numberOfHistogramBins=32)\n        bspline_reg.SetOptimizerAsGradientDescent(learningRate=0.1,\n                                                  numberOfIterations=200,\n                                                  convergenceMinimumValue=1e-6,\n                                                  convergenceWindowSize=10)\n        bspline_reg.SetInterpolator(sitk.sitkLinear)\n        bspline_reg.SetInitialTransform(bspline_tx, inPlace=False)\n    \n        bspline_tx = bspline_reg.Execute(atlas_img, moving_img)\n    \n        # ✅ Use CompositeTransform to combine affine + bspline\n        final_tx = sitk.CompositeTransform(3)\n        final_tx.AddTransform(affine_tx)\n        final_tx.AddTransform(bspline_tx)\n    else:\n        final_tx = affine_tx\n\n    return final_tx, atlas_img","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-19T15:20:30.528635Z","iopub.execute_input":"2025-08-19T15:20:30.529078Z","iopub.status.idle":"2025-08-19T15:20:30.54752Z","shell.execute_reply.started":"2025-08-19T15:20:30.529054Z","shell.execute_reply":"2025-08-19T15:20:30.546452Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Step 1: Vessel enhancement\nves = vessel_enhancement(vol, sigmas=(1,2,3))\n\n# Step 2: Skeleton\nmask, skel = vessel_skeleton(ves, thresh=0.2)\n\n# Step 3: Registration (atlas must be provided)\natlas_vol = np.random.rand(32,384,384).astype(np.float32)\nfinal_tx, atlas_img = register_to_atlas(atlas_vol, vol)\n\n# Visualization\nproj_ves = ves.max(axis=0)\nproj_skel = skel.max(axis=0)\n\nplt.figure()\nplt.subplot(1,2,1); plt.imshow(proj_ves, cmap=\"gray\"); plt.title(\"Vesselness Projection\")\nplt.subplot(1,2,2); plt.imshow(proj_skel, cmap=\"gray\"); plt.title(\"Skeleton Projection\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-19T15:20:30.548334Z","iopub.execute_input":"2025-08-19T15:20:30.548601Z","execution_failed":"2025-08-19T15:20:58.948Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}