{"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":"# EDA\n\n\nIn this notebook, we perform an EDA on the data.\nMore specifically, we try to assess what is remarkable in the scan associated to each type of injury. \n\nTo make this notebook, a variety of others were used :\n* https://www.kaggle.com/code/aritrag/kerascv-starter-notebook-train\n* https://www.kaggle.com/code/aritrag/eda-train-csv\n* https://www.kaggle.com/code/ayushs9020/understanding-the-competition-rsna#2-|-Visualization-%F0%9F%94%AC\n* https://www.kaggle.com/code/bvinning/interactive-viewing-of-scans-and-segmentations\n* https://www.kaggle.com/code/parhammostame/construct-3d-arrays-from-dcm-nii-3-view-angles#Define-function-that-takes-3D-images-and-plots-slices-from-a-specific-angle-(e.g.-Sagittal)\n\nTo their authors I say thank you !\n\n---\n---\n# Generalities\n\nBefore that, we need to know some basic information, about the number of patients, scans, injuries we are dealing with.\n\n---\n## Counts\n\nThe number of patients, and the number of patients with each of the injuries is important for future work, to know if we need to balance anything.","metadata":{}},{"cell_type":"code","source":"! pip install -U dicomsdl","metadata":{"execution":{"iopub.status.busy":"2023-09-21T15:24:03.722575Z","iopub.execute_input":"2023-09-21T15:24:03.722972Z","iopub.status.idle":"2023-09-21T15:24:19.843802Z","shell.execute_reply.started":"2023-09-21T15:24:03.722941Z","shell.execute_reply":"2023-09-21T15:24:19.842781Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport cv2\nimport pydicom\nimport dicomsdl\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport matplotlib.animation as animation\nimport ipywidgets as widgets\nimport nibabel as nib\nimport numpy as np\n\nfrom IPython.display import HTML\nfrom tqdm import tqdm","metadata":{"execution":{"iopub.status.busy":"2023-09-21T15:24:19.845841Z","iopub.execute_input":"2023-09-21T15:24:19.846253Z","iopub.status.idle":"2023-09-21T15:24:21.300679Z","shell.execute_reply.started":"2023-09-21T15:24:19.846218Z","shell.execute_reply":"2023-09-21T15:24:21.299393Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"BASE_PATH = \"/kaggle/input/rsna-2023-abdominal-trauma-detection\"\nTRAIN_CSV = f\"{BASE_PATH}/train.csv\"","metadata":{"execution":{"iopub.status.busy":"2023-09-21T15:24:21.302109Z","iopub.execute_input":"2023-09-21T15:24:21.302669Z","iopub.status.idle":"2023-09-21T15:24:21.308596Z","shell.execute_reply.started":"2023-09-21T15:24:21.302637Z","shell.execute_reply":"2023-09-21T15:24:21.307305Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data = pd.read_csv(TRAIN_CSV)\ndata.head()","metadata":{"execution":{"iopub.status.busy":"2023-09-21T15:24:21.311276Z","iopub.execute_input":"2023-09-21T15:24:21.311672Z","iopub.status.idle":"2023-09-21T15:24:21.375568Z","shell.execute_reply.started":"2023-09-21T15:24:21.311634Z","shell.execute_reply":"2023-09-21T15:24:21.374391Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"nb_patients = data['patient_id'].count()\nnb_patients","metadata":{"execution":{"iopub.status.busy":"2023-09-21T08:50:20.176602Z","iopub.execute_input":"2023-09-21T08:50:20.177098Z","iopub.status.idle":"2023-09-21T08:50:20.188451Z","shell.execute_reply.started":"2023-09-21T08:50:20.177055Z","shell.execute_reply":"2023-09-21T08:50:20.186959Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data.loc[:, data.columns != 'patient_id'].apply(lambda x: round(100*x.sum()/nb_patients, 2))\n# another way would be to use data.describe()","metadata":{"execution":{"iopub.status.busy":"2023-09-21T08:50:20.190086Z","iopub.execute_input":"2023-09-21T08:50:20.190785Z","iopub.status.idle":"2023-09-21T08:50:20.210461Z","shell.execute_reply.started":"2023-09-21T08:50:20.190743Z","shell.execute_reply":"2023-09-21T08:50:20.20931Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Data is quite imbalanced with a majority of healthy people.\n\n---\n---\n# Visualize \n\nWe have to look for the patients with specific injuries and check their CT.\n\nLet's first open a CT to check how to visualize it.\n\n---\n## CT Scans","metadata":{}},{"cell_type":"code","source":"image_file = \"/kaggle/input/rsna-2023-abdominal-trauma-detection/train_images/10004/21057/1001.dcm\"\nds = pydicom.read_file(image_file)\nds","metadata":{"execution":{"iopub.status.busy":"2023-09-20T15:40:30.295168Z","iopub.execute_input":"2023-09-20T15:40:30.295939Z","iopub.status.idle":"2023-09-20T15:40:30.323993Z","shell.execute_reply.started":"2023-09-20T15:40:30.295895Z","shell.execute_reply":"2023-09-20T15:40:30.322706Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"DICOM is the standard medical format for images. As you can see it contains a variety of metadata.\nOne interests us in particular, \"Pixel Data\" countains a numpy ndarray.","metadata":{}},{"cell_type":"code","source":"plt.imshow(ds.pixel_array, cmap=plt.cm.bone)","metadata":{"execution":{"iopub.status.busy":"2023-09-20T15:40:30.325642Z","iopub.execute_input":"2023-09-20T15:40:30.325977Z","iopub.status.idle":"2023-09-20T15:40:30.803092Z","shell.execute_reply.started":"2023-09-20T15:40:30.325949Z","shell.execute_reply":"2023-09-20T15:40:30.801741Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There is actually a lot of images in a patient scan. According to the file tree there are 1022 dcm in scan 21057 of patient 10004 and 1044 in scan 51033. Hence it is not possible for me to go through all of them and determine what is wrong.\n\nMaybe other files could be useful, from the competition data documentation we know these :\n\n    image_level_labels.csv Train only. Identifies specific images that contain either bowel or extravasation injuries.\n\n        patient_id - A unique ID code for each patient.\n        series_id - A unique ID code for each scan.\n        instance_number - The image number within the scan. The lowest instance number for many series is above zero as the original scans were cropped to the abdomen.\n        injury_name - The type of injury visible in the frame.\n\n    segmentations/ Model generated pixel-level annotations of the relevant organs and some major bones for a subset of the scans in the training set. This data is provided in the nifti file format. The filenames are series IDs. You can find a description of the source model (total segmentator) here and the data used to train that model here.\n    \nWe start with the image_level_labels.csv, as it is supposed to contain image level information. ","metadata":{}},{"cell_type":"code","source":"IMAGE_LEVEL_LABEL_CSV = f\"{BASE_PATH}/image_level_labels.csv\"\nimage_label = pd.read_csv(IMAGE_LEVEL_LABEL_CSV)\nimage_label.head(), image_label.dtypes","metadata":{"execution":{"iopub.status.busy":"2023-09-20T15:40:30.804401Z","iopub.execute_input":"2023-09-20T15:40:30.804822Z","iopub.status.idle":"2023-09-20T15:40:30.830559Z","shell.execute_reply.started":"2023-09-20T15:40:30.804783Z","shell.execute_reply":"2023-09-20T15:40:30.82938Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"For patient 10004, and its scan 21057, we have extravasation injury that should be visible. ","metadata":{}},{"cell_type":"code","source":"image_label[(image_label['patient_id'] == 10004) & (image_label['series_id'] == 21057)]","metadata":{"execution":{"iopub.status.busy":"2023-09-20T15:40:30.832398Z","iopub.execute_input":"2023-09-20T15:40:30.833181Z","iopub.status.idle":"2023-09-20T15:40:30.854751Z","shell.execute_reply.started":"2023-09-20T15:40:30.833108Z","shell.execute_reply":"2023-09-20T15:40:30.853547Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"From 362 to 416, there are active extravasation according to this table. According to this [definition from Hôpitaux Universitaires de Genève](https://www.hug.ch/procedures-de-soins/extravasation-medicaments-non-cytotoxiques), extravasation is a leakage of drug in the body. An animation with a slider is the best way to see all the scans.","metadata":{}},{"cell_type":"code","source":"image_file = \"/kaggle/input/rsna-2023-abdominal-trauma-detection/train_images/10004/21057/1001.dcm\"\nds = pydicom.read_file(image_file)","metadata":{"execution":{"iopub.status.busy":"2023-09-20T15:40:30.860919Z","iopub.execute_input":"2023-09-20T15:40:30.861455Z","iopub.status.idle":"2023-09-20T15:40:30.871343Z","shell.execute_reply.started":"2023-09-20T15:40:30.861418Z","shell.execute_reply":"2023-09-20T15:40:30.870069Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Get a list of image file paths (replace 'image_folder' with your folder path)\nimage_folder = f'{BASE_PATH}/train_images/10004/21057'\nimages = [[filename[:-4], pydicom.read_file(os.path.join(image_folder, filename)).pixel_array] for filename in os.listdir(image_folder) if filename.endswith('.dcm') and 362 <= int(filename[:-4]) <= 416]\nimages = sorted(images, key=lambda x:x[0])\n\nfig, ax = plt.subplots()\nim = ax.imshow(images[0][1], cmap=plt.cm.bone)\n\nupdate = lambda i : im.set_array(images[i][1])\n\nani = animation.FuncAnimation(fig, update, frames=range(len(images)), repeat=True)\n\nHTML(ani.to_jshtml())","metadata":{"execution":{"iopub.status.busy":"2023-09-20T15:40:30.873345Z","iopub.execute_input":"2023-09-20T15:40:30.874159Z","iopub.status.idle":"2023-09-20T15:40:42.357706Z","shell.execute_reply.started":"2023-09-20T15:40:30.874092Z","shell.execute_reply":"2023-09-20T15:40:42.356526Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# save animation as gif\noutput_video_filename = '/kaggle/working/dicom_scan.gif'\nani.save(output_video_filename, writer='pillow', fps=10)  # You can adjust the 'fps' parameter as needed","metadata":{"execution":{"iopub.status.busy":"2023-09-20T15:40:42.359221Z","iopub.execute_input":"2023-09-20T15:40:42.359544Z","iopub.status.idle":"2023-09-20T15:40:48.168838Z","shell.execute_reply.started":"2023-09-20T15:40:42.359515Z","shell.execute_reply":"2023-09-20T15:40:48.167645Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Nice we are now going through the body of patient 10004. However we still do not know what is to be seen in that scan. The segmentations mentionned earlier should contain more information.\n\n---\n## Segmentation","metadata":{}},{"cell_type":"code","source":"try:\n    nifti_file = f'{BASE_PATH}/segmentations/21057.nii'\n    img = nib.load(nifti_file)\nexcept FileNotFoundError:\n    print(\"Scan 21057 of patient 10004  is not in the segmentations files\")","metadata":{"execution":{"iopub.status.busy":"2023-09-20T15:40:48.170393Z","iopub.execute_input":"2023-09-20T15:40:48.170701Z","iopub.status.idle":"2023-09-20T15:40:48.190579Z","shell.execute_reply.started":"2023-09-20T15:40:48.170675Z","shell.execute_reply":"2023-09-20T15:40:48.189418Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"nifti_file = f'{BASE_PATH}/segmentations/21057.nii'\nimg = nib.load(nifti_file)\nsegmentation_data = img.get_fdata()\nheader = img.header\n\n# Display the shape of the data\nprint(\"Data shape:\", segmentation_data.shape)\n\n# Display header information (metadata)\nprint(header)","metadata":{"execution":{"iopub.status.busy":"2023-09-20T15:40:48.192427Z","iopub.execute_input":"2023-09-20T15:40:48.193602Z","iopub.status.idle":"2023-09-20T15:40:56.651242Z","shell.execute_reply.started":"2023-09-20T15:40:48.193562Z","shell.execute_reply":"2023-09-20T15:40:56.650041Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"A description of entries is available [here](https://brainder.org/2012/09/23/the-nifti-file-format/). There length of the third dimension is 1022, which corresponds to the number of .dcm files. ","metadata":{}},{"cell_type":"code","source":"segmentation_data.shape","metadata":{"execution":{"iopub.status.busy":"2023-09-20T15:40:56.652976Z","iopub.execute_input":"2023-09-20T15:40:56.653692Z","iopub.status.idle":"2023-09-20T15:40:56.661028Z","shell.execute_reply.started":"2023-09-20T15:40:56.653652Z","shell.execute_reply":"2023-09-20T15:40:56.659981Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.imshow(segmentation_data[:,:,362].T, cmap=\"Set1\", origin=\"lower\") ","metadata":{"execution":{"iopub.status.busy":"2023-09-20T15:40:56.662313Z","iopub.execute_input":"2023-09-20T15:40:56.663121Z","iopub.status.idle":"2023-09-20T15:40:56.99551Z","shell.execute_reply.started":"2023-09-20T15:40:56.663078Z","shell.execute_reply":"2023-09-20T15:40:56.994187Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots()\nim = ax.imshow(segmentation_data[:,:,362].T, cmap=\"Set1\", origin=\"lower\")\n\nupdate = lambda i : im.set_array(segmentation_data[:,:,i].T)\n\nani = animation.FuncAnimation(fig, update, frames=range(362, 416), repeat=True)\n\nHTML(ani.to_jshtml())","metadata":{"execution":{"iopub.status.busy":"2023-09-20T15:40:56.997553Z","iopub.execute_input":"2023-09-20T15:40:56.997994Z","iopub.status.idle":"2023-09-20T15:41:03.855579Z","shell.execute_reply.started":"2023-09-20T15:40:56.997954Z","shell.execute_reply":"2023-09-20T15:41:03.854234Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Not completely satisfying yet as we do not know how to interpret the shapes we are seeing.\nAfter going through some notebooks I found some answers [here](https://www.kaggle.com/code/bvinning/interactive-viewing-of-scans-and-segmentations/comments) and [there](https://www.kaggle.com/competitions/rsna-2023-abdominal-trauma-detection/discussion/428538).\n\n    The segmentation labels are as follows:\n    i. 0 = background\n    ii. 1 = liver\n    iii. 2 = spleen\n    iv. 3 = left kidney\n    v. 4 = right kidney\n    vi. 5 = bowel ","metadata":{}},{"cell_type":"code","source":"np.unique(segmentation_data[:,:,362])","metadata":{"execution":{"iopub.status.busy":"2023-09-20T15:41:03.8569Z","iopub.execute_input":"2023-09-20T15:41:03.857256Z","iopub.status.idle":"2023-09-20T15:41:03.872771Z","shell.execute_reply.started":"2023-09-20T15:41:03.857226Z","shell.execute_reply":"2023-09-20T15:41:03.871431Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There are bowel in the segmentation. To do the overlay with the CT scan and the segmentation is not trivial. I rely on [this notebook](https://www.kaggle.com/code/parhammostame/construct-3d-arrays-from-dcm-nii-3-view-angles#Define-function-that-takes-3D-images-and-plots-slices-from-a-specific-angle-(e.g.-Sagittal)) to do it.\n\nA preliminary step for the segmentation voxels is to correct their orientation to match the one of the CT scan. Think of it as if for CT scan the patient was laying on a table and we scan him/her from head to toe. However for the segmentation, the patient is standing up and you're observing him/her from the bottom to the top. TODO recheck that not sure","metadata":{}},{"cell_type":"code","source":"img = np.transpose(segmentation_data, [1, 0, 2]) # invert 2 axes of each image\nimg = np.rot90(img, 1, (1,2)) # , forget about images, imagine a rotation of a 3D object made of voxels\nimg = img[::-1,:,:] # mirror\nimg = np.transpose(img, [1, 0, 2]) # retour aux axes originaux","metadata":{"execution":{"iopub.status.busy":"2023-09-20T15:41:03.874196Z","iopub.execute_input":"2023-09-20T15:41:03.874576Z","iopub.status.idle":"2023-09-20T15:41:03.885505Z","shell.execute_reply.started":"2023-09-20T15:41:03.874542Z","shell.execute_reply":"2023-09-20T15:41:03.884156Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots()\nim = ax.imshow(img[:,:,362].T, cmap=\"Set1\", origin=\"lower\")\n\nupdate = lambda i : im.set_array(img[:,:,i].T)\n\nani = animation.FuncAnimation(fig, update, frames=range(362, 416), repeat=True)\n\nHTML(ani.to_jshtml())","metadata":{"execution":{"iopub.status.busy":"2023-09-20T15:41:03.887221Z","iopub.execute_input":"2023-09-20T15:41:03.887614Z","iopub.status.idle":"2023-09-20T15:41:12.319154Z","shell.execute_reply.started":"2023-09-20T15:41:03.887579Z","shell.execute_reply":"2023-09-20T15:41:12.317747Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n# Overlay CT scans and segmentation\n\nThat's a lot better, apparently there is some work to do back on the CT scans. \nRecall that we had the following fields in our dicom dataset :\n\n    (0028, 1050) Window Center                       DS: '50.0'\n    (0028, 1051) Window Width                        DS: '400.0'\n    (0028, 1052) Rescale Intercept                   DS: '-1024.0'\n    (0028, 1053) Rescale Slope                       DS: '1.0'\n    (0028, 1054) Rescale Type                        LO: 'HU'\n\nThe Rescale Type indicates that there is rescaling to perform with Rescale Intercept and Rescale Slope. This rescaling is linked to the [Hounsfield scale](https://en.wikipedia.org/wiki/Hounsfield_scale), it enables us to pay attention to specific parts of the body. ","metadata":{}},{"cell_type":"code","source":"def clip_rescale_extract_dicom_image(dicom_ds):\n    image = dicom_ds.pixel_array\n    \n    # find rescale params\n    if (\"RescaleIntercept\" in dicom_ds) and (\"RescaleSlope\" in dicom_ds):\n        intercept = float(dicom_ds.RescaleIntercept)\n        slope = float(dicom_ds.RescaleSlope)\n\n    # find clipping params\n    center = int(dicom_ds.WindowCenter)\n    width = int(dicom_ds.WindowWidth)\n    low = center - width / 2\n    high = center + width / 2    \n\n    image = (image * slope) + intercept\n    image = np.clip(image, low, high)\n    image = (image / np.max(image) * 255).astype(np.int16)\n    return image\n\n# Get a list of image file paths (replace 'image_folder' with your folder path)\nimage_folder = f'{BASE_PATH}/train_images/10004/21057'\nimages = [[filename[:-4], clip_rescale_extract_dicom_image(pydicom.read_file(os.path.join(image_folder, filename)))] for filename in os.listdir(image_folder) if filename.endswith('.dcm') and 362 <= int(filename[:-4]) <= 416]\nimages = sorted(images, key=lambda x:x[0])\n\nfig, ax = plt.subplots()\nim = ax.imshow(images[0][1], cmap=plt.cm.bone)\n\nupdate = lambda i : im.set_array(images[i][1])\n\nani = animation.FuncAnimation(fig, update, frames=range(len(images)), repeat=True)\n\nHTML(ani.to_jshtml())","metadata":{"execution":{"iopub.status.busy":"2023-09-20T15:41:12.321102Z","iopub.execute_input":"2023-09-20T15:41:12.321838Z","iopub.status.idle":"2023-09-20T15:41:21.356659Z","shell.execute_reply.started":"2023-09-20T15:41:12.321778Z","shell.execute_reply":"2023-09-20T15:41:21.355165Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# save animation as gif\noutput_video_filename = '/kaggle/working/dicom_scan_hu.gif'\nani.save(output_video_filename, writer='pillow', fps=10)  # You can adjust the 'fps' parameter as needed","metadata":{"execution":{"iopub.status.busy":"2023-09-20T15:41:21.358361Z","iopub.execute_input":"2023-09-20T15:41:21.358738Z","iopub.status.idle":"2023-09-20T15:41:27.23077Z","shell.execute_reply.started":"2023-09-20T15:41:21.358704Z","shell.execute_reply":"2023-09-20T15:41:27.22952Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Reusing the very useful functions from Parham Mostame on https://www.kaggle.com/code/parhammostame/construct-3d-arrays-from-dcm-nii-3-view-angles#Define-function-that-takes-3D-images-and-plots-slices-from-a-specific-angle-(e.g.-Sagittal)\nDS_RATE = 3\n\ndef 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 = pydicom.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-09-20T15:41:27.232524Z","iopub.execute_input":"2023-09-20T15:41:27.233008Z","iopub.status.idle":"2023-09-20T15:41:37.096892Z","shell.execute_reply.started":"2023-09-20T15:41:27.232973Z","shell.execute_reply":"2023-09-20T15:41:37.09569Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_image_with_seg(volume, volume_seg=[], orientation='Coronal'):\n    if orientation == 'Coronal':\n        volume = volume.transpose([1, 0, 2])\n        volume_seg = volume_seg.transpose([1, 0, 2])\n    elif orientation == 'Sagittal':\n        volume = volume.transpose([2, 0, 1])\n        volume_seg = volume_seg.transpose([2, 0, 1])\n    elif orientation == 'Axial':\n        pass\n    else:\n        print(\"Orientation is either 'Axial', 'Coronal' or 'Sagittal'\")\n        \n    fig, ax = plt.subplots()\n    im = ax.imshow(volume[0, :, :], cmap=plt.cm.bone)\n    mask = np.where(volume_seg[0, :, :], volume_seg[0, :, :], np.nan)\n    seg = ax.imshow(mask, cmap='Set1', alpha=0.5)\n    \n    def update(i):\n        im.set_array(volume[i, :, :])\n        mask = np.where(volume_seg[i, :, :], volume_seg[i, :, :], np.nan)\n        seg.set_array(mask)\n\n    ani = animation.FuncAnimation(fig, update, frames=range(len(volume)), repeat=True)\n    return ani","metadata":{"execution":{"iopub.status.busy":"2023-09-20T15:41:37.098363Z","iopub.execute_input":"2023-09-20T15:41:37.098687Z","iopub.status.idle":"2023-09-20T15:41:37.109302Z","shell.execute_reply.started":"2023-09-20T15:41:37.098659Z","shell.execute_reply":"2023-09-20T15:41:37.108348Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dicom_scan_seg = plot_image_with_seg(volume, volume_seg, orientation='Axial')\nHTML(dicom_scan_seg.to_jshtml())","metadata":{"execution":{"iopub.status.busy":"2023-09-20T15:41:37.110327Z","iopub.execute_input":"2023-09-20T15:41:37.11068Z","iopub.status.idle":"2023-09-20T15:42:12.890071Z","shell.execute_reply.started":"2023-09-20T15:41:37.110652Z","shell.execute_reply":"2023-09-20T15:42:12.888791Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# save animation as gif\noutput_video_filename = '/kaggle/working/dicom_scan_seg.gif'\ndicom_scan_seg.save(output_video_filename, writer='pillow', fps=10)  # You can adjust the 'fps' parameter as needed","metadata":{"execution":{"iopub.status.busy":"2023-09-20T15:42:12.891835Z","iopub.execute_input":"2023-09-20T15:42:12.892207Z","iopub.status.idle":"2023-09-20T15:42:59.348608Z","shell.execute_reply.started":"2023-09-20T15:42:12.892174Z","shell.execute_reply":"2023-09-20T15:42:59.347206Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n---\n# Other files\n\n---\n## Tags","metadata":{}},{"cell_type":"code","source":"train_tags = pd.read_parquet(f'{BASE_PATH}/train_dicom_tags.parquet')","metadata":{"execution":{"iopub.status.busy":"2023-09-21T16:08:49.509206Z","iopub.execute_input":"2023-09-21T16:08:49.509952Z","iopub.status.idle":"2023-09-21T16:08:56.168519Z","shell.execute_reply.started":"2023-09-21T16:08:49.509899Z","shell.execute_reply":"2023-09-21T16:08:56.167471Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_tags.columns","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-09-21T16:10:36.216947Z","iopub.execute_input":"2023-09-21T16:10:36.217453Z","iopub.status.idle":"2023-09-21T16:10:36.226222Z","shell.execute_reply.started":"2023-09-21T16:10:36.217413Z","shell.execute_reply":"2023-09-21T16:10:36.225048Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"One thing that will be useful soon is the Slice Thickness. Ideally, we would have all our slices evenly thick.","metadata":{}},{"cell_type":"code","source":"train_tags.SliceThickness.value_counts().sort_index()","metadata":{"execution":{"iopub.status.busy":"2023-09-21T16:46:48.332508Z","iopub.execute_input":"2023-09-21T16:46:48.333437Z","iopub.status.idle":"2023-09-21T16:46:48.378544Z","shell.execute_reply.started":"2023-09-21T16:46:48.333377Z","shell.execute_reply":"2023-09-21T16:46:48.376622Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There are duplicates certainly due to the precision, a 3 decimal places precision should be enough","metadata":{}},{"cell_type":"code","source":"thickness_data = train_tags.SliceThickness.round(3).value_counts().sort_index()\nthickness_df = pd.DataFrame({'Slice Thickness': thickness_data.index, 'Count': thickness_data.values})\n\n# Compute the cumulative sum of ordered counts\nthickness_df.sort_values(by='Count', ascending=False, inplace=True)\nthickness_df['Cumulative Sum'] = thickness_df['Count'].cumsum()\n\n# Compute the percentage of the overall total for each step\ntotal_count = thickness_df['Count'].sum()\nthickness_df['Percentage of Total'] = (thickness_df['Cumulative Sum'] / total_count) * 100\n\n# Create a bar chart\nplt.figure(figsize=(15, 3))  # Set the figure size\nplt.bar(thickness_df['Slice Thickness'], thickness_df['Count'], width=0.1)\n\n# Customize the chart\nplt.xlabel('Slice Thickness')\nplt.ylabel('Count')\nplt.title('Slice Thickness Count')\n\n# Show the chart\nplt.xticks(thickness_df['Slice Thickness'], rotation=45)  # Rotate x-axis labels for better visibility\nplt.tight_layout()  # Adjust layout for better display\nplt.show()\n\n# Display the DataFrame with cumulative sum and percentages\nprint(thickness_df[['Slice Thickness', 'Cumulative Sum', 'Percentage of Total']])","metadata":{"execution":{"iopub.status.busy":"2023-09-21T16:40:29.311435Z","iopub.execute_input":"2023-09-21T16:40:29.311906Z","iopub.status.idle":"2023-09-21T16:40:29.817235Z","shell.execute_reply.started":"2023-09-21T16:40:29.311873Z","shell.execute_reply":"2023-09-21T16:40:29.815402Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n---\n# Patient cases\n\nAfter discovering how to load and visualize the data, I am interested in looking at specific cases where patient have such and such disease, to have an idea of how to perceive them from the images.\n\nWe go through each \"category\" (bowel, extravasation, kidney, liver, spleen) and compare a selection of healthy ones vs injured ones.\nEach time, I add medical information found online as I have no special training in that domain. \n\n---\n## Bowels\n\nFirst we have to find a patient with an injured bowel and one with an healthy one.","metadata":{}},{"cell_type":"code","source":"data_patient_injury_bowel = data[data[\"bowel_injury\"]==1].head(1)\npatient_injury_bowel = data_patient_injury_bowel[\"patient_id\"].values[0]\ndata_patient_injury_bowel","metadata":{"execution":{"iopub.status.busy":"2023-09-21T15:24:32.659219Z","iopub.execute_input":"2023-09-21T15:24:32.659848Z","iopub.status.idle":"2023-09-21T15:24:32.685834Z","shell.execute_reply.started":"2023-09-21T15:24:32.65979Z","shell.execute_reply":"2023-09-21T15:24:32.684176Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We might as well try to get a completely healthy patient.","metadata":{}},{"cell_type":"code","source":"data_patient_healthy_bowel = data[data[\"any_injury\"]==0].head(1)\npatient_healthy_bowel = data_patient_healthy_bowel[\"patient_id\"].values[0]\ndata_patient_healthy_bowel","metadata":{"execution":{"iopub.status.busy":"2023-09-21T15:24:34.905951Z","iopub.execute_input":"2023-09-21T15:24:34.906705Z","iopub.status.idle":"2023-09-21T15:24:34.925704Z","shell.execute_reply.started":"2023-09-21T15:24:34.906664Z","shell.execute_reply":"2023-09-21T15:24:34.924344Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So there actually are completely healthy patients in this dataset.\nIn order to compare them we must load their data, display them side by side and align them too.\nThis should be helpful later on if we use a model that does ","metadata":{}},{"cell_type":"code","source":"# load data for both patients\n\n# How many series for each patient and how many files in each serie ?\npatients_bowel = [patient_injury_bowel, patient_healthy_bowel]\ndf_patients_bowel = pd.DataFrame(columns=[\"patient_id\", \"serie_id\", \"image_id\"])\nfor patient_id in patients_bowel:\n    series = os.listdir(f\"{BASE_PATH}/train_images/{patient_id}\")\n    for serie_id in series:\n        images = os.listdir(f\"{BASE_PATH}/train_images/{patient_id}/{serie_id}\")\n        rows = {\"patient_id\" : patient_id, \"serie_id\" : int(serie_id), \"image_id\" : images}\n        df_patients_bowel = pd.concat([df_patients_bowel, pd.DataFrame.from_dict(rows)])","metadata":{"execution":{"iopub.status.busy":"2023-09-21T15:24:37.282568Z","iopub.execute_input":"2023-09-21T15:24:37.282995Z","iopub.status.idle":"2023-09-21T15:24:37.425624Z","shell.execute_reply.started":"2023-09-21T15:24:37.282958Z","shell.execute_reply":"2023-09-21T15:24:37.424445Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_patients_bowel.groupby([\"patient_id\", \"serie_id\"]).count()","metadata":{"execution":{"iopub.status.busy":"2023-09-21T15:24:38.854109Z","iopub.execute_input":"2023-09-21T15:24:38.854821Z","iopub.status.idle":"2023-09-21T15:24:38.876632Z","shell.execute_reply.started":"2023-09-21T15:24:38.854783Z","shell.execute_reply":"2023-09-21T15:24:38.875307Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pydicom.read_file(f\"{BASE_PATH}/train_images/10005/18667/100.dcm\").SliceThickness, \\\npydicom.read_file(f\"{BASE_PATH}/train_images/10065/37324/100.dcm\").SliceThickness, \\\npydicom.read_file(f\"{BASE_PATH}/train_images/10065/46839/100.dcm\").SliceThickness,","metadata":{"execution":{"iopub.status.busy":"2023-09-21T15:55:48.547878Z","iopub.execute_input":"2023-09-21T15:55:48.54839Z","iopub.status.idle":"2023-09-21T15:55:48.567261Z","shell.execute_reply.started":"2023-09-21T15:55:48.548325Z","shell.execute_reply":"2023-09-21T15:55:48.565602Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Open the files \n# would be great to have several sliders side by side\n# even views from different angles\n\ndef open_dicom_image(image_path, dsize=(256, 256)):\n    dcm_file = dicomsdl.open(image_path)\n    \n    info = dcm_file.getPixelDataInfo() # rescaling it taken care of in dicomsdl\n    \n    # find clipping params\n    center = int(info[\"WindowCenter\"])\n    width = int(info[\"WindowWidth\"])\n    low = center - width / 2\n    high = center + width / 2\n    \n    # rescale and clip image\n    image = dcm_file.pixelData()\n    image = np.clip(image, low, high)\n    \n    image = (image - image.min()) / (image.max() - image.min())\n\n    if info['PhotometricInterpretation'] == \"MONOCHROME1\":\n        image = 1 - image\n\n    data = cv2.resize(image, dsize=dsize)\n    data = (data * 255).astype(np.uint8)\n    return data\n\n    \ndef plot_dicom(df_patients_to_compare):\n    patients_series = df_patients_bowel[[\"patient_id\", \"serie_id\"]].drop_duplicates(ignore_index=True)\n    patients_series_image_paths = []\n    for index, row in patients_series.iterrows():\n        patient_id = row[\"patient_id\"]\n        serie_id = row[\"serie_id\"]\n        image_files = df_patients_bowel.loc[(df_patients_bowel[\"patient_id\"]==patient_id) & (df_patients_bowel[\"serie_id\"]==serie_id), \"image_id\"]\n        image_files.sort_values(key = lambda x : x.str.slice(stop=-4).astype(int), inplace=True, ignore_index=True)\n        image_paths = f'{BASE_PATH}/train_images/{patient_id}/{serie_id}/' + image_files\n        patients_series_image_paths.append(image_paths)\n    fig, axs = plt.subplots(ncols=len(patients_series_image_paths))\n    ims = []\n    for j in range(len(patients_series_image_paths)):\n        im = axs[j].imshow(open_dicom_image(patients_series_image_paths[j][0]), cmap=plt.cm.bone)\n        ims.append(im)\n    \n    def update(i):\n        for j in range(len(patients_series_image_paths)):\n            ims[j].set_array(open_dicom_image(patients_series_image_paths[j][i]))\n\n    ani = animation.FuncAnimation(fig, update, frames=range(100), repeat=True)\n    return ani","metadata":{"execution":{"iopub.status.busy":"2023-09-21T15:37:24.648541Z","iopub.execute_input":"2023-09-21T15:37:24.64906Z","iopub.status.idle":"2023-09-21T15:37:24.669251Z","shell.execute_reply.started":"2023-09-21T15:37:24.649023Z","shell.execute_reply":"2023-09-21T15:37:24.667745Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"HTML(plot_dicom(df_patients_bowel).to_jshtml())","metadata":{"execution":{"iopub.status.busy":"2023-09-21T15:37:25.639788Z","iopub.execute_input":"2023-09-21T15:37:25.640272Z","iopub.status.idle":"2023-09-21T15:37:44.983229Z","shell.execute_reply.started":"2023-09-21T15:37:25.640235Z","shell.execute_reply":"2023-09-21T15:37:44.982046Z"},"trusted":true},"execution_count":null,"outputs":[]}]}