{"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":"!pip install ../input/kerasapplications \n!pip install ../input/classificationmodelsed/keras-2.6.0-py2.py3-none-any.whl\n!pip install ../input/classificationmodelsed/classification_models_3D-1.0.2-py3-none-any.whl","metadata":{"execution":{"iopub.status.busy":"2021-09-27T17:27:20.617101Z","iopub.execute_input":"2021-09-27T17:27:20.617472Z","iopub.status.idle":"2021-09-27T17:27:44.529179Z","shell.execute_reply.started":"2021-09-27T17:27:20.617434Z","shell.execute_reply":"2021-09-27T17:27:44.528309Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"@InProceedings{RSolovyev_2021_stalled,\n  author = {Solovyev, Roman and Kalinin, Alexandr A. and Gabruseva, Tatiana},\n  title = {3D Convolutional Neural Networks for Stalled Brain Capillary Detection},\n  booktitle = {Arxiv: 2104.01687},\n  month = {April},\n  year = {2021}\n}","metadata":{"execution":{"iopub.status.busy":"2021-09-26T18:35:36.392988Z","iopub.execute_input":"2021-09-26T18:35:36.394073Z","iopub.status.idle":"2021-09-26T18:35:36.520285Z","shell.execute_reply.started":"2021-09-26T18:35:36.393948Z","shell.execute_reply":"2021-09-26T18:35:36.518355Z"}}},{"cell_type":"markdown","source":"# RSNA MICCAI Brain Tumor Radiogenomic Classification using 3D Conv [TF]\nIn this notebook we will learn to train a 3D conv model using transfer learning approach. \nWe will be using all the MRI types from the dataset, and in the inference use a blending based appraoch to predict.\nPlease refer to the inference notebook (TBD) to see how the trained models are used to predict. ","metadata":{}},{"cell_type":"code","source":"import os\nimport re \nimport glob\nimport numpy as np\nimport pandas as pd\nimport cv2\nimport seaborn as sns\nfrom pathlib import Path\nfrom tqdm.notebook import tqdm\nimport warnings\nwarnings.filterwarnings('ignore')\nimport random as rn\nimport matplotlib.pyplot as plt\nimport imageio\nimport pydicom\nimport math\nfrom numpy.random import default_rng\nos.environ['TF_CPP_MIN_LOG_LEVEL'] = '3' \n\n# Deep learning packages\nimport tensorflow as tf\nfrom tensorflow.keras.callbacks import *\nfrom tensorflow.keras.layers import *\nfrom tensorflow.keras.models import *\nfrom tensorflow import keras\n\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.model_selection import StratifiedKFold\nfrom random import shuffle\n\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\nfrom classification_models_3D.tfkeras import Classifiers","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-09-27T17:27:44.532619Z","iopub.execute_input":"2021-09-27T17:27:44.532855Z","iopub.status.idle":"2021-09-27T17:27:44.544088Z","shell.execute_reply.started":"2021-09-27T17:27:44.532825Z","shell.execute_reply":"2021-09-27T17:27:44.543324Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"config = {\n  'images_source_path' : '../input/rsna-miccai-brain-tumor-radiogenomic-classification/train',\n  'test_images_source_path' : '../input/rsna-miccai-brain-tumor-radiogenomic-classification/test',\n  'csv_path': '../input/rsna-miccai-brain-tumor-radiogenomic-classification/train_labels.csv',\n  'data_path': '../input/rsna-miccai-brain-tumor-radiogenomic-classification',\n  'nfolds': 3,\n  'global_seed': 42,\n  'batch_size': 2,\n  'frames_per_seq': 16,\n  'img_size': 128,\n  'learning_rate': 0.0001,\n  'num_epochs': 15,\n  'channels': 3,\n  'scale' : 0.8\n}\n\n# mri_types = ['FLAIR'] \nmri_types = ['FLAIR','T1w','T1wCE','T2w']","metadata":{"execution":{"iopub.status.busy":"2021-09-27T17:38:55.407953Z","iopub.execute_input":"2021-09-27T17:38:55.408288Z","iopub.status.idle":"2021-09-27T17:38:55.416035Z","shell.execute_reply.started":"2021-09-27T17:38:55.408254Z","shell.execute_reply":"2021-09-27T17:38:55.413799Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def set_seed(seed):\n    rn.seed(seed)\n    np.random.seed(seed)\n    tf.compat.v1.random.set_random_seed(seed)\n\nset_seed(config['global_seed'])","metadata":{"execution":{"iopub.status.busy":"2021-09-27T17:38:55.609863Z","iopub.execute_input":"2021-09-27T17:38:55.610528Z","iopub.status.idle":"2021-09-27T17:38:55.686995Z","shell.execute_reply.started":"2021-09-27T17:38:55.610475Z","shell.execute_reply":"2021-09-27T17:38:55.68569Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_data = pd.read_csv(config['csv_path'])\ndf_data[\"folder_name\"] = [format(x, \"05d\") for x in df_data[\"BraTS21ID\"]]\ndf_data[\"folder_path\"] = [os.path.join(config['images_source_path'], x) for x in df_data[\"folder_name\"]]\nskf = StratifiedKFold(n_splits=config['nfolds'], shuffle=True, random_state=config['global_seed'])\nfor index, (train_index, val_index) in enumerate(skf.split(X=df_data.index, y=df_data.MGMT_value)):\n    df_data.loc[val_index, 'fold'] = index\n# df_data = df_data.head(30)\n# data from following patients is invalid as per organizer.\ndf_data = df_data[~df_data.folder_name.isin([\"00109\", \"00123\", \"00709\"])]\nlen(df_data)","metadata":{"execution":{"iopub.status.busy":"2021-09-27T17:38:55.78048Z","iopub.execute_input":"2021-09-27T17:38:55.780718Z","iopub.status.idle":"2021-09-27T17:38:55.80621Z","shell.execute_reply.started":"2021-09-27T17:38:55.780691Z","shell.execute_reply":"2021-09-27T17:38:55.805258Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class Dataset(tf.keras.utils.Sequence):\n    def __init__(self,df,is_train=True,batch_size=config['batch_size'],shuffle=True):\n        self.idx = df[\"BraTS21ID\"].values\n        self.paths = df[\"folder_path\"].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        self.df = df\n        \n    def __len__(self):\n        return math.ceil(len(self.idx)/self.batch_size)\n   \n\n    def rotate_image(self, image, angle):\n        image_center = tuple(np.array(image.shape[1::-1]) / 2)\n        rot_mat = cv2.getRotationMatrix2D(image_center, angle, 1.0)\n        result = cv2.warpAffine(image, rot_mat, image.shape[1::-1], flags=cv2.INTER_LINEAR)\n        return result\n    \n    def __getitem__(self,ids):\n        \n        id_path= self.paths[ids]\n        \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 =  [self.load_dicom_images_3d(x,split=\"train\") for x in batch_paths]\n            batch_X = np.stack(list_x, axis=0)\n            return batch_X,batch_y\n        else:\n            list_x =  self.load_dicom_images_3d(id_path,split=\"test\")\n            batch_X = np.stack(list_x)\n            return batch_X\n    \n    def load_dicom_images_3d(self, scan_id, num_imgs=config['frames_per_seq'], img_size=config['img_size'], \n                             mri_type=mri_types[0], split=\"train\", rotate=0):\n\n        target_file_paths = self.get_img_path_3d(scan_id, mri_type)\n        \n        img3d = np.array([self.read_mri(f) for f in target_file_paths]) # (12, 256, 256, 3)\n        \n        if img3d.shape[0] < num_imgs:\n            n_zero = np.zeros((num_imgs - img3d.shape[0],img_size, img_size, config['channels']))\n            img3d = np.concatenate((img3d,  n_zero), axis = 0)\n        \n        if np.min(img3d) < np.max(img3d):\n            img3d = img3d - np.min(img3d)\n            img3d = img3d / np.max(img3d)\n        \n        return img3d\n     \n    def crop_center_square(self, frame, scale=config['scale']):\n        y, x = frame.shape[0:2]\n        center_x, center_y = x / 2, y / 2\n        width_scaled, height_scaled = x * scale, y * 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        return frame[int(top_y):int(bottom_y), int(left_x):int(right_x)]\n    \n    def read_mri(self, path, voi_lut = True, fix_monochrome = True):\n        # Original from: https://www.kaggle.com/raddar/convert-dicom-to-np-array-the-correct-way\n        dicom = pydicom.read_file(path)\n        if voi_lut:\n            data = apply_voi_lut(dicom.pixel_array, dicom)\n        else:\n            data = dicom.pixel_array\n        if fix_monochrome and dicom.PhotometricInterpretation == \"MONOCHROME1\":\n            data = np.amax(data) - data\n        data = data - np.min(data)\n        data = data / np.max(data)\n        data = (data * 255).astype(np.uint8)\n        data = self.rotate_image(data, np.random.randint(0,20))\n        data = self.crop_center_square(data)\n        data = cv2.resize(data, (config['img_size'], config['img_size']))\n        data = np.repeat(data[..., np.newaxis], 3, -1) # 256,256,3\n        return data\n\n    def get_img_path_3d(self, scan_id, mri_type):\n        modality_path = os.path.join(scan_id, mri_type)\n        total_img_num = len(glob.glob(f\"{modality_path}/*.dcm\"))\n        files = sorted(glob.glob(f\"{modality_path}/*.dcm\"), \n                       key=lambda var:[int(x) if x.isdigit() else x for x in re.findall(r'[^0-9]|[0-9]+', var)])\n        mid_num = total_img_num // 2\n        num_3d2 = config['frames_per_seq'] // 2\n        start_idx = max(0, mid_num - num_3d2)\n        end_idx = min(len(files), mid_num + num_3d2)\n        target_file_paths = files[start_idx:end_idx]\n        return target_file_paths\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":{"execution":{"iopub.status.busy":"2021-09-27T17:38:55.952407Z","iopub.execute_input":"2021-09-27T17:38:55.952887Z","iopub.status.idle":"2021-09-27T17:38:55.976808Z","shell.execute_reply.started":"2021-09-27T17:38:55.952855Z","shell.execute_reply":"2021-09-27T17:38:55.976109Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_dataset = Dataset(df_data,batch_size=config['batch_size'])\n\nfor i in range(2):\n    images, label = train_dataset[i]\n    print(\"Dimension of the CT scan is:\", images.shape)\n    print(\"label=\",label.shape)\n    plt.imshow(images[0,0,:,:,0], cmap=\"gray\")\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2021-09-27T17:38:56.143539Z","iopub.execute_input":"2021-09-27T17:38:56.144086Z","iopub.status.idle":"2021-09-27T17:38:57.060065Z","shell.execute_reply.started":"2021-09-27T17:38:56.144049Z","shell.execute_reply":"2021-09-27T17:38:57.059264Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_3d_model(width=config['img_size'], height=config['img_size'], depth=config['frames_per_seq'], model_arch='custom'):\n    \"\"\"Build a 3D convolutional neural network model.\"\"\"\n    inputs = tf.keras.Input((width, height, depth, config['channels']))\n    if model_arch == \"custom\":\n        x = Conv3D(filters=64, kernel_size=3, padding='same', activation=\"relu\")(inputs)\n        x = MaxPool3D(pool_size=2)(x)\n        x = BatchNormalization()(x)\n    \n        x = Conv3D(filters=64, kernel_size=3, padding='same', activation=\"relu\")(inputs)\n        x = MaxPool3D(pool_size=2)(x)\n        x = BatchNormalization()(x)\n    \n        x = Conv3D(filters=128, kernel_size=3, padding='same', activation=\"relu\")(inputs)\n        x = MaxPool3D(pool_size=2)(x)\n        x = BatchNormalization()(x)\n    \n        x = Conv3D(filters=256, kernel_size=3, padding='same', activation=\"relu\")(x)\n        x = MaxPool3D(pool_size=2)(x)\n        x = BatchNormalization()(x)\n\n        x = GlobalAveragePooling3D()(x)\n        x = Dense(units=512, activation=\"relu\")(x)\n        x = Dropout(0.08)(x)\n\n        outputs = Dense(units=1, activation=\"sigmoid\")(x)\n        model = tf.keras.Model(inputs, outputs, name=\"3dcnn\")\n    else:\n        input_shape = (depth, width,height)\n        inputs = tf.keras.layers.Input((*input_shape,3), name='inputs')\n        x = Conv3D(filters=3, kernel_size = 3, strides=(1, 1, 1), padding='same', use_bias=True)(inputs)\n        net, _ = Classifiers.get(model_arch)\n        x = net(input_shape=(*input_shape,3),include_top=False, weights='imagenet')(x)\n        x = GlobalAveragePooling3D()(x)\n        x = Dropout(rate=0.5)(x)\n        outputs = Dense(1, activation='sigmoid', dtype='float32')(x)\n        model  = tf.keras.Model(inputs, outputs, name=model_arch)\n    return model\n\nmodel = get_3d_model(model_arch='seresnet50')\nmodel.summary()\n","metadata":{"execution":{"iopub.status.busy":"2021-09-27T17:38:57.062775Z","iopub.execute_input":"2021-09-27T17:38:57.063528Z","iopub.status.idle":"2021-09-27T17:38:59.920761Z","shell.execute_reply.started":"2021-09-27T17:38:57.063484Z","shell.execute_reply":"2021-09-27T17:38:59.91989Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot(history):\n    fig, ax = plt.subplots(1, 3, figsize=(20, 7))\n    ax = ax.ravel()\n    for fold in history:\n        for i, metric in enumerate([\"accuracy\",\"loss\",\"auc\"]):\n            ax[i].plot(history[fold].history[metric], label=\"train \"+str(fold))\n            ax[i].plot(history[fold].history[\"val_\" + metric], linestyle=\"dotted\", label=\"val \"+str(fold))\n            ax[i].set_title(\"Model {}\".format(metric))\n            ax[i].set_xlabel(\"epochs\")\n            ax[i].set_ylabel(metric)\n            ax[i].legend()\n    \ndef train_each_mri_type(mri_types, model_arch):\n    history = {}\n    for m_type in mri_types:\n        model = get_3d_model(model_arch=model_arch)\n        print(f\"Training for {m_type}\")\n        print('*'*100)\n        train_dataset = Dataset(df_data,batch_size=config['batch_size'])\n        valid_dataset = Dataset(df_data,batch_size=config['batch_size'])\n        optimizer = tf.keras.optimizers.Adam(learning_rate=config['learning_rate'])\n        early_stopping = tf.keras.callbacks.EarlyStopping(monitor='val_loss', patience=5, verbose=0, mode='min', restore_best_weights=True)\n        model_checkpoint = tf.keras.callbacks.ModelCheckpoint(f'{model_arch}_{m_type}.h5', save_best_only=True, save_weights_only=False)\n        LR = ReduceLROnPlateau(monitor='val_loss', factor=0.1, patience=5, min_lr=0.000001, verbose=1, mode='min')\n        model.compile(optimizer=optimizer,loss='binary_crossentropy', metrics=['accuracy', tf.keras.metrics.AUC(name='auc')])\n        log =  model.fit(train_dataset,validation_data=valid_dataset,epochs=config['num_epochs'],\n                         shuffle=True, callbacks=[LR, early_stopping, model_checkpoint])\n        history[m_type] = log\n    plot(history)\n    \ntrain_each_mri_type(mri_types, \"seresnet50\")","metadata":{"execution":{"iopub.status.busy":"2021-09-27T17:38:59.922936Z","iopub.execute_input":"2021-09-27T17:38:59.923319Z","iopub.status.idle":"2021-09-27T17:40:45.223744Z","shell.execute_reply.started":"2021-09-27T17:38:59.923279Z","shell.execute_reply":"2021-09-27T17:40:45.22263Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}