{"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":13190393,"sourceType":"competition"}],"dockerImageVersionId":31089,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"!pip install pyvista","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-07T15:44:23.972278Z","iopub.execute_input":"2025-08-07T15:44:23.972546Z","iopub.status.idle":"2025-08-07T15:44:28.976082Z","shell.execute_reply.started":"2025-08-07T15:44:23.972525Z","shell.execute_reply":"2025-08-07T15:44:28.975037Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nimport glob\nfrom pathlib import Path\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm\nfrom scipy.ndimage import binary_dilation\n\nimport pydicom\nimport nibabel as nib\nimport SimpleITK as sitk\nimport pyvista as pv\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-08-07T15:44:28.977768Z","iopub.execute_input":"2025-08-07T15:44:28.978036Z","iopub.status.idle":"2025-08-07T15:44:31.676338Z","shell.execute_reply.started":"2025-08-07T15:44:28.97801Z","shell.execute_reply":"2025-08-07T15:44:31.675763Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class config:\n    root_dicom_folder ='/kaggle/input/rsna-intracranial-aneurysm-detection/series/'\n    root_mask = '/kaggle/input/rsna-intracranial-aneurysm-detection/segmentations/'\n    ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-07T15:44:31.677054Z","iopub.execute_input":"2025-08-07T15:44:31.677387Z","iopub.status.idle":"2025-08-07T15:44:31.68165Z","shell.execute_reply.started":"2025-08-07T15:44:31.677368Z","shell.execute_reply":"2025-08-07T15:44:31.6807Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd\ntrain_localizer = pd.read_csv('/kaggle/input/rsna-intracranial-aneurysm-detection/train_localizers.csv')\ntrain_localizer['path_folder_dicom'] = train_localizer['SeriesInstanceUID'].apply(lambda uid: config.root_dicom_folder + uid )\ntrain_localizer['path_slice'] = train_localizer.apply(lambda row: config.root_dicom_folder + row.SeriesInstanceUID + '/' + row.SOPInstanceUID + '.dcm',axis=1)\ntrain_localizer['path_seg_mask'] = train_localizer['SeriesInstanceUID'].apply(lambda name: config.root_mask + name + '/' +  name + '_cowseg.nii')\n\n## only keep vessel mask\n\nlist_uid = os.listdir('/kaggle/input/rsna-intracranial-aneurysm-detection/segmentations')\ntrain_localizer = train_localizer[train_localizer['SeriesInstanceUID'].isin(list_uid)]\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-07T15:44:31.683348Z","iopub.execute_input":"2025-08-07T15:44:31.683529Z","iopub.status.idle":"2025-08-07T15:44:31.783255Z","shell.execute_reply.started":"2025-08-07T15:44:31.683513Z","shell.execute_reply":"2025-08-07T15:44:31.782686Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import SimpleITK as sitk\nimport numpy as np\nfrom pathlib import Path\n\ndef load_dicom_volume_sitk(series_dir):\n    \"\"\"\n    Load a DICOM series into a NumPy volume using SimpleITK.\n    \n    Args:\n        series_dir (str or Path): Path to folder containing a DICOM series.\n    \n    Returns:\n        volume (np.ndarray): 3D array (z, y, x) in float32.\n        spacing (tuple): Physical spacing between voxels (z, y, x) in mm.\n        origin (tuple): Physical origin (x, y, z) in mm.\n        direction (tuple): Direction cosines.\n    \"\"\"\n    series_dir = Path(series_dir)\n    reader = sitk.ImageSeriesReader()\n    \n    dicom_names = reader.GetGDCMSeriesFileNames(str(series_dir))\n    if not dicom_names:\n        raise FileNotFoundError(f\"No DICOM files found in {series_dir}\")\n    \n    reader.SetFileNames(dicom_names)\n    image = reader.Execute()\n    \n    # Convert to NumPy array (z, y, x)\n    # volume = sitk.GetArrayFromImage(image).astype(np.float32)\n    \n    return image","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-07T15:44:31.78396Z","iopub.execute_input":"2025-08-07T15:44:31.784147Z","iopub.status.idle":"2025-08-07T15:44:31.789211Z","shell.execute_reply.started":"2025-08-07T15:44:31.784131Z","shell.execute_reply":"2025-08-07T15:44:31.788367Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from pathlib import Path\nimport SimpleITK as sitk\nimport torch\nimport torch.nn.functional as F\nimport numpy as np\nfrom scipy.ndimage import generate_binary_structure\n\nimport numpy as np\nimport pydicom\nfrom pathlib import Path\nfrom scipy.ndimage import binary_dilation\n\ndef spherical_struct(radius):\n    \"\"\"\n    Create a spherical structuring element for 3D dilation.\n\n    Args:\n        radius (int): Radius in voxels.\n\n    Returns:\n        np.ndarray: 3D spherical structuring element (bool array).\n    \"\"\"\n    L = np.arange(-radius, radius + 1)\n    X, Y, Z = np.meshgrid(L, L, L, indexing=\"ij\")\n    return (X**2 + Y**2 + Z**2) <= radius**2\n\ndef build_mask_for_volume(path_dcm, path_slices, coords, dilation_radius=10):\n    \"\"\"\n    Build a 3D binary mask for a DICOM volume given labeled slice coordinates.\n    \n    Parameters:\n        path_dcm (str): Path to DICOM folder for the volume.\n        path_slices (list[str]): Paths to the DICOM slice files containing labels.\n        coords (list[tuple]): (x, y) coordinates for each labeled slice.\n        dilation_radius (int): Dilation radius in voxels.\n    \n    Returns:\n        np.ndarray: 3D mask (z, y, x)\n    \"\"\"\n    path_dcm = Path(path_dcm)\n    dicom_files = sorted(path_dcm.glob(\"*.dcm\"))\n    if not dicom_files:\n        raise FileNotFoundError(f\"No DICOM files found in {path_dcm}\")\n\n    # Read metadata for size\n    ds0 = pydicom.dcmread(str(dicom_files[0]), stop_before_pixels=True)\n    rows, cols = ds0.Rows, ds0.Columns\n    num_slices = len(dicom_files)\n\n    # Map SOPInstanceUID to slice index\n    sop_uid_to_index = {}\n    for idx, f in enumerate(dicom_files):\n        ds = pydicom.dcmread(str(f), stop_before_pixels=True)\n        sop_uid_to_index[ds.SOPInstanceUID] = idx\n\n    # Initialize empty mask\n    mask = np.zeros((num_slices, rows, cols), dtype=np.uint8)\n\n    # Assign points to mask\n    for path_slice, point in zip(path_slices, coords):\n        ds = pydicom.dcmread(str(path_slice), stop_before_pixels=True)\n        # print(ds.SOPInstanceUID,sop_uid_to_index)\n        sop_uid = ds.SOPInstanceUID\n        if sop_uid not in sop_uid_to_index:\n            continue\n        z_idx = sop_uid_to_index[sop_uid]\n\n        # Ensure point is tuple\n        if isinstance(point, str):\n            point = eval(point)\n        x, y = int(point[0]), int(point[1])\n        if 0 <= y < rows and 0 <= x < cols:\n            mask[z_idx, y, x] = 1\n        # mask[z_idx, ...] = dilate_mask_2d(mask[z_idx, ...],radius=dilation_radius)\n\n    # Dilate mask in 3D\n    if dilation_radius > 0:\n        from scipy.ndimage import generate_binary_structure\n        # struct = generate_binary_structure(3, 1)\n        struct = spherical_struct(dilation_radius)\n        mask = binary_dilation(mask, structure=struct, iterations=dilation_radius) #.astype(np.uint8)\n\n    return mask\n\nimport numpy as np\nimport plotly.graph_objects as go\n\ndef vis_3d_gpu(vol_3d: np.ndarray, mask_3d: np.ndarray, sample_rate: float = 1.0):\n    \"\"\"\n    GPU-accelerated 3D visualization using Plotly WebGL backend.\n    Shows volume and mask using voxel coordinates.\n    \n    Parameters:\n        vol_3d (np.ndarray): Binary vessel volume.\n        mask_3d (np.ndarray): Binary abnormal region mask.\n        sample_rate (float): Fraction of voxels to visualize for performance (0.1–0.5 is typical).\n    \"\"\"\n    assert vol_3d.shape == mask_3d.shape, \"Volume and mask must have the same shape\"\n\n    # Get voxel positions\n    vol_coords = np.argwhere(vol_3d > 0)\n    mask_coords = np.argwhere(mask_3d > 0)\n\n    # Optionally subsample for speed\n    if sample_rate < 1.0:\n        vol_coords = vol_coords[np.random.choice(len(vol_coords), int(len(vol_coords) * sample_rate), replace=False)]\n        mask_coords = mask_coords[np.random.choice(len(mask_coords), int(len(mask_coords) * sample_rate), replace=False)]\n\n    fig = go.Figure()\n\n    # Vessel volume points\n    if len(vol_coords) > 0:\n        fig.add_trace(go.Scatter3d(\n            x=vol_coords[:, 2], y=vol_coords[:, 1], z=vol_coords[:, 0],\n            mode='markers',\n            marker=dict(size=2, color='lightgray', opacity=0.2),\n            name='Vessel Volume'\n        ))\n\n    # Abnormal mask points\n    if len(mask_coords) > 0:\n        fig.add_trace(go.Scatter3d(\n            x=mask_coords[:, 2], y=mask_coords[:, 1], z=mask_coords[:, 0],\n            mode='markers',\n            marker=dict(size=2, color='red', opacity=0.7),\n            name='Abnormal Mask'\n        ))\n\n    fig.update_layout(\n        scene=dict(aspectmode='data'),\n        margin=dict(l=0, r=0, b=0, t=30),\n        title=\"GPU-accelerated 3D visualization (Plotly WebGL)\",\n        showlegend=True\n    )\n\n    fig.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-07T15:44:31.790061Z","iopub.execute_input":"2025-08-07T15:44:31.790524Z","iopub.status.idle":"2025-08-07T15:44:35.651067Z","shell.execute_reply.started":"2025-08-07T15:44:31.790507Z","shell.execute_reply":"2025-08-07T15:44:35.650299Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# seg_mask = nii.get_fdata()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-07T15:44:35.652178Z","iopub.execute_input":"2025-08-07T15:44:35.652593Z","iopub.status.idle":"2025-08-07T15:44:35.665827Z","shell.execute_reply.started":"2025-08-07T15:44:35.652564Z","shell.execute_reply":"2025-08-07T15:44:35.664499Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# seg_mask.shape","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-07T15:44:35.666236Z","iopub.status.idle":"2025-08-07T15:44:35.666447Z","shell.execute_reply.started":"2025-08-07T15:44:35.666344Z","shell.execute_reply":"2025-08-07T15:44:35.666353Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"mask_3d = build_mask_for_volume(path_dcm, path_slices, coords,dilation_radius=2)\nprint(mask_3d.shape)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-07T15:44:35.667634Z","iopub.status.idle":"2025-08-07T15:44:35.667893Z","shell.execute_reply.started":"2025-08-07T15:44:35.667781Z","shell.execute_reply":"2025-08-07T15:44:35.667792Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"mask_3d = mask_3d[::-1, :, :]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-07T15:44:35.669025Z","iopub.status.idle":"2025-08-07T15:44:35.669291Z","shell.execute_reply.started":"2025-08-07T15:44:35.669159Z","shell.execute_reply":"2025-08-07T15:44:35.669168Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"mask_3d.shape","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-07T15:44:35.670238Z","iopub.status.idle":"2025-08-07T15:44:35.670561Z","shell.execute_reply.started":"2025-08-07T15:44:35.670391Z","shell.execute_reply":"2025-08-07T15:44:35.670404Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":" seg_mask = np.transpose(seg_mask, (2,0,1))\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-07T15:44:35.671593Z","iopub.status.idle":"2025-08-07T15:44:35.67219Z","shell.execute_reply.started":"2025-08-07T15:44:35.672006Z","shell.execute_reply":"2025-08-07T15:44:35.672023Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"seg_mask.shape","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-07T15:44:35.673247Z","iopub.status.idle":"2025-08-07T15:44:35.673493Z","shell.execute_reply.started":"2025-08-07T15:44:35.67335Z","shell.execute_reply":"2025-08-07T15:44:35.673358Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":" vis_3d_gpu(seg_mask, mask_3d)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-07T15:44:35.674784Z","iopub.status.idle":"2025-08-07T15:44:35.675083Z","shell.execute_reply.started":"2025-08-07T15:44:35.674935Z","shell.execute_reply":"2025-08-07T15:44:35.674955Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def numpy_to_sitk_mask(mask_np, reference_image):\n    sitk_mask = sitk.GetImageFromArray(mask_np.astype(np.uint8))  # (z, y, x)\n    sitk_mask.CopyInformation(reference_image)\n    return sitk_mask","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-07T15:45:47.615499Z","iopub.execute_input":"2025-08-07T15:45:47.615829Z","iopub.status.idle":"2025-08-07T15:45:47.619953Z","shell.execute_reply.started":"2025-08-07T15:45:47.615793Z","shell.execute_reply":"2025-08-07T15:45:47.619222Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from tqdm import tqdm\n\nfor uid, sub in tqdm(train_localizer.groupby('path_folder_dicom')):\n    \n    coords = sub['coordinates'].tolist()\n    coords = [(int(eval(cood)['x']),int(eval(cood)['y'])) for  cood in coords]\n    \n    path_dcm = sub['path_folder_dicom'].iloc[0]\n    path_slices = sub['path_slice'].tolist()\n\n    # numpy 3d mask\n    mask_3d = build_mask_for_volume(path_dcm, path_slices, coords,dilation_radius=2)\n\n\n    volume_3d = load_dicom_volume_sitk(path_dcm)\n    \n    \n    # nii = nib.load(sub['path_seg_mask'].values[0])\n\n    seg_sitk = sitk.ReadImage(sub['path_seg_mask'].values[0])\n\n    sitk_mask_3d = numpy_to_sitk_mask(mask_3d, seg_sitk)\n    resampled_mask = sitk.Resample(sitk_mask_3d, seg_sitk, sitk.Transform(), sitk.sitkNearestNeighbor, 0)\n\n    vol_3d = sitk.Resample(volume_3d, seg_sitk, sitk.Transform(), sitk.sitkNearestNeighbor, 0)\n    volume_3d = sitk.GetArrayFromImage(volume_3d).astype(np.float32)\n    \n    # Convert back to numpy for visualization\n    abnormal_mask = sitk.GetArrayFromImage(resampled_mask)\n    \n    seg_mask = sitk.GetArrayFromImage(seg_sitk)\n    \n    # seg_mask = nii.get_fdata()\n    # seg_mask = np.transpose(seg_mask, (2,1,0))\n    seg_mask = seg_mask[:, ::-1, :]\n\n\n    # vis_3d_gpu(seg_mask.astype(np.uint8), volume_3d.astype(np.uint8))\n\n    vis_3d_gpu(seg_mask, abnormal_mask)\n    \n    print('done')\n    break","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"volume_3d = (volume_3d-volume_3d.min())/(volume_3d.max()-volume_3d.min())*255","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-07T15:46:01.580795Z","iopub.execute_input":"2025-08-07T15:46:01.581004Z","iopub.status.idle":"2025-08-07T15:46:01.833077Z","shell.execute_reply.started":"2025-08-07T15:46:01.580988Z","shell.execute_reply":"2025-08-07T15:46:01.832186Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"vis_3d_gpu(seg_mask, volume_3d)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-07T15:47:59.024258Z","iopub.execute_input":"2025-08-07T15:47:59.02501Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}