{"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":"# Genetic Biomarker prediction with CRNN - train\nA lot of this work is inspired from work shared by fellow kagglers, thank everyone for teaching me good stuff !!\n\n### Credits\n- https://keras.io/examples/vision/video_classification/\n- https://www.kaggle.com/ayuraj/train-brain-tumor-as-video-classification-w-b\n- https://www.kaggle.com/sergeybulanov/tf-simple-prediction-with-vgg16#2.-DataLoader\n- https://www.kaggle.com/michaelfumery/brain-tumor-transfert-learning-flair-kfold","metadata":{"id":"YhgooTH2rrEF"}},{"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\n\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 tensorflow.keras import backend as K, optimizers, regularizers\n\nimport wandb\nfrom wandb.keras import WandbCallback\n\nrn.seed(30)\nnp.random.seed(30)\ntf.compat.v1.random.set_random_seed(30)\nprint('W&B version: ', wandb.__version__)\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut","metadata":{"execution":{"iopub.status.busy":"2021-09-20T10:02:11.35705Z","iopub.execute_input":"2021-09-20T10:02:11.357317Z","iopub.status.idle":"2021-09-20T10:02:11.52941Z","shell.execute_reply.started":"2021-09-20T10:02:11.357292Z","shell.execute_reply":"2021-09-20T10:02:11.527875Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Min Images for MRI Type\n- T2w = 15\n- FLAIR = 14\n- T1w = 14\n- T1wCE = 15","metadata":{}},{"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  'output_path': './crnn/',\n  'nfolds': 3,\n  'global_seed': 42,\n  'batch_size': 4,\n  'frames_per_seq': 24,\n  'img_size': 224,\n  'learning_rate': 0.0001,\n  'rnn_cells': 16,  \n  'num_epochs': 10,\n  'channels': 3,\n  'scale' : 0.75\n}\n\nmri_types = ['T2w'] \n# mri_types = ['FLAIR','T1w','T1wCE','T2w']","metadata":{"id":"5EytvHEo4aAv","outputId":"04be0bdb-8991-45ba-bed4-7d6b55047394","execution":{"iopub.status.busy":"2021-09-20T10:02:11.536589Z","iopub.execute_input":"2021-09-20T10:02:11.540624Z","iopub.status.idle":"2021-09-20T10:02:11.553276Z","shell.execute_reply.started":"2021-09-20T10:02:11.540587Z","shell.execute_reply":"2021-09-20T10:02:11.547765Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data Preprocessing","metadata":{}},{"cell_type":"code","source":"class BrainTumor_GeneticSequence():\n    \"\"\"Prepares the train and the validation data pipeline for mri_type, for ex: mri_type = FLAIR\"\"\"\n    mri_type = \"FLAIR\"\n    df_data = None\n    df_train_labels = pd.read_csv(config['csv_path'])\n    \n    def __init__(self, mri_type):\n        self.mri_type = mri_type\n        self.df_data = pd.DataFrame(columns=['BraTS21ID'] + mri_types)\n        for key in mri_types:\n            self.df_data[key] = self.df_data[key].astype(int)\n        self.df_data['BraTS21ID'] = self.df_data['BraTS21ID'].astype(int)\n\n    def prepare_dataframe(self, mode='train'):\n        train_folders = ''\n        if mode == 'test':\n            folders_path = \"test_images_source_path\"\n        else:\n            folders_path = \"images_source_path\"\n        train_folders = config[folders_path] + '/'\n        for f in tqdm(os.listdir(train_folders)):\n            if f in [\"00109\", \"00123\", \"00709\"]: \n                continue\n            BraTS21ID = int(f)\n            self.df_data = self.df_data.append({'BraTS21ID': BraTS21ID, 'FLAIR': 0, 'T1w': 0, 'T1wCE': 0, 'T2w' : 0}, ignore_index=True)\n            BraTS21ID_key_path = f'{config[folders_path]}/{format(BraTS21ID, \"05d\")}/{self.mri_type}/*.dcm'\n            files_len = len(glob.glob(BraTS21ID_key_path))\n            # update file count or remove the patient from the dataset for the mri_type chosen\n            if files_len > 0:\n                self.df_data.loc[self.df_data['BraTS21ID'] == BraTS21ID, self.mri_type] = files_len\n            else:\n                self.df_data = self.df_data.loc[self.df_data.BraTS21ID!=BraTS21ID]\n        self.df_data[\"folder_name\"] = [format(x, '05d') for x in self.df_data[\"BraTS21ID\"]]\n        self.df_data[\"folder_path\"] = [os.path.join(config[folders_path], x) for x in self.df_data[\"folder_name\"]]\n#         self.df_data = self.df_data.head(30) # for testing\n#         print(self.df_data.head())\n        self.df_data = pd.merge(self.df_data, self.df_train_labels,how='left',on='BraTS21ID')\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 normalize(self, image):\n        result = cv2.normalize(image, None, alpha=0, beta=1, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_32F)\n        return result\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.normalize(data)\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#         (thresh, im_bw) = cv2.threshold(data, 128, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)\n#         data = cv2.threshold(data, thresh, 255, cv2.THRESH_BINARY)[1]\n        data = np.repeat(data[..., np.newaxis], 3, -1)\n        return data\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 get_img_path_3d(self, dir_path):\n        modality_path = os.path.join(dir_path.decode('utf8'), self.mri_type)\n        files = sorted(glob.glob(f\"{modality_path}/*.dcm\"), key=lambda var:[int(x) if x.isdigit() else x for x in re.findall(r'[^0-9]|[0-9]+', var)])\n        total_img_num = len(files) \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 = tf.convert_to_tensor(files[start_idx:end_idx], dtype=tf.string) \n        \n        def get_frames(path):\n            file_path = path.numpy().decode('UTF-8')\n            image = self.read_mri(file_path)    \n            return image\n    \n        mri_images = tf.nest.map_structure(tf.stop_gradient, tf.map_fn(fn=get_frames, elems=target_file_paths, fn_output_signature=tf.float32))\n\n        # padding null images \n        if mri_images.shape[0] < config['frames_per_seq']:\n            n_zero = tf.zeros((config['frames_per_seq'] - mri_images.shape[0], config['img_size'], config['img_size'], config['channels']), dtype=tf.dtypes.float32)\n            mri_images = np.concatenate((mri_images,  n_zero), axis = 0)\n        return mri_images\n    \n    def load_frame(self, df_dict):\n        dirname = df_dict['folder_path']\n        paths = tf.numpy_function(self.get_img_path_3d, [dirname], tf.float32)\n        label = df_dict['MGMT_value']\n        label = tf.cast(label, tf.float32)\n        return paths, label\n    \n    \n    def train(self, epochs=config['num_epochs'], output_path = config['output_path'], cnn_model_arch='custom', wandblogging = False):\n        tf.keras.backend.clear_session() \n        self.prepare_dataframe()\n        skf = StratifiedKFold(n_splits=config['nfolds'], shuffle=True)\n        fold = 1\n        history = {}\n        run =  None\n        \n        if wandblogging:\n            wandb.login()    \n            run = wandb.init(project='brain-tumor-video', job_type='dataloader-viz')\n            os.makedirs('gifs/', exist_ok=True)\n        \n        for train_index, valid_index in skf.split(self.df_data, self.df_data.MGMT_value.values):\n    \n            train_df = self.df_data.loc[train_index,:]\n            valid_df = self.df_data.loc[valid_index,:]\n            print(f'Size of train_df: {len(train_df)}; valid_df: {len(valid_df)}')\n\n            AUTOTUNE = tf.data.AUTOTUNE\n            \n            trainloader = tf.data.Dataset.from_tensor_slices(dict(train_df))\n            validloader = tf.data.Dataset.from_tensor_slices(dict(valid_df))\n            \n            trainloader = (\n                trainloader\n                .shuffle(1024)\n                .map(self.load_frame, num_parallel_calls=AUTOTUNE)\n                .batch(config['batch_size'])\n                .prefetch(AUTOTUNE)\n            )\n\n            validloader = (\n                validloader\n                .map(self.load_frame, num_parallel_calls=AUTOTUNE)\n                .batch(config['batch_size'])\n                .prefetch(AUTOTUNE)\n            )\n            \n            if wandblogging:\n                frames, labels = next(iter(trainloader))\n                for i, frame in enumerate(frames):\n                    imageio.mimsave(f'gifs/out_{i}.gif', (frame).numpy().astype('uint8')) \n                    wandb.log({'examples': [wandb.Image(f'gifs/out_{i}.gif', caption=f'{label.numpy()}') for i, label in enumerate(labels)]})\n        \n            model = create_custom_crnn_model(cnn_model_arch)\n#             model.summary()\n            optimizer = tf.keras.optimizers.Adam(learning_rate=config['learning_rate'])\n            loss_fn = tf.keras.losses.BinaryCrossentropy()\n            early_stopping = tf.keras.callbacks.EarlyStopping(monitor='val_accuracy', patience=5, verbose=0, mode='max', restore_best_weights=True)\n            model_checkpoint = tf.keras.callbacks.ModelCheckpoint(f'{cnn_model_arch}_{fold}.h5', save_best_only=True, save_weights_only=False)\n            LR = ReduceLROnPlateau(monitor='val_accuracy', factor=0.1, patience=5, min_lr=0.000001, verbose=1, mode='max')\n            model.compile(optimizer=optimizer,loss=loss_fn, metrics=['accuracy', tf.keras.metrics.AUC(name='auc')])\n            print(f'Training for fold...{fold}')\n            history[fold] = model.fit(trainloader, epochs=epochs,validation_data=validloader, callbacks=[LR, early_stopping, model_checkpoint])\n            fold+=1\n        \n        if wandblogging:\n            run.finish()\n        return history\n    \n    def predit(self, model):\n        AUTOTUNE = tf.data.AUTOTUNE\n        testloader = tf.data.Dataset.from_tensor_slices(dict(self.df_data))\n        testloader = (testloader\n                    .map(self.load_frame, num_parallel_calls=AUTOTUNE)\n                    .batch(config['batch_size'])\n                    .prefetch(AUTOTUNE)\n                    )\n        proba = model.predict(testloader, verbose=1)\n        return proba ","metadata":{"execution":{"iopub.status.busy":"2021-09-20T10:02:11.557604Z","iopub.execute_input":"2021-09-20T10:02:11.558609Z","iopub.status.idle":"2021-09-20T10:02:11.646643Z","shell.execute_reply.started":"2021-09-20T10:02:11.558573Z","shell.execute_reply":"2021-09-20T10:02:11.645843Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Test Transformations\nHere let us see difference between original image and transformed image","metadata":{}},{"cell_type":"code","source":"dp = BrainTumor_GeneticSequence('FLAIR')\nsample_img = pydicom.read_file(f'{config[\"images_source_path\"]}/00068/FLAIR/Image-132.dcm').pixel_array\nprocessed_img = dp.read_mri(f'{config[\"images_source_path\"]}/00068/FLAIR/Image-132.dcm')\n\nfig = plt.figure(figsize=(12,8))\nax1 = plt.subplot(1,2,1)\nax1.imshow(sample_img, cmap=\"gray\")\nax1.set_title(f\"Original image shape = {sample_img.shape}\")\n\nax2 = plt.subplot(1,2,2)\nax2.imshow(processed_img[:,:,0], cmap=\"gray\")\nax2.set_title(f\"Preproc image shape = {processed_img.shape}\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-09-20T10:02:11.650527Z","iopub.execute_input":"2021-09-20T10:02:11.650871Z","iopub.status.idle":"2021-09-20T10:02:12.311884Z","shell.execute_reply.started":"2021-09-20T10:02:11.650835Z","shell.execute_reply":"2021-09-20T10:02:12.311196Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from tensorflow.keras.applications import *\n\ndef FeatureExtractor(model_arch='custom'):\n    shape=(config['img_size'],config['img_size'],config['channels'])\n    model = Sequential()\n    if model_arch == 'custom':\n        model.add(Conv2D(64, (3, 3), padding='same', activation='relu', input_shape=shape))\n        model.add(MaxPooling2D(pool_size=(2, 2)))\n        model.add(BatchNormalization())\n\n        model.add(Conv2D(64, (3, 3), padding='same', activation='relu',))\n        model.add(MaxPooling2D(pool_size=(2, 2)))\n        model.add(BatchNormalization())\n\n        model.add(Conv2D(128, (3, 3), padding='same', activation='relu'))\n        model.add(MaxPooling2D(pool_size=(2, 2)))\n        model.add(BatchNormalization())\n        \n        model.add(Conv2D(128, (3, 3), padding='same', activation='relu'))\n        model.add(MaxPooling2D(pool_size=(2, 2)))\n\n        model.add(Flatten())\n        model.add(Dense(512, activation='relu'))\n        model.add(Dense(256))\n    elif model_arch == 'resnet50':\n        model = ResNet50(include_top=False, weights='imagenet', input_shape=shape)\n        for layer in model.layers[:-10]:\n            layer.trainable=True\n    elif model_arch == 'efb7':\n        model = EfficientNetB7(include_top=False,weights=\"imagenet\", input_shape=shape)\n        for layer in model.layers[:-10]:\n            layer.trainable=True\n    elif model_arch == 'xception':\n        model = Xception(include_top=False,weights=\"imagenet\", input_shape=shape)\n        for layer in model.layers[:-10]:\n            layer.trainable=True\n    return model\n\ndef create_custom_crnn_model(cnn_model_arch='custom'):\n    feature_extractor = FeatureExtractor(cnn_model_arch)\n    model = Sequential()\n    model.add(TimeDistributed(feature_extractor,input_shape=(config['frames_per_seq'], config['img_size'], config['img_size'], config['channels'])))\n    model.add(TimeDistributed(Flatten()))\n    model.add(GRU(config['rnn_cells'],return_sequences=True)) # using GRU here because they have less parameters compared to LSTM\n    model.add(GRU(config['rnn_cells']//2))\n    model.add(Dense(256, activation='relu', kernel_regularizer=regularizers.l2(0.01)))\n    model.add(Dropout(0.4))\n    model.add(Dense(1,activation='sigmoid'))\n    return model\n\ndef 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, cnn_model_arch):\n    for m_type in mri_types:\n        dp = BrainTumor_GeneticSequence(m_type)\n        print(\"*\"*100)\n        print(f\"Training for {m_type} with {cnn_model_arch} as feature extractor, configuration: {config}\")\n        print('*'*100)\n        history = dp.train(epochs=10,cnn_model_arch=cnn_model_arch)\n        plot(history)","metadata":{"id":"7wGy5yx9fTps","outputId":"a0d1e04c-06aa-43ff-e47b-c6319ad80373","execution":{"iopub.status.busy":"2021-09-20T10:02:12.313194Z","iopub.execute_input":"2021-09-20T10:02:12.314053Z","iopub.status.idle":"2021-09-20T10:02:12.339495Z","shell.execute_reply.started":"2021-09-20T10:02:12.314017Z","shell.execute_reply":"2021-09-20T10:02:12.33744Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Training","metadata":{}},{"cell_type":"code","source":"%%time\nmodel_arch = 'xception'\ntrain_each_mri_type(mri_types, model_arch)","metadata":{"execution":{"iopub.status.busy":"2021-09-20T10:02:12.34557Z","iopub.execute_input":"2021-09-20T10:02:12.34648Z","iopub.status.idle":"2021-09-20T10:03:27.585309Z","shell.execute_reply.started":"2021-09-20T10:02:12.346445Z","shell.execute_reply":"2021-09-20T10:03:27.584631Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Prediction","metadata":{}},{"cell_type":"code","source":"model = tf.keras.models.load_model(f'./{model_arch}_1.h5')\ndp = BrainTumor_GeneticSequence(mri_types[0])\ndp.prepare_dataframe(mode='test')\ndp.df_data['MGMT_pred'] = 0\nproba = dp.predit(model)\nsample_submission_path = os.path.join(config['data_path'], 'sample_submission.csv')\nsample_df = pd.read_csv(sample_submission_path); \ndp.df_data['MGMT_pred'] += proba.squeeze()    \ndp.df_data[['BraTS21ID','MGMT_pred']].head()","metadata":{"execution":{"iopub.status.busy":"2021-09-20T10:03:27.586548Z","iopub.execute_input":"2021-09-20T10:03:27.586867Z","iopub.status.idle":"2021-09-20T10:03:27.614748Z","shell.execute_reply.started":"2021-09-20T10:03:27.58683Z","shell.execute_reply":"2021-09-20T10:03:27.613592Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}