{"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":"# MAGNETIC RESONANCE IMAGING (MRI)\n\n**N.B.** This excerpt is taken from https://case.edu/med/neurology/NR/MRI%20Basics.htm \n\nMagnetic resonance imaging (MRI) is one of the most commonly used tests in neurology and neurosurgery. MRI provides exquisite detail of brain, spinal cord and vascular anatomy, and has the advantage of being able to visualize anatomy in all three planes: \n#### 1. axial\n#### 2. sagittal &\n#### 3. coronal \n\n(see the example image below). \n\n![](https://case.edu/med/neurology/NR/mri%20slices%20new.jpg \"axial, sagittal and coronal plane\")\n\n# MRI IMAGING SEQUENCES\nThere are different types of MRI sequences. The most commons are the \n#### 1. T1-weighted (T1w)\n#### 2. T2-weighted (T2w)\n#### 3. Fluid Attenuated Inversion Recovery (FLAIR)\n\nThese sequences vary from each other in their TE (Time to Echo) and TR (Repetition Time) e.g. T1-weighted images are produced by using short TE and TR times while T2-weighted images are produced by using longer TE and TR times. \nT1-weighted imaging can also be performed while infusing Gadolinium (Gad). And that is the 4th type of sequence in this dataset:\n#### 4. T1-weighted Gadolinium Post Contrast (T1wCE/T1Gd)","metadata":{}},{"cell_type":"markdown","source":"# Comparison of Different Sequences\n\nN.B. image taken from this paper https://link.springer.com/article/10.1007/s10278-019-00282-4\n\n![](https://media.springernature.com/lw685/springer-static/image/art%3A10.1007%2Fs10278-019-00282-4/MediaObjects/10278_2019_282_Fig1_HTML.png)\n\n\n**In a T1 weighted sequence:**\n* fluid (CSF): black/dark - low signal intensity\n* fat: white/light - high signal intensity\n* grey matter: grey\n* white matter: white-ish\n\n**In a T2 weighted sequence:**\n* fluid (CSF): white/bright\n* fat: white/light\n* grey matter: grey\n* white matter: black/dark\n\nFLAIR images appear similar to T1 (CSF is dark). The best way to tell the two apart is to look at the grey-white matter. T1 sequences will have grey matter being darker than white matter.\n\nT1wGd images has bright signals in the blood vessels due to the Gadolinium.","metadata":{}},{"cell_type":"code","source":"import glob\nimport pandas as pd\npd.set_option('display.max_colwidth', -1)  \nimport numpy as np\nfrom tqdm import tqdm\n\nimport matplotlib.pyplot as plt\nimport matplotlib\nmatplotlib.rcParams['animation.html'] = 'jshtml'\nimport seaborn as sns\n\nimport tensorflow as tf\nimport tensorflow.keras as keras\n\nimport pydicom","metadata":{"execution":{"iopub.status.busy":"2021-08-09T17:06:04.679236Z","iopub.execute_input":"2021-08-09T17:06:04.679615Z","iopub.status.idle":"2021-08-09T17:06:04.686052Z","shell.execute_reply.started":"2021-08-09T17:06:04.679584Z","shell.execute_reply":"2021-08-09T17:06:04.685168Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"root_dir = '../input/rsna-miccai-brain-tumor-radiogenomic-classification/'\n\ndf = pd.read_csv(root_dir+'train_labels.csv')\n\nsns.countplot(data=df, x='MGMT_value')","metadata":{"execution":{"iopub.status.busy":"2021-08-09T17:06:04.687604Z","iopub.execute_input":"2021-08-09T17:06:04.68804Z","iopub.status.idle":"2021-08-09T17:06:04.827771Z","shell.execute_reply.started":"2021-08-09T17:06:04.688009Z","shell.execute_reply":"2021-08-09T17:06:04.827014Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# What is the 'MGMT_value'?\n\nO[6]-methylguanine-DNA methyltransferase (MGMT) is a protein in cells, including tumour cells, that repairs damage to the cell’s DNA. For example, the damage caused by chemotherapy drugs to tumour cells. The more MGMT protein that the tumour produces, the less effective the chemotherapy drug is expected to be, as the protein will repair the damage to the tumour. Thus, determination of MGMT promoter methylation status in newly diagnosed GBM can influence treatment decision making.\n\nIn this dataset the MGMT promoter methylation status data is defined as a binary label (0: unmethylated, 1: methylated)\n\n","metadata":{}},{"cell_type":"code","source":"# Add the full paths for each id for different types of sequences to the csv \ndef full_ids(data):\n    zeros = 5 - len(str(data))\n    if zeros > 0:\n        prefix = ''.join(['0' for i in range(zeros)])\n    \n    return prefix+str(data)\n        \n\ndf['BraTS21ID_full'] = df['BraTS21ID'].apply(full_ids)\n\n# Add all the paths to the df for easy access\ndf['flair'] = df['BraTS21ID_full'].apply(lambda file_id : root_dir+'train/'+file_id+'/FLAIR/')\ndf['t1w'] = df['BraTS21ID_full'].apply(lambda file_id : root_dir+'train/'+file_id+'/T1w/')\ndf['t1wce'] = df['BraTS21ID_full'].apply(lambda file_id : root_dir+'train/'+file_id+'/T1wCE/')\ndf['t2w'] = df['BraTS21ID_full'].apply(lambda file_id : root_dir+'train/'+file_id+'/T2w/')\ndf","metadata":{"execution":{"iopub.status.busy":"2021-08-09T17:06:04.829087Z","iopub.execute_input":"2021-08-09T17:06:04.829468Z","iopub.status.idle":"2021-08-09T17:06:04.86274Z","shell.execute_reply.started":"2021-08-09T17:06:04.829439Z","shell.execute_reply":"2021-08-09T17:06:04.861604Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# DICOM\nDICOM is the international standard to communicate and manage medical images and data. Its mission is to ensure the interoperability of systems used to produce, store, share, display, send, query, process, retrieve and print medical images, as well as to manage related workflows. \n\nLet's see how to read a dicom file with pydicom and what it looks like.\n\nDICOM holds all MRI data such as: sequence type, MR acquisition type, patient orientation etc. along with the main MRI image.\n","metadata":{}},{"cell_type":"code","source":"data = pydicom.dcmread('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00000/T1w/Image-24.dcm')\ndata","metadata":{"execution":{"iopub.status.busy":"2021-08-09T17:06:04.880153Z","iopub.execute_input":"2021-08-09T17:06:04.880651Z","iopub.status.idle":"2021-08-09T17:06:04.925612Z","shell.execute_reply.started":"2021-08-09T17:06:04.880546Z","shell.execute_reply":"2021-08-09T17:06:04.924663Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_image(data):\n    '''\n    Returns the image data as a numpy array.\n    '''  \n    if np.max(data.pixel_array)==0:\n        img = data.pixel_array\n    else:\n        img = data.pixel_array/np.max(data.pixel_array)\n        img = (img * 255).astype(np.uint8)\n        \n    return img\n\ndata = pydicom.dcmread('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00000/T1w/Image-24.dcm')\nimg = get_image(data)\nplt.imshow(img, cmap='gray')","metadata":{"execution":{"iopub.status.busy":"2021-08-09T17:06:04.927223Z","iopub.execute_input":"2021-08-09T17:06:04.927654Z","iopub.status.idle":"2021-08-09T17:06:05.140166Z","shell.execute_reply.started":"2021-08-09T17:06:04.927605Z","shell.execute_reply":"2021-08-09T17:06:05.139109Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Display Sequence of MRIs as Animation\nEach subfolder within individual case corresponds to a MRI Imaging Sequence. I am gonna use Matplotlib here to display the sequence of images as animation. This animations can also be saved in video format.","metadata":{}},{"cell_type":"code","source":"def sorted_image_dirs(path: str):\n    '''\n    Sorts the list of image directories by image number in a path\n    '''\n    dirs = glob.glob(path+'*')\n    dirs.sort(key=lambda x: int(x.split('/')[-1].split('-')[-1].split('.')[0]))\n    \n    return dirs\n\n\ndef get_all_images(path: str):\n    '''\n    Returns a list of (non blank) images from a given path (of shape [non_blank_image_count, 512, 512])\n    '''\n    image_dirs = sorted_image_dirs(path)\n    images = []\n    \n    for directory in image_dirs:\n        data = pydicom.dcmread(directory)\n        img = get_image(data)\n        \n        # Exclude the blank images\n        if np.max(img)!=0:\n            images.append(img)\n        else:\n            pass\n    \n    return images\n    \ndef show_animation(images: list):\n    '''\n    Displays an animation from the list of images.\n    \n    set: matplotlib.rcParams['animation.html'] = 'jshtml'\n    \n    '''\n    fig = plt.figure(figsize=(6, 6))\n    plt.axis('off')\n    im = plt.imshow(images[0], cmap='gray')\n    \n    def animate_func(i):\n        im.set_array(images[i])\n        return [im]\n    \n    return matplotlib.animation.FuncAnimation(fig, animate_func, frames = len(images), interval = 20)","metadata":{"execution":{"iopub.status.busy":"2021-08-09T17:06:05.142232Z","iopub.execute_input":"2021-08-09T17:06:05.142526Z","iopub.status.idle":"2021-08-09T17:06:05.152926Z","shell.execute_reply.started":"2021-08-09T17:06:05.142496Z","shell.execute_reply":"2021-08-09T17:06:05.15186Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Show FLAIR Sequence of a Patient as Images","metadata":{}},{"cell_type":"code","source":"# A random patient - patient 10\npatient = 10\n\nflair_images = get_all_images(df['flair'][patient])\nprint('No of images:', len(flair_images))\nprint('MGMT: ', df['MGMT_value'][patient])\n\nfig = plt.figure(figsize=(30,30))\n\nc = 1\nfor image in flair_images:\n    ax = fig.add_subplot(len(flair_images)//10+1, 10, c)\n    ax.imshow(image, cmap='gray')\n    c+=1\n    \n    plt.axis('off')\n    \nfig.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2021-08-09T17:06:05.155092Z","iopub.execute_input":"2021-08-09T17:06:05.155444Z","iopub.status.idle":"2021-08-09T17:06:15.134499Z","shell.execute_reply.started":"2021-08-09T17:06:05.155415Z","shell.execute_reply":"2021-08-09T17:06:15.133466Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Show Flair Sequence as Animation","metadata":{}},{"cell_type":"code","source":"print('No of images:', len(flair_images))\n\nflair_animation = show_animation(flair_images)\nflair_animation\n# flair_animation.save('./a.mp4')","metadata":{"execution":{"iopub.status.busy":"2021-08-09T17:06:15.135774Z","iopub.execute_input":"2021-08-09T17:06:15.136045Z","iopub.status.idle":"2021-08-09T17:06:19.527976Z","shell.execute_reply.started":"2021-08-09T17:06:15.136019Z","shell.execute_reply":"2021-08-09T17:06:19.526892Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Show T1w Sequence as Animation","metadata":{}},{"cell_type":"code","source":"t1w_images = get_all_images(df['t1w'][patient])\n    \nprint('No of images:', len(t1w_images))\nshow_animation(t1w_images)","metadata":{"execution":{"iopub.status.busy":"2021-08-09T17:06:19.529251Z","iopub.execute_input":"2021-08-09T17:06:19.52954Z","iopub.status.idle":"2021-08-09T17:06:21.838737Z","shell.execute_reply.started":"2021-08-09T17:06:19.529511Z","shell.execute_reply":"2021-08-09T17:06:21.837832Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# MRI Plane\nAfter reading a dicom file we can get all the information we need like patient orientation, pixel data (Each of the image is a 512x512 array) etc.\nAs mentioned at the begining of this notebook, MRI can be visualized in 3 different planes - Axial(Transverse), Coronal, Sagittal. The plane of an MR image relative to the patient's body can be calculated using the DICOM Image Orientation Patient tag. \n\n![](https://my-ms.org/images/mri_planes_gnu.jpg)","metadata":{}},{"cell_type":"code","source":"# https://www.kaggle.com/davidbroberts/determining-mr-image-planes\ndef get_image_plane(data):\n    '''\n    Returns the MRI's plane from the dicom data.\n    \n    '''\n    x1,y1,_,x2,y2,_ = [round(j) for j in data.ImageOrientationPatient]\n    cords = [x1,y1,x2,y2]\n\n    if cords == [1,0,0,0]:\n        return 'coronal'\n    if cords == [1,0,0,1]:\n        return 'axial'\n    if cords == [0,1,0,0]:\n        return 'sagittal'","metadata":{"execution":{"iopub.status.busy":"2021-08-09T17:06:21.840159Z","iopub.execute_input":"2021-08-09T17:06:21.840567Z","iopub.status.idle":"2021-08-09T17:06:21.847769Z","shell.execute_reply.started":"2021-08-09T17:06:21.840523Z","shell.execute_reply":"2021-08-09T17:06:21.846539Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Plot of Some Random MRIs with Their Respective Sequences, Planes & MGMT_values\nAs each of the sequence contains multiple images, I am gonna display the middle image of the sequence.","metadata":{}},{"cell_type":"code","source":"def get_middle_image(path: str, plane=False):\n    '''\n    Returns the middle image from the path\n    \n    if plane=True returns the plane\n    '''\n    image_dirs = sorted_image_dirs(path)\n\n    dicom = pydicom.dcmread(image_dirs[len(image_dirs)//2])\n    img = get_image(dicom)\n    plane = get_image_plane(dicom)\n    \n    if plane:\n        return img, plane\n    else:\n        return img","metadata":{"execution":{"iopub.status.busy":"2021-08-09T17:06:21.849294Z","iopub.execute_input":"2021-08-09T17:06:21.849759Z","iopub.status.idle":"2021-08-09T17:06:21.866677Z","shell.execute_reply.started":"2021-08-09T17:06:21.849716Z","shell.execute_reply":"2021-08-09T17:06:21.865389Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(35,20))\n\nseq_types = ['flair', 't1w', 't1wce', 't2w']\n\nfor i in range(20):\n    \n    patient = np.random.randint(low=0, high=len(df))\n    seq_type = np.random.choice(seq_types)\n\n    # path for the randomly selected image and sequence type\n    seq_path = df[seq_type][patient]\n\n    # Get the middle image and plane from the path\n    img, plane = get_middle_image(seq_path, plane=True)\n    \n    patient_id, mgmt = df['BraTS21ID_full'][patient], df['MGMT_value'][patient]\n    \n    ax = fig.add_subplot(4,5,i+1)\n    ax.imshow(img, cmap='gray')\n    plt.title(f'ID: {patient_id}, MGMT_value: {mgmt}, Plane: {plane}, Seq_type: {seq_type}')\n    plt.axis('off')","metadata":{"execution":{"iopub.status.busy":"2021-08-09T17:06:21.868088Z","iopub.execute_input":"2021-08-09T17:06:21.868426Z","iopub.status.idle":"2021-08-09T17:06:25.289138Z","shell.execute_reply.started":"2021-08-09T17:06:21.868395Z","shell.execute_reply":"2021-08-09T17:06:25.2884Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# The Planes of 4 type of Sequences for the Same Patient are Not the Same!\n\nAfter plotting one image (the middle image) from each sequence for the first patient (ID 00000) we can see that, not all the axis are the same. Different sequences planes for one patient can be different.","metadata":{}},{"cell_type":"code","source":"fig = plt.figure(figsize=(35,20))\n\nseq_types = ['flair', 't1w', 't1wce', 't2w']\n\npatient = 0  # The first patient\n\nfor i in range(4):\n    \n    seq_type = seq_types[i]\n    seq_path = df[seq_type][patient]\n\n    # Get the middle image and plane from the path\n    img, plane = get_middle_image(seq_path, plane=True)\n    \n    patient_id, mgmt= df['BraTS21ID_full'][patient], df['MGMT_value'][patient]\n    \n    ax = fig.add_subplot(1,4,i+1)\n    ax.imshow(img, cmap='gray')\n    plt.title(f'ID: {patient_id}, MGMT_value: {mgmt}, Plane: {plane}, Seq_type: {seq_type}')\n    plt.axis('off')","metadata":{"execution":{"iopub.status.busy":"2021-08-09T17:06:25.290147Z","iopub.execute_input":"2021-08-09T17:06:25.290567Z","iopub.status.idle":"2021-08-09T17:06:26.187607Z","shell.execute_reply.started":"2021-08-09T17:06:25.290525Z","shell.execute_reply":"2021-08-09T17:06:26.186874Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Add the plane data for each sequence to the exisitng dataframe ","metadata":{}},{"cell_type":"code","source":"seq_types = ['flair', 't1w', 't1wce', 't2w']\n\nseq_axes = {'flair':[], 't1w':[], 't1wce':[], 't2w':[]}\n\nfor index in tqdm(range(len(df))):\n    \n    for seq in seq_types:\n        \n        seq_path = df[seq][index]\n\n        # Get the middle image from the path\n        img, plane = get_middle_image(seq_path, plane=True)\n    \n        seq_axes[seq].append(plane)\n        \n# Add the axes to the dataframe\ndf['flair_axis'] = seq_axes['flair']\ndf['t1w_axis'] = seq_axes['t1w']\ndf['t1wce_axis'] = seq_axes['t1wce']\ndf['t2w_axis'] = seq_axes['t2w']\n\n# Print the axes with the respective patient ID\ndf[['BraTS21ID_full', 'flair_axis', 't1w_axis', 't1wce_axis', 't2w_axis']]","metadata":{"execution":{"iopub.status.busy":"2021-08-09T17:06:26.188693Z","iopub.execute_input":"2021-08-09T17:06:26.189233Z","iopub.status.idle":"2021-08-09T17:08:49.265135Z","shell.execute_reply.started":"2021-08-09T17:06:26.189172Z","shell.execute_reply":"2021-08-09T17:08:49.264014Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots(1,4, figsize=(30,10))\n\nsns.countplot(x='flair_axis', data=df, ax=ax[0])\nsns.countplot(x='t1w_axis', data=df, ax=ax[1])\nsns.countplot(x='t1wce_axis', data=df, ax=ax[2])\nsns.countplot(x='t2w_axis', data=df, ax=ax[3])\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-09T17:08:49.266748Z","iopub.execute_input":"2021-08-09T17:08:49.267157Z","iopub.status.idle":"2021-08-09T17:08:49.71787Z","shell.execute_reply.started":"2021-08-09T17:08:49.267115Z","shell.execute_reply":"2021-08-09T17:08:49.717167Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Important: \nIn this [discussion](https://www.kaggle.com/c/rsna-miccai-brain-tumor-radiogenomic-classification/discussion/262046) a competition host has notified that there are some issues with these 3 cases - \n\nPatient IDs -  \n* 00109 (FLAIR images are blank) \n* 00123 (T1w images are blank) \n* 00709 (FLAIR images are blank)\n\nHence these can be excluded.","metadata":{}},{"cell_type":"code","source":"to_exclude = [109, 123, 709]\n\ndf = df[~df['BraTS21ID'].isin(to_exclude)]","metadata":{"execution":{"iopub.status.busy":"2021-08-09T17:14:10.934544Z","iopub.execute_input":"2021-08-09T17:14:10.934935Z","iopub.status.idle":"2021-08-09T17:14:10.940748Z","shell.execute_reply.started":"2021-08-09T17:14:10.934904Z","shell.execute_reply":"2021-08-09T17:14:10.939876Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Work In Progress","metadata":{}}]}