{"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 Tumor Segmentation BRATS2021**\n\n**Problems:**\n\n*1.- Due to memory comsumption the max. size of dataset are 20 nii files, and then I have checked the results are not quite good in predictions with rsna dataset. So, I guess that ensemble is not working quite well as well...*\n\n*2.- TPU are not working quite well, so I have used GPU for running the model, I you can give some advise about what is going on, I definetelly thank you!*\n\nLoading classes for segmentation from Brats 2020:\n\n-[1] https://www.kaggle.com/frlemarchand/brain-tumour-segmentation-in-mri-slices/\n\n-[2] https://www.kaggle.com/arashmehrzadi/brain-tumor-segmentation-unet/output\n\n\n***After testing the performance for segmentation the rsna data with [2], the accuracy is not working very well, so, we are trying to run the U-Net with Brats 2021 data, and finally doing the ensemble with 2020, and check again:\n\n- https://www.kaggle.com/dschettler8845/how-to-load-basic-data-exploration","metadata":{}},{"cell_type":"code","source":"###Libraries and imports\nimport numpy as np\nimport math\nimport random\nimport os\nimport shutil\nimport gzip\nimport shutil\nimport glob\nimport gc\nimport cv2 as cv\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport tensorflow as tf\nimport numpy as np\nimport pandas as pd\nimport tarfile\nimport PIL\nimport scipy.misc\nimport skimage\nimport nibabel as nib\nimport matplotlib.pyplot as plt\n\nfrom PIL import Image\nfrom tqdm import tqdm\nfrom tensorflow import keras\nfrom keras.models import Model, load_model\nfrom keras.layers import Input ,BatchNormalization , Activation \nfrom keras.layers.convolutional import Conv2D, UpSampling2D\nfrom keras.layers.pooling import MaxPooling2D\nfrom keras.layers.merge import concatenate\nfrom keras.callbacks import EarlyStopping, ModelCheckpoint\nfrom keras import optimizers \nfrom keras import backend as K\nfrom sklearn.model_selection import train_test_split\nfrom skimage.morphology import ball, disk, dilation, binary_erosion, remove_small_objects, erosion, closing, reconstruction, binary_closing\nfrom skimage.measure import label,regionprops, perimeter\nfrom skimage.morphology import binary_dilation, binary_opening\nfrom skimage.filters import roberts, sobel\nfrom skimage import measure, feature\nfrom skimage.segmentation import clear_border\nfrom skimage import data\nfrom skimage.io import imread\nfrom scipy import ndimage as ndi\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\nfrom glob import glob\n\n\nDEVICE = \"GPU\"","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-08-25T10:02:18.419891Z","iopub.execute_input":"2021-08-25T10:02:18.420306Z","iopub.status.idle":"2021-08-25T10:02:28.572545Z","shell.execute_reply.started":"2021-08-25T10:02:18.420211Z","shell.execute_reply":"2021-08-25T10:02:28.571375Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if DEVICE == \"TPU\":\n    print(\"connecting to TPU...\")\n    try:\n        tpu = tf.distribute.cluster_resolver.TPUClusterResolver()\n        print('Running on TPU ', tpu.master())\n    except ValueError:\n        print(\"Could not connect to TPU\")\n        tpu = None\n\n    if tpu:\n        try:\n            print(\"initializing  TPU ...\")\n            tf.config.experimental_connect_to_cluster(tpu)\n            tf.tpu.experimental.initialize_tpu_system(tpu)\n            strategy = tf.distribute.experimental.TPUStrategy(tpu)\n            print(\"TPU initialized\")\n        except _:\n            print(\"failed to initialize TPU\")\n    else:\n        DEVICE = \"GPU\"\n\nif DEVICE != \"TPU\":\n    print(\"Using default strategy for CPU and single GPU\")\n    strategy = tf.distribute.get_strategy()\n\nif DEVICE == \"GPU\":\n    print(\"Num GPUs Available: \", len(tf.config.experimental.list_physical_devices('GPU')))\n    \n\nAUTO     = tf.data.experimental.AUTOTUNE\nREPLICAS = strategy.num_replicas_in_sync\nprint(f'REPLICAS: {REPLICAS}')","metadata":{"execution":{"iopub.status.busy":"2021-08-25T10:02:28.574402Z","iopub.execute_input":"2021-08-25T10:02:28.57485Z","iopub.status.idle":"2021-08-25T10:02:28.599843Z","shell.execute_reply.started":"2021-08-25T10:02:28.574804Z","shell.execute_reply":"2021-08-25T10:02:28.598466Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def Convolution(input_tensor,filters):\n    x = Conv2D(filters=filters,kernel_size=(3, 3),padding = 'same',strides=(1, 1))(input_tensor)\n    x = BatchNormalization()(x)\n    x = Activation('relu')(x) \n    return x\n\ndef model(input_shape):\n    \n    inputs = Input((input_shape))\n    \n    conv_1 = Convolution(inputs,32)\n    maxp_1 = MaxPooling2D(pool_size = (2, 2), strides = (2, 2), padding = 'same') (conv_1)\n    \n    conv_2 = Convolution(maxp_1,64)\n    maxp_2 = MaxPooling2D(pool_size = (2, 2), strides = (2, 2), padding = 'same') (conv_2)\n    \n    conv_3 = Convolution(maxp_2,128)\n    maxp_3 = MaxPooling2D(pool_size = (2, 2), strides = (2, 2), padding = 'same') (conv_3)\n    \n    conv_4 = Convolution(maxp_3,256)\n    maxp_4 = MaxPooling2D(pool_size = (2, 2), strides = (2, 2), padding = 'same') (conv_4)\n    \n    conv_5 = Convolution(maxp_4,512)\n    upsample_6 = UpSampling2D((2, 2)) (conv_5)\n    \n    conv_6 = Convolution(upsample_6,256)\n    upsample_7 = UpSampling2D((2, 2)) (conv_6)\n    \n    upsample_7 = concatenate([upsample_7, conv_3])\n    \n    conv_7 = Convolution(upsample_7,128)\n    upsample_8 = UpSampling2D((2, 2)) (conv_7)\n    \n    conv_8 = Convolution(upsample_8,64)\n    upsample_9 = UpSampling2D((2, 2)) (conv_8)\n    \n    upsample_9 = concatenate([upsample_9, conv_1])\n    \n    conv_9 = Convolution(upsample_9,32)\n    outputs = Conv2D(1, (1, 1), activation='sigmoid') (conv_9)\n    \n    model = Model(inputs=[inputs], outputs=[outputs]) \n    \n    return model\n\n# Computing Dice_Coefficient\ndef dice_coef(y_true, y_pred, smooth=1.0):\n    y_true_f = K.flatten(y_true)\n    y_pred_f = K.flatten(y_pred)\n    intersection = K.sum(y_true_f * y_pred_f)\n    return (2. * intersection + smooth) / (K.sum(y_true_f) + K.sum(y_pred_f) + smooth)\n\n# Computing Precision \ndef precision(y_true, y_pred):\n        true_positives = K.sum(K.round(K.clip(y_true * y_pred, 0, 1)))\n        predicted_positives = K.sum(K.round(K.clip(y_pred, 0, 1)))\n        precision = true_positives / (predicted_positives + K.epsilon())\n        return precision\n\n# Computing Sensitivity      \ndef sensitivity(y_true, y_pred):\n    true_positives = K.sum(K.round(K.clip(y_true * y_pred, 0, 1)))\n    possible_positives = K.sum(K.round(K.clip(y_true, 0, 1)))\n    return true_positives / (possible_positives + K.epsilon())\n\n# Computing Specificity\ndef specificity(y_true, y_pred):\n    true_negatives = K.sum(K.round(K.clip((1-y_true) * (1-y_pred), 0, 1)))\n    possible_negatives = K.sum(K.round(K.clip(1-y_true, 0, 1)))\n    return true_negatives / (possible_negatives + K.epsilon())\n\n# Accuracy vs Epoch\ndef Accuracy_Graph(history):\n    plt.plot(history.history['accuracy'])\n    plt.plot(history.history['val_accuracy'])\n    #plt.title('Model accuracy')\n    plt.ylabel('Accuracy')\n    plt.xlabel('Epoch')\n    plt.legend(['Train', 'Validation'], loc='upper left')\n    plt.subplots_adjust(top=1.00, bottom=0.0, left=0.0, right=0.95, hspace=0.25,\n                        wspace=0.35)\n    plt.show()\n    \n# Dice Similarity Coefficient vs Epoch\ndef Dice_coefficient_Graph(history):\n\n    plt.plot(history.history['dice_coef'])\n    plt.plot(history.history['val_dice_coef'])\n    #plt.title('Dice_Coefficient')\n    plt.ylabel('Dice_Coefficient')\n    plt.xlabel('Epoch')\n    plt.legend(['Train', 'Validation'], loc='upper left')\n    plt.subplots_adjust(top=1.00, bottom=0.0, left=0.0, right=0.95, hspace=0.25,\n                        wspace=0.35)\n    plt.show()\n# Loss vs Epoch\ndef Loss_Graph(history):\n\n    plt.plot(history.history['loss'])\n    plt.plot(history.history['val_loss'])\n    #plt.title('Model loss')\n    plt.ylabel('Loss')\n    plt.xlabel('Epoch')\n    plt.legend(['Train', 'Validation'], loc='upper left')\n    plt.subplots_adjust(top=1.00, bottom=0.0, left=0.0, right=0.95, hspace=0.25,\n                        wspace=0.35)\n    plt.show()\n\nprint(\"-- Compiling the model with: \", DEVICE)\ninput_size = 240\nif DEVICE=='TPU':\n    tpu = tf.distribute.cluster_resolver.TPUClusterResolver.connect()\n    tpu_strategy = tf.distribute.experimental.TPUStrategy(tpu)\n    with tpu_strategy.scope():\n        model = model(input_shape = (input_size, input_size, 1))\n        Adam=optimizers.Adam(lr=0.001)\n        model.compile(optimizer=Adam, loss='binary_crossentropy', metrics=['accuracy',dice_coef,precision,sensitivity,specificity])        \nelse:   \n    model = model(input_shape = (input_size, input_size, 1))\n    Adam=optimizers.Adam(lr=0.001)\n    model.compile(optimizer=Adam, loss='binary_crossentropy', metrics=['accuracy',dice_coef,precision,sensitivity,specificity])","metadata":{"execution":{"iopub.status.busy":"2021-08-25T10:02:28.602458Z","iopub.execute_input":"2021-08-25T10:02:28.602987Z","iopub.status.idle":"2021-08-25T10:02:29.006354Z","shell.execute_reply.started":"2021-08-25T10:02:28.602934Z","shell.execute_reply":"2021-08-25T10:02:29.00494Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!rm -Rf ./post_process_data\n!rm -Rf ./data\n!rm -Rf ./final_models","metadata":{"execution":{"iopub.status.busy":"2021-08-25T10:02:29.008344Z","iopub.execute_input":"2021-08-25T10:02:29.008839Z","iopub.status.idle":"2021-08-25T10:02:31.292057Z","shell.execute_reply.started":"2021-08-25T10:02:29.00879Z","shell.execute_reply":"2021-08-25T10:02:31.290773Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!mkdir data\n\ntar = tarfile.open(\"/kaggle/input/brats-2021-task1/BraTS2021_Training_Data.tar\")\ntar.extractall(\"./data\")\ntar.close()\n\ntar = tarfile.open(\"/kaggle/input/brats-2021-task1/BraTS2021_00495.tar\")\ntar.extractall(\"./data\")\ntar.close()\n\ntar = tarfile.open(\"/kaggle/input/brats-2021-task1/BraTS2021_00621.tar\")\ntar.extractall(\"./data\", )\ntar.close()\n\nimg_id = \"01281\"\nplt.figure(figsize=(18, 5))\nfor i, nii in enumerate([f'./data/BraTS2021_{img_id}/BraTS2021_{img_id}_{s_type}.nii.gz' for s_type in [\"flair\", \"t1\", \"t1ce\", \"t2\", \"seg\"]]):\n    plt.subplot(1,5,i+1)\n    image=nib.load(nii).get_fdata()\n    plt.title(nii.rsplit(\"_\", 1)[1].split(\".\", 1)[0], fontweight=\"bold\")\n    plt.axis(False)\n    plt.imshow(image[:, :, 80], cmap=\"bone\")\nplt.tight_layout()    \nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-25T10:03:59.038287Z","iopub.execute_input":"2021-08-25T10:03:59.038754Z","iopub.status.idle":"2021-08-25T10:06:48.836122Z","shell.execute_reply.started":"2021-08-25T10:03:59.038711Z","shell.execute_reply":"2021-08-25T10:06:48.83483Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!ls ./data/BraTS2021_00000/","metadata":{"execution":{"iopub.status.busy":"2021-08-25T10:06:48.837768Z","iopub.execute_input":"2021-08-25T10:06:48.838087Z","iopub.status.idle":"2021-08-25T10:06:49.603113Z","shell.execute_reply.started":"2021-08-25T10:06:48.838054Z","shell.execute_reply":"2021-08-25T10:06:49.601977Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Dataset organization in subforders**","metadata":{}},{"cell_type":"code","source":"import glob\nPath_brats= './data/BraTS2021_'\npost_data='./post_process_data'\nfinal_models='./final_models'\nos.mkdir(post_data)\nos.mkdir(final_models)\n\n!ls ./data/ | awk -F'[_]' '{print $2}' > test.txt\npost_data='./post_process_data/'\nwith open('./test.txt') as f:\n    for linea in f:\n        path_dest=post_data+linea\n        path_dest=path_dest.replace('\\n','')\n        if os.path.isdir(path_dest) != True:\n            os.mkdir(path_dest)\n        create_path_brats=Path_brats+linea+'/*'\n        for name in sorted(glob.glob(create_path_brats.replace('\\n',''))):        \n            shutil.move(name,path_dest)","metadata":{"execution":{"iopub.status.busy":"2021-08-25T10:06:49.605406Z","iopub.execute_input":"2021-08-25T10:06:49.605748Z","iopub.status.idle":"2021-08-25T10:06:54.004434Z","shell.execute_reply.started":"2021-08-25T10:06:49.60571Z","shell.execute_reply":"2021-08-25T10:06:54.003098Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Preparing dataset and model**","metadata":{}},{"cell_type":"code","source":"def Data_Concatenate(Input_Data):\n    counter=0\n    Output= []\n    for i in range(5):\n        print('$')\n        c=0\n        counter=0\n        for ii in range(len(Input_Data)):\n            if (counter != len(Input_Data)):\n                a= Input_Data[counter][:,:,:,i]\n                #print('a={}'.format(a.shape))\n                b= Input_Data[counter+1][:,:,:,i]\n                #print('b={}'.format(b.shape))\n                if(counter==0):\n                    c= np.concatenate((a, b), axis=0)\n                    print('c1={}'.format(c.shape))\n                    counter= counter+2\n                else:\n                    c1= np.concatenate((a, b), axis=0)\n                    c= np.concatenate((c, c1), axis=0)\n                    print('c2={}'.format(c.shape))\n                    counter= counter+2\n        c= c[:,:,:,np.newaxis]\n        Output.append(c)\n    return Output\n    \ndef launch_model(Input_Data,code,model):\n    \n    InData= Data_Concatenate(Input_Data)\n    AIO= concatenate(InData, axis=3)\n    AIO=np.array(AIO,dtype='float32')\n    TR=np.array(AIO[:,:,:,1],dtype='float32')\n    TRL=np.array(AIO[:,:,:,4],dtype='float32')\n    X_train , X_test, Y_train, Y_test = train_test_split(TR, TRL, test_size=0.15, random_state=32)\n    AIO=TRL=0\n\n    # Fitting the model over the data\n    print(\"-- Fitting the model over the data --\")\n    history = model.fit(X_train,Y_train,batch_size=32,epochs=20,validation_split=0.20,verbose=1,initial_epoch=0)\n    \n    # Evaluating the model on the training and testing data \n    print(\"-- Evaluating the model on the training and testing data --\")\n    model.evaluate(x=X_train, y=Y_train, batch_size=32 , verbose=1, sample_weight=None, steps=None)\n    model.evaluate(x=X_test, y=Y_test, batch_size=32, verbose=1, sample_weight=None, steps=None)    \n    \n    # Plotting the Graphs of Accuracy, Dice_coefficient, Loss at each epoch on Training and Testing data\n    print(\"-- Plotting the Graphs of Accuracy, Dice_coefficient, Loss at each epoch on Training and Testing data --\")\n    Accuracy_Graph(history)\n    Dice_coefficient_Graph(history)\n    Loss_Graph(history)\n               \n    model.save('./final_models/BraTs2021_'+code+'.h5')","metadata":{"execution":{"iopub.status.busy":"2021-08-25T10:06:54.007585Z","iopub.execute_input":"2021-08-25T10:06:54.007935Z","iopub.status.idle":"2021-08-25T10:06:54.031702Z","shell.execute_reply.started":"2021-08-25T10:06:54.0079Z","shell.execute_reply":"2021-08-25T10:06:54.029309Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Launch model**\n\n* top_limit_number = Number of loops that we take.\n* split_number = Number of elements that we take for the look\n* init_counter = Loop local counter\n* inside_split_countert = Loop global counter\n","metadata":{}},{"cell_type":"code","source":"import glob\nPath= './post_process_data'\np=os.listdir(Path)\nInput_Data= []\n\n\ndef Data_Preprocessing(modalities_dir):\n    all_modalities = []    \n    for modality in modalities_dir:      \n        nifti_file   = nib.load(modality)\n        brain_numpy  = np.asarray(nifti_file.dataobj)    \n        all_modalities.append(brain_numpy)\n    brain_affine   = nifti_file.affine\n    all_modalities = np.array(all_modalities)\n    all_modalities = np.rint(all_modalities).astype(np.int16)\n    all_modalities = all_modalities[:, :, :, :]\n    all_modalities = np.transpose(all_modalities)\n    return all_modalities\n\ntop_limit_number = 15\nsplit_number = 20\ninit_counter=0\ninside_split_countert=1\ntotal_count_img, partial_count_img = len(p), int(len(p)/split_number)\n\nfor i in tqdm(p):\n    if (int(init_counter*inside_split_countert) == total_count_img) or (int(top_limit_number*split_number)==int(init_counter*inside_split_countert)):\n        print(\"Launch Final model.\")\n        launch_model(Input_Data,str(init_counter), model)\n        del(Input_Data)\n        gc.collect()\n        break    \n    if (init_counter == split_number):\n        print(\"Launch model :\"+str(init_counter*inside_split_countert))\n        launch_model(Input_Data,str(init_counter*inside_split_countert),model)\n        del(Input_Data)\n        gc.collect()\n        Input_Data= []\n        inside_split_countert=inside_split_countert+1\n        init_counter=0\n    create_path_post=Path+'/'+i+'/*'\n    \n    for name in sorted(glob.glob(create_path_post.replace('\\n',''))):\n        os.system('gunzip ' + name)\n    brain_dir = os.path.normpath(Path+'/'+i+'/')\n    flair     = glob.glob(os.path.join(brain_dir, '*_flair*.nii'))\n    t1        = glob.glob(os.path.join(brain_dir, '*_t1*.nii'))\n    t1ce      = glob.glob(os.path.join(brain_dir, '*_t1ce*.nii'))\n    t2        = glob.glob(os.path.join(brain_dir, '*_t2*.nii'))\n    gt        = glob.glob( os.path.join(brain_dir, '*_seg*.nii'))\n    modalities_dir = [flair[0], t1[0], t1ce[0], t2[0], gt[0]]\n    P_Data = Data_Preprocessing(modalities_dir)\n    Input_Data.append(P_Data)\n    shutil.rmtree(Path+'/'+i)\n    init_counter = init_counter + 1\n     ","metadata":{"execution":{"iopub.status.busy":"2021-08-25T10:06:54.033509Z","iopub.execute_input":"2021-08-25T10:06:54.033987Z","iopub.status.idle":"2021-08-25T10:13:11.730817Z","shell.execute_reply.started":"2021-08-25T10:06:54.03395Z","shell.execute_reply":"2021-08-25T10:13:11.729731Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!rm -Rf ./post_process_data\n!rm -Rf ./data","metadata":{},"execution_count":null,"outputs":[]}]}