{"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":"![brain_baner](http://www.mf-data-science.fr/images/projects/brain_baner.jpg)","metadata":{"papermill":{"duration":0.020401,"end_time":"2021-08-27T14:49:47.254332","exception":false,"start_time":"2021-08-27T14:49:47.233931","status":"completed"},"tags":[]}},{"cell_type":"markdown","source":"

Context

\n\nThe goal of this competition, initiated by the **Radiological Society of North America *(RSNA)*** in partnership with the **Medical Image Computing and Computer Assisted Intervention Society *(the MICCAI Society)*** is to predict the methylation of the **MGMT promoter**, which is an important gene biomarker for treatment of brain tumors.\n\nThese predictions will be based on a database of **MRI *(magnetic resonance imaging)*** scans of several hundred patients.\n\n

Data

\n\nEach independent case has a dedicated folder identified by a five-digit number. Within each of these “case” folders, there are four sub-folders, each of them corresponding to each of the structural multi-parametric MRI (mpMRI) scans, in DICOM format. The exact mpMRI scans included are:\n\n- Fluid Attenuated Inversion Recovery (FLAIR)\n- T1-weighted pre-contrast (T1w)\n- T1-weighted post-contrast (T1Gd)\n- T2-weighted (T2)\n\n| ![brain_baner](http://www.mf-data-science.fr/images/projects/brain_tumor_types.png) | \n|:--:| \n| *Examples of the four MR sequence types included in this work* |","metadata":{}},{"cell_type":"markdown","source":"

Acknowledgement

\n\nThis Notebook is copied and edited from *Ammar Alhaj Ali* work :\n- [🧠Brain Tumor 3D [Training]](https://www.kaggle.com/ammarnassanalhajali/brain-tumor-3d-training)\n- [🧠Brain Tumor 3D [Inference]](https://www.kaggle.com/ammarnassanalhajali/brain-tumor-3d-inference)\n\nIn the first part, we will add to these Notebooks **a complete preprocessing of the images** (crop, equalization, denoising ...) in order to try to improve the results and **reduce the size of images**.","metadata":{}},{"cell_type":"markdown","source":"# Exploratory data analysis (EDA)\n\nFirst, we have to load the usefull Python libraries :","metadata":{"papermill":{"duration":0.018755,"end_time":"2021-08-27T14:49:47.330074","exception":false,"start_time":"2021-08-27T14:49:47.311319","status":"completed"},"tags":[]}},{"cell_type":"code","source":"import os\nimport sys \nimport json\nimport glob\nimport random\nimport collections\nimport time\nimport re\nimport math\nimport numpy as np\nimport pandas as pd\nimport cv2\n\nimport matplotlib.pyplot as plt\nfrom matplotlib import animation, rc\nimport seaborn as sns\nimport pydicom\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\n\nfrom random import shuffle\nfrom sklearn import model_selection as sk_model_selection\n\nfrom tensorflow import keras\nfrom tensorflow.keras import layers\nfrom tensorflow.keras.callbacks import ModelCheckpoint, EarlyStopping, ReduceLROnPlateau\nfrom keras.metrics import AUC\n\nimport tensorflow as tf","metadata":{"papermill":{"duration":5.717903,"end_time":"2021-08-27T14:49:53.066958","exception":false,"start_time":"2021-08-27T14:49:47.349055","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-09-01T08:27:56.468642Z","iopub.execute_input":"2021-09-01T08:27:56.469062Z","iopub.status.idle":"2021-09-01T08:28:03.613265Z","shell.execute_reply.started":"2021-09-01T08:27:56.468975Z","shell.execute_reply":"2021-09-01T08:28:03.611954Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Global params for animations\nrc('animation', html='jshtml')","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:28:03.619331Z","iopub.execute_input":"2021-09-01T08:28:03.622056Z","iopub.status.idle":"2021-09-01T08:28:03.629720Z","shell.execute_reply.started":"2021-09-01T08:28:03.622010Z","shell.execute_reply":"2021-09-01T08:28:03.628204Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"and configure data path :","metadata":{}},{"cell_type":"code","source":"data_directory = '../input/rsna-miccai-brain-tumor-radiogenomic-classification'\npytorch3dpath = \"../input/efficientnetpyttorch3d/EfficientNet-PyTorch-3D\"","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:28:03.636571Z","iopub.execute_input":"2021-09-01T08:28:03.640002Z","iopub.status.idle":"2021-09-01T08:28:03.647683Z","shell.execute_reply.started":"2021-09-01T08:28:03.639955Z","shell.execute_reply":"2021-09-01T08:28:03.646177Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Loading data\n\nLet's load the **dataframe with training data and testing data**. We are adding a column that will contain the patient record ID of 5 characters.","metadata":{"papermill":{"duration":0.020721,"end_time":"2021-08-27T14:49:53.107826","exception":false,"start_time":"2021-08-27T14:49:53.087105","status":"completed"},"tags":[]}},{"cell_type":"code","source":"train_df = pd.read_csv(data_directory+\"/train_labels.csv\")\ntrain_df['BraTS21ID5'] = [format(x, '05d') for x in train_df.BraTS21ID]\ntrain_df.head(3)","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:28:03.653368Z","iopub.execute_input":"2021-09-01T08:28:03.655578Z","iopub.status.idle":"2021-09-01T08:28:03.711779Z","shell.execute_reply.started":"2021-09-01T08:28:03.655531Z","shell.execute_reply":"2021-09-01T08:28:03.710766Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test = pd.read_csv(\n data_directory+'/sample_submission.csv')\n\ntest['BraTS21ID5'] = [format(x, '05d') for x in test.BraTS21ID]\ntest.head(3)","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:28:03.716542Z","iopub.execute_input":"2021-09-01T08:28:03.719037Z","iopub.status.idle":"2021-09-01T08:28:03.747682Z","shell.execute_reply.started":"2021-09-01T08:28:03.718992Z","shell.execute_reply":"2021-09-01T08:28:03.746684Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"All of the subjects in this dataset appear to have a brain tumor. MGMT_Class = 0 refers to people who do not have the MGMT promoter methylation. MGMT_Class = 1 appears to be someone who has the MGMT promoter methylation. \n\nIt is a competition which gives **the probability of it in `MGMT_value`** feature. `BraTS21ID` is the patient's identification. We can see that the data structure is the same as for the submission file than the train file, knowing that here **MGMT_value is indeed equal to 0 or 1 and no longer a probability**.\n\nLet's look at the distribution of the values of this variable in the train set :","metadata":{}},{"cell_type":"code","source":"sns.set_style(\"whitegrid\")\nfig = plt.figure(figsize=(8,6))\n# Countplot with Seaborn\nax = sns.countplot(data=train_df,\n x=\"MGMT_value\")\n# Annotating bars\nfor p in ax.patches:\n ax.annotate(\n format(p.get_height(), '.0f'), \n (p.get_x() + p.get_width() / 2., p.get_height()),\n ha = 'center', va = 'center', \n xytext = (0, 10), \n textcoords = 'offset points')\n\nsns.despine(left=True, bottom=True)\nplt.title(\"MGMT value distribution in train labels\\n\",\n fontsize=18, color=\"#0b0a2d\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:28:03.752347Z","iopub.execute_input":"2021-09-01T08:28:03.755035Z","iopub.status.idle":"2021-09-01T08:28:04.044418Z","shell.execute_reply.started":"2021-09-01T08:28:03.754976Z","shell.execute_reply":"2021-09-01T08:28:04.043333Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## MRI data\n\nTrain data contains one record per patient. For each patient, four sub-files are available *(FLAIR, T1w, T1wCE and T2w)* in which the MRI image sequences are distributed.\n\n![data_structure](http://www.mf-data-science.fr/images/projects/data_structure.jpg)\n\nWe are going to take a look at what an MRI image looks like :","metadata":{}},{"cell_type":"code","source":"def load_dicom(path):\n dicom = pydicom.read_file(path)\n data = dicom.pixel_array\n data = data - np.min(data)\n if np.max(data) != 0:\n data = data / np.max(data)\n data = (data * 255).astype(np.uint8)\n return data\n\ndef visualize_sample(\n brats21id, \n slice_i,\n mgmt_value,\n types=(\"FLAIR\", \"T1w\", \"T1wCE\", \"T2w\")):\n \n plt.figure(figsize=(20, 6))\n patient_path = os.path.join(\n data_directory+\"/train/\", \n str(brats21id).zfill(5),\n )\n for i, t in enumerate(types, 1):\n t_paths = sorted(\n glob.glob(os.path.join(patient_path, t, \"*\")), \n key=lambda x: int(x[:-4].split(\"-\")[-1]),\n )\n data = load_dicom(t_paths[int(len(t_paths) * slice_i)])\n plt.subplot(1, 4, i)\n plt.imshow(data, cmap=\"gray\")\n plt.title(f\"{t}\", fontsize=10)\n plt.axis(\"off\")\n plt.show()","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:28:04.049359Z","iopub.execute_input":"2021-09-01T08:28:04.052168Z","iopub.status.idle":"2021-09-01T08:28:04.273492Z","shell.execute_reply.started":"2021-09-01T08:28:04.052121Z","shell.execute_reply":"2021-09-01T08:28:04.272228Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### With MGMT = 0","metadata":{}},{"cell_type":"code","source":"list0=[315,176,153,164]\nfor i in list0:\n _brats21id = train_df.iloc[i][\"BraTS21ID\"]\n _mgmt_value = train_df.iloc[i][\"MGMT_value\"]\n visualize_sample(brats21id=_brats21id, mgmt_value=_mgmt_value, slice_i=0.55)","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:28:04.277982Z","iopub.execute_input":"2021-09-01T08:28:04.278677Z","iopub.status.idle":"2021-09-01T08:28:07.018663Z","shell.execute_reply.started":"2021-09-01T08:28:04.278629Z","shell.execute_reply":"2021-09-01T08:28:07.017509Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### With MGMT = 1","metadata":{}},{"cell_type":"code","source":"list1=[184,315,155,228]\nfor i in list1:\n _brats21id = train_df.iloc[i][\"BraTS21ID\"]\n _mgmt_value = train_df.iloc[i][\"MGMT_value\"]\n visualize_sample(brats21id=_brats21id, mgmt_value=_mgmt_value, slice_i=0.55)","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:28:07.021412Z","iopub.execute_input":"2021-09-01T08:28:07.021928Z","iopub.status.idle":"2021-09-01T08:28:09.726299Z","shell.execute_reply.started":"2021-09-01T08:28:07.021880Z","shell.execute_reply":"2021-09-01T08:28:09.725248Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To better understand the representation of these MRI images, we can also **create an animation to visualize the sequence of images** of a certain category for a given patient.","metadata":{}},{"cell_type":"code","source":"def create_animation(ims):\n fig = plt.figure(figsize=(10, 10))\n plt.axis('off')\n im = plt.imshow(ims[0], cmap=\"gray\")\n\n def animate_func(i):\n im.set_array(ims[i])\n #return [im]\n return animation.FuncAnimation(fig, animate_func, frames = len(ims), interval = 1000//4)\n\ndef load_dicom_line(path):\n t_paths = sorted(\n glob.glob(os.path.join(path, \"*\")), \n key=lambda x: int(x[:-4].split(\"-\")[-1]),\n )\n images = []\n for filename in t_paths:\n data = load_dicom(filename)\n if data.max() == 0:\n continue\n images.append(data)\n \n return images","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-09-01T08:28:09.727991Z","iopub.execute_input":"2021-09-01T08:28:09.728499Z","iopub.status.idle":"2021-09-01T08:28:09.738137Z","shell.execute_reply.started":"2021-09-01T08:28:09.728437Z","shell.execute_reply":"2021-09-01T08:28:09.736923Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"images = load_dicom_line(data_directory+\"/train/00176/FLAIR\")\nanm_FLAIR0=create_animation(images)\nimages = load_dicom_line(data_directory+\"/train/00176/T1w\")\nanm_T1W0=create_animation(images)\n\n\nimages = load_dicom_line(data_directory+\"/train/00184/FLAIR\")\nanm_FLAIR1=create_animation(images)\nimages = load_dicom_line(data_directory+\"/train/00184/T1w\")\nanm_T1W1=create_animation(images)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-09-01T08:28:09.740226Z","iopub.execute_input":"2021-09-01T08:28:09.741191Z","iopub.status.idle":"2021-09-01T08:28:14.186622Z","shell.execute_reply.started":"2021-09-01T08:28:09.741144Z","shell.execute_reply":"2021-09-01T08:28:14.185592Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### With MGMT = 0\n### *FLAIR type :*","metadata":{}},{"cell_type":"code","source":"anm_FLAIR0","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:28:14.190273Z","iopub.execute_input":"2021-09-01T08:28:14.190639Z","iopub.status.idle":"2021-09-01T08:28:20.215213Z","shell.execute_reply.started":"2021-09-01T08:28:14.190607Z","shell.execute_reply":"2021-09-01T08:28:20.208618Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### *T1w type :*","metadata":{}},{"cell_type":"code","source":"anm_T1W0","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:28:20.217099Z","iopub.execute_input":"2021-09-01T08:28:20.217542Z","iopub.status.idle":"2021-09-01T08:28:38.764716Z","shell.execute_reply.started":"2021-09-01T08:28:20.217499Z","shell.execute_reply":"2021-09-01T08:28:38.763029Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### With MGMT = 1\n### *FLAIR type :*","metadata":{}},{"cell_type":"code","source":"anm_FLAIR1","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:28:38.766544Z","iopub.execute_input":"2021-09-01T08:28:38.767188Z","iopub.status.idle":"2021-09-01T08:28:44.853831Z","shell.execute_reply.started":"2021-09-01T08:28:38.767142Z","shell.execute_reply":"2021-09-01T08:28:44.846577Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### *T1w type :*","metadata":{}},{"cell_type":"code","source":"anm_T1W1","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:28:44.855726Z","iopub.execute_input":"2021-09-01T08:28:44.856139Z","iopub.status.idle":"2021-09-01T08:29:02.872919Z","shell.execute_reply.started":"2021-09-01T08:28:44.856096Z","shell.execute_reply":"2021-09-01T08:29:02.868589Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Distribution of images in MRI types :\nWe are now going to check the number of DCM files to check if their number is the same for each category and for each patient. For that, we will complete a copy of train.csv with the calculated informations:","metadata":{}},{"cell_type":"code","source":"mri_types = ['FLAIR','T1w','T1wCE','T2w']\n\ntrain_dataset = train_df.copy()\nfor scan in mri_types:\n train_dataset[scan + \"_count\"] = [\n len(os.listdir(data_directory + \"/train/\"\n + str(p) \n + \"/\" + scan))\n for p in train_dataset.BraTS21ID5]","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:29:02.875151Z","iopub.execute_input":"2021-09-01T08:29:02.875656Z","iopub.status.idle":"2021-09-01T08:29:57.148186Z","shell.execute_reply.started":"2021-09-01T08:29:02.875607Z","shell.execute_reply":"2021-09-01T08:29:57.147114Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize = (25,40))\nfor i, scan in enumerate(mri_types):\n ax = plt.subplot(4,1,i+1)\n plt.xticks(rotation=70)\n sns.countplot(x=train_dataset[scan + \"_count\"], ax=ax)\n ax.set_xlabel(\"Mean of number of MRI in scan folder\")\n ax.set_ylabel(\"Count of patients\")\n ax.set_title(\"Distribution of number of DCM file in {} scans\".format(scan),\n fontsize=18, color=\"#0b0a2d\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:29:57.150103Z","iopub.execute_input":"2021-09-01T08:29:57.150555Z","iopub.status.idle":"2021-09-01T08:30:08.881959Z","shell.execute_reply.started":"2021-09-01T08:29:57.150512Z","shell.execute_reply":"2021-09-01T08:30:08.880553Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Note that some values for each scan category are over-represented. On the other hand, the span ranges of the counters are important. This may be due, for example, to the use of different X-ray machines ...\n\nHowever, if we only consider patients with maximum values, the amount of data available may be too low for a complex machine learning algorithm.","metadata":{}},{"cell_type":"code","source":"train_dataset[(train_dataset[\"FLAIR_count\"] == int(train_dataset[\"FLAIR_count\"].mode()))\n & (train_dataset[\"T1w_count\"] == int(train_dataset[\"T1w_count\"].mode()))\n & (train_dataset[\"T1wCE_count\"] == int(train_dataset[\"T1wCE_count\"].mode()))\n & (train_dataset[\"T2w_count\"] == int(train_dataset[\"T2w_count\"].mode()))]","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:30:08.883974Z","iopub.execute_input":"2021-09-01T08:30:08.884509Z","iopub.status.idle":"2021-09-01T08:30:08.930770Z","shell.execute_reply.started":"2021-09-01T08:30:08.884458Z","shell.execute_reply":"2021-09-01T08:30:08.929593Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There are only 151 patients whose scans contain the maximum of DCM images out of the 585 at the start. We therefore keep the entire dataset for the moment. **Lets have a look to the test dataset** :","metadata":{}},{"cell_type":"code","source":"test_dataset = test.copy()\nfor scan in mri_types:\n test_dataset[scan + \"_count\"] = [\n len(os.listdir(data_directory + \"/test/\"\n + str(p)\n + \"/\" + scan))\n for p in test_dataset.BraTS21ID5]","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:30:08.932369Z","iopub.execute_input":"2021-09-01T08:30:08.933099Z","iopub.status.idle":"2021-09-01T08:30:12.495884Z","shell.execute_reply.started":"2021-09-01T08:30:08.933023Z","shell.execute_reply":"2021-09-01T08:30:12.494816Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize = (25,40))\nfor i, scan in enumerate(mri_types):\n ax = plt.subplot(4,1,i+1)\n plt.xticks(rotation=70)\n sns.countplot(x=test_dataset[scan + \"_count\"], ax=ax)\n ax.set_xlabel(\"Mean of number of MRI in scan folder\")\n ax.set_ylabel(\"Count of patients\")\n ax.set_title(\"Distribution of number of DCM file in TEST {} scans\".format(scan),\n fontsize=18, color=\"#0b0a2d\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:30:12.497579Z","iopub.execute_input":"2021-09-01T08:30:12.498125Z","iopub.status.idle":"2021-09-01T08:30:16.299204Z","shell.execute_reply.started":"2021-09-01T08:30:12.498059Z","shell.execute_reply":"2021-09-01T08:30:16.298143Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3D CNN from scratch on unique MRI Type\n\nIn this section, a convolutional neural network will be trained on a single type of MRI scan. We will take a **sequence of 64 consecutive images in 80x80 pixels** which will allow the scanners to be treated like videos.","metadata":{}},{"cell_type":"code","source":"# Fixe MRI type to FLAIR\nmri_types_id = 0 # 0,1,2,3\n\n# Initial parameters\nIMAGE_SIZE = 80\nNUM_IMAGES = 64\nBATCH_SIZE = 4\n\nnum_folds = 5\nSelected_fold = 1 #1,2,3,4,5 ","metadata":{"papermill":{"duration":0.056368,"end_time":"2021-08-27T14:49:53.183399","exception":false,"start_time":"2021-08-27T14:49:53.127031","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-09-01T08:30:16.300960Z","iopub.execute_input":"2021-09-01T08:30:16.301458Z","iopub.status.idle":"2021-09-01T08:30:16.307515Z","shell.execute_reply.started":"2021-09-01T08:30:16.301381Z","shell.execute_reply":"2021-09-01T08:30:16.305991Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Preprocessing params\n# Scale for image crop\nSCALE = .8\n# Tile size for CLAHE equalizer\nTILE_SIZE = (8,8)\n# H param for denoising filter\nH_PARAM = 10","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:30:16.309513Z","iopub.execute_input":"2021-09-01T08:30:16.309982Z","iopub.status.idle":"2021-09-01T08:30:16.320744Z","shell.execute_reply.started":"2021-09-01T08:30:16.309938Z","shell.execute_reply":"2021-09-01T08:30:16.319621Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df[\"Fold\"]=\"train\"\ntrain_df.head(5)","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:30:16.327213Z","iopub.execute_input":"2021-09-01T08:30:16.327538Z","iopub.status.idle":"2021-09-01T08:30:16.343551Z","shell.execute_reply.started":"2021-09-01T08:30:16.327507Z","shell.execute_reply":"2021-09-01T08:30:16.342321Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test.head(3)","metadata":{"papermill":{"duration":0.038589,"end_time":"2021-08-27T14:49:53.242402","exception":false,"start_time":"2021-08-27T14:49:53.203813","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-09-01T08:30:16.345835Z","iopub.execute_input":"2021-09-01T08:30:16.346154Z","iopub.status.idle":"2021-09-01T08:30:16.359415Z","shell.execute_reply.started":"2021-09-01T08:30:16.346124Z","shell.execute_reply":"2021-09-01T08:30:16.357823Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Images preprocessing","metadata":{}},{"cell_type":"code","source":"def mri_preprocessing(img, scale=SCALE, img_size=IMAGE_SIZE, tile_size=TILE_SIZE, h_param=H_PARAM):\n # Crop image\n center_x, center_y = img.shape[1] / 2, img.shape[0] / 2\n width_scaled, height_scaled = img.shape[1] * scale, img.shape[0] * scale\n left_x, right_x = center_x - width_scaled / 2, center_x + width_scaled / 2\n top_y, bottom_y = center_y - height_scaled / 2, center_y + height_scaled / 2\n img = img[int(top_y):int(bottom_y), int(left_x):int(right_x)]\n img = cv2.resize(img, (img_size, img_size))\n \n # CLAHE Equalizer\n img = np.uint8(cv2.normalize(img, None, 0, 255, cv2.NORM_MINMAX))\n img = cv2.equalizeHist(img)\n clahe = cv2.createCLAHE(clipLimit=2.0, tileGridSize=tile_size)\n img = clahe.apply(np.uint8(img))\n \n # Non-local mean filter denoising\n img = cv2.fastNlMeansDenoising(\n src=img,\n dst=None,\n h=h_param,\n templateWindowSize=7,\n searchWindowSize=21)\n \n return img","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:30:16.361962Z","iopub.execute_input":"2021-09-01T08:30:16.362657Z","iopub.status.idle":"2021-09-01T08:30:16.376747Z","shell.execute_reply.started":"2021-09-01T08:30:16.362611Z","shell.execute_reply":"2021-09-01T08:30:16.375517Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_img_path = ''.join([data_directory, '/train/00176/FLAIR/Image-25.dcm'])\nsample_img = pydicom.dcmread(sample_img_path)\nsample_img = sample_img.pixel_array\n\nfig = plt.figure(figsize=(12,6))\nax = plt.subplot(1,2,1)\nax.imshow(sample_img, cmap=\"gray\")\nax.set_title(\"Original image\")\nax1 = plt.subplot(1,2,2)\nimg_preproc = mri_preprocessing(sample_img)\nax1.imshow(img_preproc, cmap=\"gray\")\nax1.set_title(\"Preproceced image\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:30:16.378179Z","iopub.execute_input":"2021-09-01T08:30:16.378725Z","iopub.status.idle":"2021-09-01T08:30:16.947775Z","shell.execute_reply.started":"2021-09-01T08:30:16.378680Z","shell.execute_reply":"2021-09-01T08:30:16.946488Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Functions to load images\nWe are going to use 2 functions: one to load the Dicom images, then a second to load a sequence of x images distributed evenly on each side of the central image. These functions will then be used in a custom image generator.","metadata":{"papermill":{"duration":0.019507,"end_time":"2021-08-27T14:49:53.283203","exception":false,"start_time":"2021-08-27T14:49:53.263696","status":"completed"},"tags":[]}},{"cell_type":"code","source":"def load_dicom_image(path, img_size=IMAGE_SIZE, voi_lut=True, rotate=0, preproc=True):\n dicom = pydicom.read_file(path)\n data = dicom.pixel_array\n if voi_lut:\n data = apply_voi_lut(dicom.pixel_array, dicom)\n else:\n data = dicom.pixel_array\n \n if rotate > 0:\n rot_choices = [0, cv2.ROTATE_90_CLOCKWISE, cv2.ROTATE_90_COUNTERCLOCKWISE, cv2.ROTATE_180]\n data = cv2.rotate(data, rot_choices[rotate])\n \n if (preproc==True):\n data = mri_preprocessing(data)\n \n else:\n data = cv2.resize(data, (img_size, img_size))\n \n return data\n\n\ndef load_dicom_images_3d(scan_id, num_imgs=NUM_IMAGES, img_size=IMAGE_SIZE, mri_type=mri_types[mri_types_id], split=\"train\", rotate=0):\n\n files = sorted(glob.glob(f\"{data_directory}/{split}/{scan_id}/{mri_type}/*.dcm\"), \n key=lambda var:[int(x) if x.isdigit() else x for x in re.findall(r'[^0-9]|[0-9]+', var)])\n\n middle = len(files)//2\n num_imgs2 = num_imgs//2\n p1 = max(0, middle - num_imgs2)\n p2 = min(len(files), middle + num_imgs2)\n img3d = np.stack([load_dicom_image(f, rotate=rotate) for f in files[p1:p2]]).T \n if img3d.shape[-1] < num_imgs:\n n_zero = np.zeros((img_size, img_size, num_imgs - img3d.shape[-1]))\n img3d = np.concatenate((img3d, n_zero), axis = -1)\n \n if np.min(img3d) < np.max(img3d):\n img3d = img3d - np.min(img3d)\n img3d = img3d / np.max(img3d)\n \n return np.expand_dims(img3d,0)","metadata":{"papermill":{"duration":1.274245,"end_time":"2021-08-27T14:49:54.577219","exception":false,"start_time":"2021-08-27T14:49:53.302974","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-09-01T08:30:16.949396Z","iopub.execute_input":"2021-09-01T08:30:16.949891Z","iopub.status.idle":"2021-09-01T08:30:16.964251Z","shell.execute_reply.started":"2021-09-01T08:30:16.949832Z","shell.execute_reply":"2021-09-01T08:30:16.963055Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can test these functions on a test sequence :","metadata":{}},{"cell_type":"code","source":"a = load_dicom_images_3d(\"00000\")\nprint(a.shape)\nprint(np.min(a), np.max(a), np.mean(a), np.median(a))\nimage = a[0]\nprint(\"Dimension of the CT scan is:\", image.shape)\nplt.imshow(np.squeeze(image[:, :, 31]), cmap=\"gray\")","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:30:16.966139Z","iopub.execute_input":"2021-09-01T08:30:16.966935Z","iopub.status.idle":"2021-09-01T08:30:18.947874Z","shell.execute_reply.started":"2021-09-01T08:30:16.966884Z","shell.execute_reply":"2021-09-01T08:30:18.946758Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_slices(num_rows, num_columns, width, height, data):\n \"\"\"Plot a montage of 50 CT slices\"\"\"\n data = np.rot90(np.array(data)) \n data = np.transpose(data)\n data = np.reshape(data, (num_rows, num_columns, width, height))\n rows_data, columns_data = data.shape[0], data.shape[1]\n heights = [slc[0].shape[0] for slc in data]\n widths = [slc.shape[1] for slc in data[0]]\n fig_width = 12.0\n fig_height = fig_width * sum(heights) / sum(widths)\n f, axarr = plt.subplots(\n rows_data,\n columns_data,\n figsize=(fig_width, fig_height),\n gridspec_kw={\"height_ratios\": heights},\n )\n for i in range(rows_data):\n for j in range(columns_data):\n axarr[i, j].imshow(data[i][j], cmap=\"gray\")\n axarr[i, j].axis(\"off\")\n plt.subplots_adjust(wspace=0, hspace=0, left=0, right=1, bottom=0, top=1)\n plt.show()","metadata":{"papermill":{"duration":2.354203,"end_time":"2021-08-27T14:49:56.952892","exception":false,"start_time":"2021-08-27T14:49:54.598689","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-09-01T08:30:18.949701Z","iopub.execute_input":"2021-09-01T08:30:18.950126Z","iopub.status.idle":"2021-09-01T08:30:18.962228Z","shell.execute_reply.started":"2021-09-01T08:30:18.950082Z","shell.execute_reply":"2021-09-01T08:30:18.960866Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Visualize montage of slices.\n# 5 rows and 10 columns for 50 slices of the CT scan.\nplot_slices(5, 10, IMAGE_SIZE, IMAGE_SIZE, image[:, :, :50])","metadata":{"execution":{"iopub.status.busy":"2021-09-01T08:37:19.080557Z","iopub.execute_input":"2021-09-01T08:37:19.080966Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Define Folds","metadata":{"papermill":{"duration":0.026541,"end_time":"2021-08-27T14:49:57.006206","exception":false,"start_time":"2021-08-27T14:49:56.979665","status":"completed"},"tags":[]}},{"cell_type":"code","source":"from sklearn.model_selection import KFold,StratifiedKFold\nsfolder = StratifiedKFold(n_splits=5,random_state=13,shuffle=True)\nX = train_df[['BraTS21ID']]\ny = train_df[['MGMT_value']]\n\nfold_no = 1\nfor train, valid in sfolder.split(X,y):\n if fold_no==Selected_fold:\n train_df.loc[valid, \"Fold\"] = \"valid\"\n fold_no += 1","metadata":{"papermill":{"duration":0.043641,"end_time":"2021-08-27T14:49:57.078955","exception":false,"start_time":"2021-08-27T14:49:57.035314","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_train=train_df[train_df.Fold==\"train\"]\ndf_valid=train_df[train_df.Fold==\"valid\"]\nprint(\"df_train=\",len(df_train),\"-- df_valid=\",len(df_valid))","metadata":{"papermill":{"duration":0.038848,"end_time":"2021-08-27T14:49:57.144203","exception":false,"start_time":"2021-08-27T14:49:57.105355","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Keras custom Data Generator\nThanks to the Sequence module of the Keras library, we are going to create a personalized image generator. This will prevent us from creating Numpy arrays or Tensors containing all the sequences which would quickly overload the memory.","metadata":{"papermill":{"duration":0.027252,"end_time":"2021-08-27T14:49:57.198451","exception":false,"start_time":"2021-08-27T14:49:57.171199","status":"completed"},"tags":[]}},{"cell_type":"code","source":"from keras.utils import Sequence\nclass Dataset(Sequence):\n def __init__(self,df,is_train=True,batch_size=BATCH_SIZE,shuffle=True):\n self.idx = df[\"BraTS21ID\"].values\n self.paths = df[\"BraTS21ID5\"].values\n self.y = df[\"MGMT_value\"].values\n self.is_train = is_train\n self.batch_size = batch_size\n self.shuffle = shuffle\n def __len__(self):\n return math.ceil(len(self.idx)/self.batch_size)\n \n def __getitem__(self,ids):\n id_path= self.paths[ids]\n batch_paths = self.paths[ids * self.batch_size:(ids + 1) * self.batch_size]\n \n if self.y is not None:\n batch_y = self.y[ids * self.batch_size: (ids + 1) * self.batch_size]\n \n if self.is_train:\n list_x = [load_dicom_images_3d(x,split=\"train\") for x in batch_paths]\n batch_X = np.stack(list_x, axis=4)\n return batch_X,batch_y\n else:\n list_x = load_dicom_images_3d(id_path,split=\"test\")#str(scan_id).zfill(5)\n batch_X = np.stack(list_x)\n return batch_X\n \n def on_epoch_end(self):\n if self.shuffle and self.is_train:\n ids_y = list(zip(self.idx, self.y))\n shuffle(ids_y)\n self.idx, self.y = list(zip(*ids_y))","metadata":{"papermill":{"duration":0.040884,"end_time":"2021-08-27T14:49:57.268911","exception":false,"start_time":"2021-08-27T14:49:57.228027","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_dataset = Dataset(df_train,batch_size=BATCH_SIZE)\nvalid_dataset = Dataset(df_valid,batch_size=BATCH_SIZE)","metadata":{"papermill":{"duration":0.035076,"end_time":"2021-08-27T14:49:57.331697","exception":false,"start_time":"2021-08-27T14:49:57.296621","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Once the generators are created, we can project an image to check:","metadata":{}},{"cell_type":"code","source":"for i in range(1):\n images, label = train_dataset[i]\n print(\"Dimension of the CT scan is:\", images.shape)\n print(\"label=\",label)\n plt.imshow(images[0,:,:,32,0], cmap=\"gray\")\n plt.show()","metadata":{"papermill":{"duration":3.132042,"end_time":"2021-08-27T14:50:00.490834","exception":false,"start_time":"2021-08-27T14:49:57.358792","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Define CNN Model","metadata":{"papermill":{"duration":0.028384,"end_time":"2021-08-27T14:50:00.54865","exception":false,"start_time":"2021-08-27T14:50:00.520266","status":"completed"},"tags":[]}},{"cell_type":"code","source":"def get_model(width=IMAGE_SIZE, height=IMAGE_SIZE, depth=NUM_IMAGES):\n \"\"\"Build a 3D convolutional neural network model.\"\"\"\n\n inputs = keras.Input((width, height, depth, 1))\n \n x = layers.Conv3D(filters=32, kernel_size=3, activation=\"relu\")(inputs)\n x = layers.MaxPool3D(pool_size=2)(x)\n x = layers.BatchNormalization()(x)\n \n x = layers.Conv3D(filters=32, kernel_size=3, activation=\"relu\")(inputs)\n x = layers.MaxPool3D(pool_size=2)(x)\n x = layers.BatchNormalization()(x)\n \n x = layers.Conv3D(filters=64, kernel_size=3, activation=\"relu\")(inputs)\n x = layers.MaxPool3D(pool_size=2)(x)\n x = layers.BatchNormalization()(x)\n x = layers.Dropout(0.01)(x)\n \n x = layers.Conv3D(filters=128, kernel_size=3, activation=\"relu\")(x)\n x = layers.MaxPool3D(pool_size=2)(x)\n x = layers.BatchNormalization()(x)\n x = layers.Dropout(0.02)(x)\n\n x = layers.Conv3D(filters=256, kernel_size=3, activation=\"relu\")(x)\n x = layers.MaxPool3D(pool_size=2)(x)\n x = layers.BatchNormalization()(x)\n x = layers.Dropout(0.03)(x)\n\n x = layers.Conv3D(filters=512, kernel_size=3, activation=\"relu\")(x)\n x = layers.MaxPool3D(pool_size=2)(x)\n x = layers.BatchNormalization()(x)\n x = layers.Dropout(0.04)(x)\n\n x = layers.GlobalAveragePooling3D()(x)\n x = layers.Dense(units=1024, activation=\"relu\")(x)\n x = layers.Dropout(0.08)(x)\n\n outputs = layers.Dense(units=1, activation=\"sigmoid\")(x)\n\n # Define the model.\n model = keras.Model(inputs, outputs, name=\"3dcnn\")\n\n return model","metadata":{"papermill":{"duration":2.224738,"end_time":"2021-08-27T14:50:02.801731","exception":false,"start_time":"2021-08-27T14:50:00.576993","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Build model.\nmodel = get_model(width=IMAGE_SIZE, height=IMAGE_SIZE, depth=NUM_IMAGES)\nmodel.summary()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Training of the CNN","metadata":{"papermill":{"duration":0.028424,"end_time":"2021-08-27T14:50:02.859855","exception":false,"start_time":"2021-08-27T14:50:02.831431","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# Compile model.\ninitial_learning_rate = 0.0001\nlr_schedule = keras.optimizers.schedules.ExponentialDecay(\n initial_learning_rate, decay_steps=100000, decay_rate=0.96, staircase=True\n)\nmodel.compile(\n loss=\"binary_crossentropy\",\n optimizer=keras.optimizers.Adam(learning_rate=lr_schedule),\n metrics=[AUC(name='auc'),\"acc\"],\n)\n# Define callbacks.\nmodel_save = ModelCheckpoint(f'Brain_3d_cls_{mri_types[mri_types_id]}_Fold_{Selected_fold}.h5', \n save_best_only = True, \n monitor = 'val_auc', \n mode = 'max', verbose = 1)\nearly_stop = EarlyStopping(monitor = 'val_auc', \n patience = 5, mode = 'max', verbose = 1,\n restore_best_weights = True)","metadata":{"papermill":{"duration":6856.572755,"end_time":"2021-08-27T16:44:19.46109","exception":false,"start_time":"2021-08-27T14:50:02.888335","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Train the model, doing validation at the end of each epoch\nepochs = 50\nmodel.fit(\n train_dataset,\n validation_data=valid_dataset,\n epochs=epochs,\n shuffle=True,\n verbose=1,\n callbacks = [model_save, early_stop],\n)","metadata":{"execution":{"iopub.status.idle":"2021-09-01T08:36:47.802880Z","shell.execute_reply.started":"2021-09-01T08:30:32.808663Z","shell.execute_reply":"2021-09-01T08:36:47.799238Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Visualizing model performance","metadata":{"papermill":{"duration":0.818772,"end_time":"2021-08-27T16:44:21.111094","exception":false,"start_time":"2021-08-27T16:44:20.292322","status":"completed"},"tags":[]}},{"cell_type":"code","source":"fig, ax = plt.subplots(1, 3, figsize=(20, 7))\nax = ax.ravel()\n\nfor i, metric in enumerate([\"acc\",\"auc\",\"loss\"]):\n ax[i].plot(model.history.history[metric])\n ax[i].plot(model.history.history[\"val_\" + metric])\n ax[i].set_title(\"Model {}\".format(metric))\n ax[i].set_xlabel(\"epochs\")\n ax[i].set_ylabel(metric)\n ax[i].legend([\"train\", \"val\"])","metadata":{"papermill":{"duration":1.393824,"end_time":"2021-08-27T16:44:23.320494","exception":false,"start_time":"2021-08-27T16:44:21.92667","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-09-01T08:36:47.804705Z","iopub.status.idle":"2021-09-01T08:36:47.806127Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Make predictions with trained CNN Model","metadata":{"papermill":{"duration":0.823394,"end_time":"2021-08-27T16:44:24.972668","exception":false,"start_time":"2021-08-27T16:44:24.149274","status":"completed"},"tags":[]}},{"cell_type":"code","source":"test_dataset = Dataset(test,is_train=False,batch_size=1)\n\n\nfor i in range(1):\n image = test_dataset[i]\n print(\"Dimension of the CT scan is:\", image.shape)\n plt.imshow(image[0,:,:, 32], cmap=\"gray\")\n plt.show()","metadata":{"papermill":{"duration":1.773563,"end_time":"2021-08-27T16:44:27.807698","exception":false,"start_time":"2021-08-27T16:44:26.034135","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-09-01T08:36:47.810974Z","iopub.status.idle":"2021-09-01T08:36:47.811762Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"preds = model.predict(test_dataset)\npreds = preds.reshape(-1)","metadata":{"papermill":{"duration":58.673124,"end_time":"2021-08-27T16:45:27.589709","exception":false,"start_time":"2021-08-27T16:44:28.916585","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-09-01T08:36:47.814400Z","iopub.status.idle":"2021-09-01T08:36:47.815388Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"preds","metadata":{"papermill":{"duration":0.827581,"end_time":"2021-08-27T16:45:29.245228","exception":false,"start_time":"2021-08-27T16:45:28.417647","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-09-01T08:36:47.819010Z","iopub.status.idle":"2021-09-01T08:36:47.819876Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Submission of result predictions","metadata":{}},{"cell_type":"code","source":"submission = pd.DataFrame({'BraTS21ID':test['BraTS21ID'],'MGMT_value':preds})\nsubmission","metadata":{"papermill":{"duration":0.868896,"end_time":"2021-08-27T16:45:30.993674","exception":false,"start_time":"2021-08-27T16:45:30.124778","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-09-01T08:36:47.830653Z","iopub.status.idle":"2021-09-01T08:36:47.831507Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission.to_csv('submission.csv',index=False)","metadata":{"papermill":{"duration":1.305978,"end_time":"2021-08-27T16:45:33.396465","exception":false,"start_time":"2021-08-27T16:45:32.090487","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-09-01T08:36:47.838535Z","iopub.status.idle":"2021-09-01T08:36:47.841296Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(5, 5))\nplt.hist(submission[\"MGMT_value\"]);","metadata":{"papermill":{"duration":0.99614,"end_time":"2021-08-27T16:45:35.239675","exception":false,"start_time":"2021-08-27T16:45:34.243535","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-09-01T08:36:47.842836Z","iopub.status.idle":"2021-09-01T08:36:47.848191Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# References","metadata":{"papermill":{"duration":0.87261,"end_time":"2021-08-27T16:45:38.578167","exception":false,"start_time":"2021-08-27T16:45:37.705557","status":"completed"},"tags":[]}},{"cell_type":"markdown","source":"1. https://keras.io/examples/vision/3D_image_classification/\n1. https://www.kaggle.com/rluethy/efficientnet3d-with-one-mri-type\n","metadata":{"papermill":{"duration":0.873657,"end_time":"2021-08-27T16:45:40.299542","exception":false,"start_time":"2021-08-27T16:45:39.425885","status":"completed"},"tags":[]}}]}