{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"In this notebook, you can find a complete set of customizable functions for:\n* Loading NIFTI and DICOM files into a 3D numpy array. \n* Plotting a range of slices of the 3D image from any angle (Coronal / Sagittal / Axial):\n    * With a mask overlay (segmentation mask)\n    * Without masks\n    \n","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nimport os\nimport random\nimport re\n\nfrom tqdm import tqdm\n\nimport pydicom as dicom\nimport nibabel as nib\n","metadata":{"execution":{"iopub.status.busy":"2023-08-05T07:03:47.896396Z","iopub.execute_input":"2023-08-05T07:03:47.896782Z","iopub.status.idle":"2023-08-05T07:03:47.903165Z","shell.execute_reply.started":"2023-08-05T07:03:47.896751Z","shell.execute_reply":"2023-08-05T07:03:47.901988Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Define hyperparameters","metadata":{}},{"cell_type":"code","source":"DS_RATE = 2","metadata":{"execution":{"iopub.status.busy":"2023-08-05T07:03:49.511428Z","iopub.execute_input":"2023-08-05T07:03:49.511822Z","iopub.status.idle":"2023-08-05T07:03:49.516901Z","shell.execute_reply.started":"2023-08-05T07:03:49.51179Z","shell.execute_reply":"2023-08-05T07:03:49.515725Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Define functions to load slice images (DCM or NII format) to contruct 3D images","metadata":{}},{"cell_type":"code","source":"def create_3D_scans(folder, downsample_rate=1): \n    filenames = os.listdir(folder)\n    filenames = [int(filename.split('.')[0]) for filename in filenames]\n    filenames = sorted(filenames)\n    filenames = [str(filename) + '.dcm' for filename in filenames]\n        \n    volume = []\n    for filename in tqdm(filenames[::downsample_rate]):\n        filepath = os.path.join(folder, filename)\n        ds = dicom.dcmread(filepath)\n        image = ds.pixel_array\n        \n        # find rescale params\n        if (\"RescaleIntercept\" in ds) and (\"RescaleSlope\" in ds):\n            intercept = float(ds.RescaleIntercept)\n            slope = float(ds.RescaleSlope)\n    \n        # find clipping params\n        center = int(ds.WindowCenter)\n        width = int(ds.WindowWidth)\n        low = center - width / 2\n        high = center + width / 2    \n        \n        \n        image = (image * slope) + intercept\n        image = np.clip(image, low, high)\n\n        image = (image / np.max(image) * 255).astype(np.int16)\n        image = image[::downsample_rate, ::downsample_rate]\n        volume.append( image )\n    \n    volume = np.stack(volume, axis=0)\n    return volume\n\n\ndef create_3D_segmentations(filepath, downsample_rate=1):\n    img = nib.load(filepath).get_fdata()\n    img = np.transpose(img, [1, 0, 2])\n    img = np.rot90(img, 1, (1,2))\n    img = img[::-1,:,:]\n    img = np.transpose(img, [1, 0, 2])\n    img = img[::downsample_rate, ::downsample_rate, ::downsample_rate]\n    return img\n\n\n\nfilepath = '/kaggle/input/rsna-2023-abdominal-trauma-detection/segmentations/21057.nii'\nvolume_seg = create_3D_segmentations(filepath, downsample_rate=DS_RATE)\nprint(f'3D segmentation file shape: {volume_seg.shape}')\n\nfilepath = '/kaggle/input/rsna-2023-abdominal-trauma-detection/train_images/10004/21057'\nvolume = create_3D_scans(filepath, downsample_rate=DS_RATE)\nprint(f'3D Image file shape: {volume.shape}')","metadata":{"execution":{"iopub.status.busy":"2023-08-05T07:03:50.307323Z","iopub.execute_input":"2023-08-05T07:03:50.307699Z","iopub.status.idle":"2023-08-05T07:04:08.355117Z","shell.execute_reply.started":"2023-08-05T07:03:50.30767Z","shell.execute_reply":"2023-08-05T07:04:08.353984Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Define function that takes 3D images and plots slices from a specific angle (e.g. Sagittal)","metadata":{}},{"cell_type":"code","source":"def plot_image_with_seg(volume, volume_seg=[], orientation='Coronal', num_subplots=20):\n    # simply copy\n    if len(volume_seg) == 0:\n        plot_mask = 0\n    else:\n        plot_mask = 1\n        \n    if orientation == 'Coronal':\n        slices = np.linspace(0, volume.shape[2]-1, num_subplots).astype(np.int16)\n        volume = volume.transpose([1, 0, 2])\n        if plot_mask:\n            volume_seg = volume_seg.transpose([1, 0, 2])\n        \n    elif orientation == 'Sagittal':\n        slices = np.linspace(0, volume.shape[2]-1, num_subplots).astype(np.int16)\n        volume = volume.transpose([2, 0, 1])\n        if plot_mask:\n            volume_seg = volume_seg.transpose([2, 0, 1])\n\n    elif orientation == 'Axial':\n        slices = np.linspace(0, volume.shape[0]-1, num_subplots).astype(np.int16)\n           \n    rows = np.max( [np.floor(np.sqrt(num_subplots)).astype(int) - 2, 1])\n    cols = np.ceil(num_subplots/rows).astype(int)\n    \n    fig, ax = plt.subplots(rows, cols, figsize=(cols * 2, rows * 4))\n    fig.tight_layout(h_pad=0.01, w_pad=0)\n    \n    ax = ax.ravel()\n    for this_ax in ax:\n        this_ax.axis('off')\n\n    for counter, this_slice in enumerate( slices ):\n        plt.sca(ax[counter])\n        \n        image = volume[this_slice, :, :]\n        plt.imshow(image, cmap='gray')\n        \n        if plot_mask:\n            mask = np.where(volume_seg[this_slice, :, :], volume_seg[this_slice, :, :], np.nan)\n            plt.imshow(mask, cmap='Set1', alpha=0.5)        \n        \n        \n        \n        \nplot_image_with_seg(volume, volume_seg, orientation='Coronal', num_subplots=10)\nplot_image_with_seg(volume, volume_seg, orientation='Sagittal', num_subplots=10)\nplot_image_with_seg(volume, volume_seg, orientation='Axial', num_subplots=10)","metadata":{"execution":{"iopub.status.busy":"2023-08-05T07:04:24.807847Z","iopub.execute_input":"2023-08-05T07:04:24.808294Z","iopub.status.idle":"2023-08-05T07:04:28.836758Z","shell.execute_reply.started":"2023-08-05T07:04:24.808259Z","shell.execute_reply":"2023-08-05T07:04:28.835583Z"},"trusted":true},"execution_count":null,"outputs":[]}]}