{"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":"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-15T14:28:41.23441Z","iopub.execute_input":"2021-10-15T14:28:41.234655Z","iopub.status.idle":"2021-10-15T14:28:41.241619Z","shell.execute_reply.started":"2021-10-15T14:28:41.234626Z","shell.execute_reply":"2021-10-15T14:28:41.240704Z"},"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-15T14:28:41.243574Z","iopub.execute_input":"2021-10-15T14:28:41.243899Z","iopub.status.idle":"2021-10-15T14:28:41.25238Z","shell.execute_reply.started":"2021-10-15T14:28:41.243848Z","shell.execute_reply":"2021-10-15T14:28:41.251652Z"},"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-15T14:28:41.254139Z","iopub.execute_input":"2021-10-15T14:28:41.254465Z","iopub.status.idle":"2021-10-15T14:28:41.280735Z","shell.execute_reply.started":"2021-10-15T14:28:41.254433Z","shell.execute_reply":"2021-10-15T14:28:41.279825Z"},"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-15T14:28:41.282909Z","iopub.execute_input":"2021-10-15T14:28:41.283171Z","iopub.status.idle":"2021-10-15T14:28:41.290086Z","shell.execute_reply.started":"2021-10-15T14:28:41.28314Z","shell.execute_reply":"2021-10-15T14:28:41.289143Z"},"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-15T14:28:41.292027Z","iopub.execute_input":"2021-10-15T14:28:41.292303Z","iopub.status.idle":"2021-10-15T14:28:41.306775Z","shell.execute_reply.started":"2021-10-15T14:28:41.292271Z","shell.execute_reply":"2021-10-15T14:28:41.306069Z"},"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-15T14:28:41.308784Z","iopub.execute_input":"2021-10-15T14:28:41.309002Z","iopub.status.idle":"2021-10-15T14:28:41.831608Z","shell.execute_reply.started":"2021-10-15T14:28:41.308979Z","shell.execute_reply":"2021-10-15T14:28:41.830917Z"},"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-15T14:28:41.833048Z","iopub.execute_input":"2021-10-15T14:28:41.833499Z","iopub.status.idle":"2021-10-15T14:28:41.841919Z","shell.execute_reply.started":"2021-10-15T14:28:41.833462Z","shell.execute_reply":"2021-10-15T14:28:41.84121Z"},"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-15T14:28:41.844932Z","iopub.execute_input":"2021-10-15T14:28:41.845211Z","iopub.status.idle":"2021-10-15T14:28:41.85869Z","shell.execute_reply.started":"2021-10-15T14:28:41.845176Z","shell.execute_reply":"2021-10-15T14:28:41.857647Z"},"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-15T14:28:41.860253Z","iopub.execute_input":"2021-10-15T14:28:41.860634Z","iopub.status.idle":"2021-10-15T14:28:41.870467Z","shell.execute_reply.started":"2021-10-15T14:28:41.860599Z","shell.execute_reply":"2021-10-15T14:28:41.869783Z"},"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-15T14:28:41.87182Z","iopub.execute_input":"2021-10-15T14:28:41.872336Z","iopub.status.idle":"2021-10-15T14:28:44.407944Z","shell.execute_reply.started":"2021-10-15T14:28:41.8723Z","shell.execute_reply":"2021-10-15T14:28:44.407271Z"},"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    \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-15T14:28:44.411671Z","iopub.execute_input":"2021-10-15T14:28:44.413673Z","iopub.status.idle":"2021-10-15T14:28:44.541254Z","shell.execute_reply.started":"2021-10-15T14:28:44.413624Z","shell.execute_reply":"2021-10-15T14:28:44.540572Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Compile model.\ninitial_learning_rate = 0.00001\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.\n#checkpoint_dir=\"../input/btw\"\n#checkpoint_dir = \"bt-weights\"\ncheckpoint_path =\"Brain_3d_classification.h5\" \n \n\ncheckpoint_cb = keras.callbacks.ModelCheckpoint(\n    checkpoint_path , save_best_only=True,monitor = 'val_acc', \n                             mode = 'max', verbose = 1\n)\nearly_stopping_cb = keras.callbacks.EarlyStopping(monitor=\"val_acc\", patience=20,mode = 'max', verbose = 1,\n                           restore_best_weights = True)\n\n# Train the model, doing validation at the end of each epoch\nepochs = 60\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)\n\nfig, ax = plt.subplots(1, 2, figsize=(20, 7))\nax = ax.ravel()\n\nfor i, metric in enumerate([\"acc\",\"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":{"execution":{"iopub.status.busy":"2021-10-15T14:28:44.54251Z","iopub.execute_input":"2021-10-15T14:28:44.542772Z","iopub.status.idle":"2021-10-15T15:45:16.71936Z","shell.execute_reply.started":"2021-10-15T14:28:44.542735Z","shell.execute_reply":"2021-10-15T15:45:16.718606Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#!ls {checkpoint_dir}\n\n#preds = model.predict(test_dataset)\n#preds = preds.reshape(-1)\n#preds\nos.get","metadata":{"execution":{"iopub.status.busy":"2021-10-15T15:45:16.721546Z","iopub.execute_input":"2021-10-15T15:45:16.721793Z","iopub.status.idle":"2021-10-15T15:45:16.727914Z","shell.execute_reply.started":"2021-10-15T15:45:16.721759Z","shell.execute_reply":"2021-10-15T15:45:16.726993Z"},"trusted":true},"execution_count":null,"outputs":[]}]}