{"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":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-10-07T13:18:52.878586Z","iopub.execute_input":"2021-10-07T13:18:52.879196Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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\nimport tensorflow as tf\nimport matplotlib.pyplot as plt\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 tensorflow.keras.metrics import AUC","metadata":{"execution":{"iopub.status.busy":"2021-10-07T14:17:10.062172Z","iopub.execute_input":"2021-10-07T14:17:10.062431Z","iopub.status.idle":"2021-10-07T14:17:15.536568Z","shell.execute_reply.started":"2021-10-07T14:17:10.062404Z","shell.execute_reply":"2021-10-07T14:17:15.535865Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#global variables,initialisations\n\ndata_directory = '../input/rsna-miccai-brain-tumor-radiogenomic-classification'\n\n \nmri_types = ['FLAIR','T1w','T1wCE','T2w']\nmri_types_id=0 # 0,1,2,3\n\n#3-D image parametrs\nIMAGE_SIZE = 128\nNUM_IMAGES = 64\nBATCH_SIZE= 4\n","metadata":{"execution":{"iopub.status.busy":"2021-10-07T14:17:19.789168Z","iopub.execute_input":"2021-10-07T14:17:19.789436Z","iopub.status.idle":"2021-10-07T14:17:19.793709Z","shell.execute_reply.started":"2021-10-07T14:17:19.789408Z","shell.execute_reply":"2021-10-07T14:17:19.792988Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#dataframe, formatting columns\ntrain_df = pd.read_csv(data_directory+\"/train_labels.csv\")\ntrain_df['BraTS21ID5'] = [format(x, '05d') for x in train_df.BraTS21ID] #formatting to display \"a\" as \"00000a\"\ntrain_df[\"Fold\"]=\"train\"\ntrain_df.head(20) #visualising 20 entries","metadata":{"execution":{"iopub.status.busy":"2021-10-07T14:18:07.140467Z","iopub.execute_input":"2021-10-07T14:18:07.141465Z","iopub.status.idle":"2021-10-07T14:18:07.181371Z","shell.execute_reply.started":"2021-10-07T14:18:07.141415Z","shell.execute_reply":"2021-10-07T14:18:07.18056Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#loading dicom images, resizing, processing them\n\ndef load_dicom_image(path, img_size=IMAGE_SIZE, voi_lut=True, rotate=0):\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    data = cv2.resize(data, (img_size, img_size)) #resizing images\n    \n    return data #pixel array\n","metadata":{"execution":{"iopub.status.busy":"2021-10-07T14:17:22.854279Z","iopub.execute_input":"2021-10-07T14:17:22.854871Z","iopub.status.idle":"2021-10-07T14:17:22.861082Z","shell.execute_reply.started":"2021-10-07T14:17:22.854834Z","shell.execute_reply":"2021-10-07T14:17:22.860206Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#creating 3-D images\n\ndef load_dicom_images_3d(scan_id, num_imgs=NUM_IMAGES, img_size=IMAGE_SIZE, mri_type=mri_types[mri_types_id], \n                         split=\"train\", rotate=0):\n    \n    list_paths=glob.glob(f\"{data_directory}/{split}/{scan_id}/{mri_type}/*.dcm\")#formatted string\n                                                                        #all files with path ../input../train/{id}/FLAIR.dcm\n                                                                        #accessed, paths stored in list_paths\n            \n    k=lambda var:[int(x) if x.isdigit() else x for x in re.findall(r'[^0-9]|[0-9]+', var)]\n    \n    files = sorted(list_paths,key=k)\n\n    middle = len(files)//2\n    ni2= num_imgs//2\n    \n    p1 = max(0, middle - ni2)\n    p2 = min(len(files), middle + ni2)\n    \n    #values selected to select 64 images from the middle of the files.\n    #for eg, for patient 000000, len(files)=400\n    #so p1=168, p2= 232, and the next line, files[168:232] are selected, i.e 64 images , from the middle\n    \n    img3d = np.stack([load_dicom_image(f, rotate=rotate) for f in files[p1:p2]])\n    img3d=img3d.T\n    #stacking along axis=0 (default), then transposed. img3d is a (256,256,64) array\n    \n    if img3d.shape[-1] < num_imgs: #if we cant afford to have 64 images. Suppose we have x. Then our array is of size \n                                   #(256,256,x)\n        n_zero = np.zeros((img_size, img_size, num_imgs - img3d.shape[-1])) #zeroes array of size (256,256,64-x)\n        \n        img3d = np.concatenate((img3d,  n_zero), axis = -1) #concatenating makes the array size (256,256,64)\n        \n    if np.min(img3d) < np.max(img3d):\n        #standardisation (making every image pixel having a value in [0,1])\n        #\"min\" and \"max\" pixel mapped to 0 and 1.\n        \n        img3d = img3d - np.min(img3d)\n        img3d = img3d / np.max(img3d)\n            \n    return np.expand_dims(img3d,0) #adding an extra dimension to img3d\n                                   #3-D images are represented using 4D arrays\n                                   #3image dimensions(height,width,depth(no.of images)) + 1 for no. of colour channels \n                        \n","metadata":{"execution":{"iopub.status.busy":"2021-10-07T14:17:25.97426Z","iopub.execute_input":"2021-10-07T14:17:25.974545Z","iopub.status.idle":"2021-10-07T14:17:25.987417Z","shell.execute_reply.started":"2021-10-07T14:17:25.974517Z","shell.execute_reply":"2021-10-07T14:17:25.986618Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#testing the function for one example\n\na = load_dicom_images_3d(\"00000\") #patient scan ID='00000'\nprint(a.shape)\n#(1,256,256,64) : array describing series of 64 images (256x256), having only 1 colour channel\nimage = a[0] #is the array (256,256,64), i.e series of 64 images\nprint(\"Dimension of the CT scan is:\", image.shape)\n\nplt.imshow(np.squeeze(image[:, :, 31]), cmap=\"gray\") #printing the 32nd image of the series as an example","metadata":{"execution":{"iopub.status.busy":"2021-10-07T14:17:29.676108Z","iopub.execute_input":"2021-10-07T14:17:29.676403Z","iopub.status.idle":"2021-10-07T14:17:31.226549Z","shell.execute_reply.started":"2021-10-07T14:17:29.676374Z","shell.execute_reply":"2021-10-07T14:17:31.225834Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#splitting into train and validation datasets, stratifying.\n\ndf_train, df_valid = sk_model_selection.train_test_split(\n    train_df, \n    test_size=0.2, #valid dataset is 1/5 of total, i.e 1/5*585=117\n    random_state=12, \n    stratify=train_df[\"MGMT_value\"], #stratifying.\n                            #the train and valid sets will have the same proportion of \"1s\" and \"0s\" as \n                            #the original. For eg, if there were 75% 1s and 25% 0s, then the splits (train and valid)\n                            #will each have 75% 1s and 25% 0s.\n)\n#checking the data frames\n#len(df_train) #468\n#len(df_valid) #117\n#df_valid\n#df_train \n","metadata":{"execution":{"iopub.status.busy":"2021-10-07T14:18:12.098401Z","iopub.execute_input":"2021-10-07T14:18:12.098692Z","iopub.status.idle":"2021-10-07T14:18:12.111746Z","shell.execute_reply.started":"2021-10-07T14:18:12.098663Z","shell.execute_reply":"2021-10-07T14:18:12.110522Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#use of sequences is helpful as it ensures that, per epoch, each input will be trained only once.\n\nfrom tensorflow.keras.utils import Sequence\n\nclass Dataset(Sequence):\n    def __init__(self,df,is_train=True,batch_size=4 #BATCH_SIZE\n                 ,shuffle=True):\n        \n        self.idx = df[\"BraTS21ID\"].values #values 0,2,3..\n        self.paths = df[\"BraTS21ID5\"].values #values 000000,00002,00003..\n        self.y =  df[\"MGMT_value\"].values #MGMT values 0/1\n        \n        #is_train and shuffle are set to true\n        self.is_train = is_train\n        self.shuffle = shuffle\n        \n        self.batch_size = batch_size\n        \n    def __len__(self):\n        \n        #eg if self.idx is the list [0,3,7.....](10 values)\n        #then len attribute will be 10/5=2.\n        \n        return math.ceil(len(self.idx)/self.batch_size)\n   \n    def __getitem__(self,ids):\n        \n        id_path= self.paths[ids] \n        \n      \n        start= ids*self.batch_size\n        end= (ids+1)*self.batch_size\n        \n        batch_paths = self.paths[start:end]   #eg if ids=0, then batch_paths will be self.paths[0:4]\n        \n        if self.y is not None: \n            batch_y = self.y[start: end] #eg if ids=0, then batch_y will be self.y[0:4], i.e MGMT values for \n                                         #first 4 patients\n        \n        if self.is_train: #set to true\n            list_x =  [load_dicom_images_3d(x,split=\"train\") for x in batch_paths] #list of all the arrays corresponding to\n                                                                                   #the images having IDs in batch_paths\n            batch_X = np.stack(list_x, axis=4) #stacking the arrays along a 4th axis\n            return batch_X,batch_y #returning the arrays for the batch, alog with MGMT values\n        \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): #after one epoch completed, we will shuffle around the values for the next epoch\n        \n        if self.shuffle and self.is_train: #both set to true\n            ids_y = list(zip(self.idx, self.y))  #will \"zip\" the two lists into one list , containg entries like\n                                        #[(0,1),(2,1),(3,0),(5,1)...]\n                                        #[(index,MGMT value),(index,MGMT value)..]\n            shuffle(ids_y) #shuffle the contents\n            self.idx, self.y = list(zip(*ids_y)) #unzips. so now the order of values in idx ,y lists will be changed","metadata":{"execution":{"iopub.status.busy":"2021-10-07T14:18:21.442837Z","iopub.execute_input":"2021-10-07T14:18:21.443258Z","iopub.status.idle":"2021-10-07T14:18:21.456833Z","shell.execute_reply.started":"2021-10-07T14:18:21.443224Z","shell.execute_reply":"2021-10-07T14:18:21.454315Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#train and valid datasets\n\n#declaring them as objects of class Dataset.\n#these two will now have all the attributes discussed earlier, like len, getitem,y.etc\n#train and valid dataframes established earlier (by splitting)\n\ntrain_dataset = Dataset(df_train)\nvalid_dataset = Dataset(df_valid)\n\n#exploring:\n#train_dataset.idx\n#valid_dataset.idx\n","metadata":{"execution":{"iopub.status.busy":"2021-10-07T14:18:26.321145Z","iopub.execute_input":"2021-10-07T14:18:26.321855Z","iopub.status.idle":"2021-10-07T14:18:26.329038Z","shell.execute_reply.started":"2021-10-07T14:18:26.321819Z","shell.execute_reply":"2021-10-07T14:18:26.328197Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#testing the datasets , i.e if they succesfully create the 3D image array\n\nfor i in range(2):\n    images, label = train_dataset[i] #calls train_dataset.__getitem__(self, i)\n                                     #images will be the array corresponding to the 4 paths: path no. 4i to path no. 4(i+1),\n                                     #(i.e 4 patients)\n            \n                                     #labels will be the MGMT values for these 4 patients\n            \n    print(\"Dimension of the CT scan is:\", images.shape) #array dimensions\n    print(\"label=\",label)\n    \n    plt.imshow(images[0,:,:,32,0], cmap=\"gray\") #color channel value: 0 (only 1 value allowed,grayscale)\n                                                #height, width: 128,128\n                                                #slice no. 32\n                                                #0 (last value) is the index value for batch_paths list\n                                                \n    plt.show() #first iteration will show 1st image of batch_paths list , 2nd iteration will again show\n                 #first image of batch_paths list, but this list will now be different, since i is now 1","metadata":{"execution":{"iopub.status.busy":"2021-10-07T14:18:30.174092Z","iopub.execute_input":"2021-10-07T14:18:30.176441Z","iopub.status.idle":"2021-10-07T14:18:37.809924Z","shell.execute_reply.started":"2021-10-07T14:18:30.176397Z","shell.execute_reply":"2021-10-07T14:18:37.809223Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#build model based on 3D CNN architecture\n\n\ndef get_model(width=128, height=128, depth=64):\n    \"\"\"Build a 3D convolutional neural network model.\"\"\"\n\n    inputs = keras.Input((width, height, depth, 1))\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\n    x = layers.Conv3D(filters=64, kernel_size=3, activation=\"relu\")(x)\n    x = layers.MaxPool3D(pool_size=2)(x)\n    x = layers.BatchNormalization()(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\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\n    x = layers.GlobalAveragePooling3D()(x)\n    x = layers.Dense(units=512, activation=\"relu\")(x)\n    x = layers.Dropout(0.3)(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    return model\n\n\n# Build model.\nmodel = get_model(width=128, height=128, depth=64)\nmodel.summary()","metadata":{"execution":{"iopub.status.busy":"2021-10-07T14:18:43.745704Z","iopub.execute_input":"2021-10-07T14:18:43.746292Z","iopub.status.idle":"2021-10-07T14:18:45.933384Z","shell.execute_reply.started":"2021-10-07T14:18:43.746258Z","shell.execute_reply":"2021-10-07T14:18:45.931979Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Compile model.\ninitial_learning_rate = 0.0005\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=[\"acc\"],\n)\n\n# Define callbacks.\ncheckpoint_cb = keras.callbacks.ModelCheckpoint(\n    \"Brain_3d_classification.h5\", save_best_only=True\n)\nearly_stopping_cb = keras.callbacks.EarlyStopping(monitor=\"val_acc\", patience=15)\n\n# 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=2,\n    callbacks=[checkpoint_cb, early_stopping_cb],\n)","metadata":{"execution":{"iopub.status.busy":"2021-10-07T14:19:02.079471Z","iopub.execute_input":"2021-10-07T14:19:02.080293Z","iopub.status.idle":"2021-10-07T15:17:50.931831Z","shell.execute_reply.started":"2021-10-07T14:19:02.080255Z","shell.execute_reply":"2021-10-07T15:17:50.928373Z"},"trusted":true},"execution_count":null,"outputs":[]}]}