{"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 numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os\nimport gc\nimport time\nfrom IPython.display import clear_output\nfrom tensorflow.keras.models import load_model\nfrom tensorflow.keras.callbacks import ModelCheckpoint as MC\nfrom tensorflow.keras import backend as K\nfrom sklearn.metrics import confusion_matrix, roc_curve, auc, recall_score, accuracy_score, balanced_accuracy_score, precision_score\nimport pickle\nfrom IPython.display import FileLink\n\n\n\nroot = '/kaggle/input/rsna-str-pulmonary-embolism-detection'\nfor item in os.listdir(root):\n    path = os.path.join(root, item)\n    if os.path.isfile(path):\n        print(path)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-01-31T22:15:52.648638Z","iopub.execute_input":"2022-01-31T22:15:52.649436Z","iopub.status.idle":"2022-01-31T22:15:58.345857Z","shell.execute_reply.started":"2022-01-31T22:15:52.64932Z","shell.execute_reply":"2022-01-31T22:15:58.3447Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load Data\n\nNext, we load all '.csv' files into memory and peek into their makeup.","metadata":{}},{"cell_type":"code","source":"print('Reading train data...')\ndata = pd.read_csv(\"../input/rsna-str-pulmonary-embolism-detection/train.csv\")\nprint(data.shape)\ndata.head()","metadata":{"execution":{"iopub.status.busy":"2022-01-31T22:16:07.606256Z","iopub.execute_input":"2022-01-31T22:16:07.607056Z","iopub.status.idle":"2022-01-31T22:16:11.189643Z","shell.execute_reply.started":"2022-01-31T22:16:07.607012Z","shell.execute_reply":"2022-01-31T22:16:11.188302Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Negative exam selection\n\nOnly the positive exams containing pulmonary embolism have information on rv_lv_ratio_gte_1 and rv_lv_ratio_lt_1. So next we select only the images from the positive exams.","metadata":{}},{"cell_type":"code","source":"#Selecting only positive exams\ndata = data.loc[data['negative_exam_for_pe'] == 0] \n\n# Agrupar por col_groupby \ncol_index = 'SOPInstanceUID'\ncol_groupby = 'StudyInstanceUID'\ndata = data[data[col_groupby].duplicated()==False].reset_index(drop=True)\ndata","metadata":{"execution":{"iopub.status.busy":"2022-01-31T22:16:20.884708Z","iopub.execute_input":"2022-01-31T22:16:20.885244Z","iopub.status.idle":"2022-01-31T22:16:21.007594Z","shell.execute_reply.started":"2022-01-31T22:16:20.885204Z","shell.execute_reply":"2022-01-31T22:16:21.000766Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#For testing a smaller dataset\ndata_small = data[:250]\ndata_small.shape","metadata":{"execution":{"iopub.status.busy":"2022-01-31T21:06:11.140315Z","iopub.execute_input":"2022-01-31T21:06:11.140801Z","iopub.status.idle":"2022-01-31T21:06:11.147046Z","shell.execute_reply.started":"2022-01-31T21:06:11.140765Z","shell.execute_reply":"2022-01-31T21:06:11.146287Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train/Test Split\n\nSince this project is for academic purposes, we will split the given train set into train+test set, to evaluate performance of the model better. Since the data is already structured by exam, a simple train_test_split(80% train, 20% test) is employed.","metadata":{}},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split\ntrain,test=train_test_split(data, test_size=0.2, random_state=42, shuffle=False)\ntrain","metadata":{"execution":{"iopub.status.busy":"2022-01-31T22:16:31.767962Z","iopub.execute_input":"2022-01-31T22:16:31.768249Z","iopub.status.idle":"2022-01-31T22:16:31.807006Z","shell.execute_reply.started":"2022-01-31T22:16:31.768217Z","shell.execute_reply":"2022-01-31T22:16:31.806275Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data pre-processing functions\n\n- normalize_image: Converts the CT scan data from Hounsfield scale to a 0-1 scale.\n- resize_volume: Resizes the volume of each exam to lower computational costs and ensure training uniformity. We resized to 128x128 on the slice level and 64 on exam depth (number of slices).\n\nThese functions were adapted from '3D image classification from CT scans' example in the Keras package documentation","metadata":{}},{"cell_type":"code","source":"def normalize_image(image):\n    min = -1000\n    max = 400\n    image[image < min] = min\n    image[image > max] = max\n    image = (image - min) / (max - min)\n    image = image.astype(\"float32\")\n    return image","metadata":{"execution":{"iopub.status.busy":"2022-01-31T22:16:50.014877Z","iopub.execute_input":"2022-01-31T22:16:50.015468Z","iopub.status.idle":"2022-01-31T22:16:50.02099Z","shell.execute_reply.started":"2022-01-31T22:16:50.015422Z","shell.execute_reply":"2022-01-31T22:16:50.020294Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from scipy import ndimage\ndef resize_volume(img):\n    \"\"\"Resize across z-axis\"\"\"\n    # Set the desired depth\n    desired_depth = 64\n    desired_width = 128\n    desired_height = 128\n    # Get current depth\n    current_depth = img.shape[-1]\n    current_width = img.shape[0]\n    current_height = img.shape[1]\n    # Compute depth factor\n    depth = current_depth / desired_depth\n    width = current_width / desired_width\n    height = current_height / desired_height\n    depth_factor = 1 / depth\n    width_factor = 1 / width\n    height_factor = 1 / height\n    # Rotate\n    #img = ndimage.rotate(img, 90, reshape=False)\n    # Resize across z-axis\n    img = ndimage.zoom(img, (width_factor, height_factor, depth_factor), order=1)\n    return img","metadata":{"execution":{"iopub.status.busy":"2022-01-31T22:16:52.174534Z","iopub.execute_input":"2022-01-31T22:16:52.175142Z","iopub.status.idle":"2022-01-31T22:16:52.181289Z","shell.execute_reply.started":"2022-01-31T22:16:52.175105Z","shell.execute_reply":"2022-01-31T22:16:52.180608Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Exam reader\n\nThe function get_exam gets the image arrays from a Dicom image and groups them by exam. This function was adapted from a function developed by [eladwar](https://www.kaggle.com/eladwar) in this notebook [here](https://www.kaggle.com/eladwar/20-seconds-or-less#VTK-is-mostly-written-in-C++-making-it-incredibly-efficient.-By-using-this-library-you-can-save-loads-of-memory-and-time.).","metadata":{}},{"cell_type":"code","source":"import vtk\nfrom vtk.util import numpy_support\nimport cv2\n\ndef get_exam(path):\n    x = 0\n    #scan=np.array([])\n    n_slices = len([name for name in os.listdir(directory)])  #number of slices in exam = number of images in exam file\n    exam = np.zeros((512, 512, n_slices))\n    for file in os.listdir(directory):\n        f = os.path.join(directory, file)\n\n        reader = vtk.vtkDICOMImageReader()\n        reader.SetFileName(f)\n        reader.Update()\n        _extent = reader.GetDataExtent()\n        ConstPixelDims = [_extent[1]-_extent[0]+1, _extent[3]-_extent[2]+1, _extent[5]-_extent[4]+1]\n\n        ConstPixelSpacing = reader.GetPixelSpacing()\n        imageData = reader.GetOutput()\n        pointData = imageData.GetPointData()\n        arrayData = pointData.GetArray(0)\n        ArrayDicom = numpy_support.vtk_to_numpy(arrayData)\n        ArrayDicom = ArrayDicom.reshape(ConstPixelDims, order='F')\n        ArrayDicom = cv2.resize(ArrayDicom,(512,512))\n        #ArrayDicom = ArrayDicom.reshape((1,ArrayDicom.shape[0],ArrayDicom.shape[1]))  #For np.concatenate method\n\n        ArrayDicom = normalize_image(ArrayDicom)\n\n\n        exam[:, :, x] = ArrayDicom  #Simply allocating each slice to the vector is more efficient than concatenating\n\n        #if scan.size==0:\n            #scan = ArrayDicom\n        #else:\n            #scan = np.concatenate((scan, ArrayDicom))\n        x+=1\n    exam[:, :, 1].shape\n    exam = resize_volume(exam)\n    \n    exam = exam.reshape((exam.shape[0],exam.shape[1],exam.shape[2],1))\n    \n    return exam\n\n\n#Check exam shape (width, height, depth)\ndirectory = os.fsencode(\"../input/rsna-str-pulmonary-embolism-detection/train/017320409cc6/fbd7d4bb6cb5\")\nexam = get_exam(directory)\nexam.shape","metadata":{"execution":{"iopub.status.busy":"2022-01-31T22:16:54.801901Z","iopub.execute_input":"2022-01-31T22:16:54.802148Z","iopub.status.idle":"2022-01-31T22:17:04.132452Z","shell.execute_reply.started":"2022-01-31T22:16:54.802121Z","shell.execute_reply":"2022-01-31T22:17:04.131764Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#plot a slice from the example exam\nimport matplotlib.pyplot as plt\n\nplt.imshow(np.squeeze(exam[:, :, 20]), cmap=\"gray\")","metadata":{"execution":{"iopub.status.busy":"2022-01-31T19:37:04.076758Z","iopub.execute_input":"2022-01-31T19:37:04.077015Z","iopub.status.idle":"2022-01-31T19:37:04.282749Z","shell.execute_reply.started":"2022-01-31T19:37:04.076986Z","shell.execute_reply":"2022-01-31T19:37:04.282059Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model creation\n\nA 3-dimensional Convolution Neural Network is proposed, consisting of several modules of 3D conv, maxpool and batch normalization layers. The model was adapted from this [paper](https://arxiv.org/pdf/2007.13224.pdf) by Zunair et al.","metadata":{}},{"cell_type":"code","source":"import tensorflow as tf\nfrom tensorflow import keras\nfrom tensorflow.keras import layers\nfrom tensorflow.keras import initializers\n\n#initializer = tf.keras.initializers.RandomNormal(mean=3., stddev=1.)\n\ninputs = keras.Input((128, 128, 64, 1))\n\nx = layers.Conv3D(filters=64, kernel_size=3, activation=\"relu\")(inputs)\nx = layers.MaxPool3D(pool_size=2)(x)\nx = layers.BatchNormalization()(x)\n\nx = layers.Conv3D(filters=64, kernel_size=3, activation=\"relu\")(x)\nx = layers.MaxPool3D(pool_size=2)(x)\nx = layers.BatchNormalization()(x)\n\nx = layers.Conv3D(filters=128, kernel_size=3, activation=\"relu\")(x)\nx = layers.MaxPool3D(pool_size=2)(x)\nx = layers.BatchNormalization()(x)\n\nx = layers.Conv3D(filters=256, kernel_size=3, activation=\"relu\")(x)\nx = layers.MaxPool3D(pool_size=2)(x)\nx = layers.BatchNormalization()(x)\n\nx = layers.GlobalAveragePooling3D()(x)\nx = layers.Dense(units=512, activation=\"relu\")(x)\nx = layers.Dropout(0.3)(x)\n\noutputs = layers.Dense(units=1, activation=\"sigmoid\")(x)\n\n# Define the model.\nmodel = keras.Model(inputs=inputs, outputs=outputs, name=\"3dcnn\")\n\n#initial_learning_rate = 0.0001\n#lr_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=0.001),\n    #metrics=[\"acc\"],\n    metrics=tf.keras.metrics.AUC()\n)\nmodel.summary()\nmodel.save('3d_CNN_rvlv.h5')\ndel model\nK.clear_session()\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-01-31T22:19:03.323534Z","iopub.execute_input":"2022-01-31T22:19:03.32408Z","iopub.status.idle":"2022-01-31T22:19:05.922492Z","shell.execute_reply.started":"2022-01-31T22:19:03.324044Z","shell.execute_reply":"2022-01-31T22:19:05.921661Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Training Model\n\ncustom_dcom_image_generator is a function that constructs the batch on which we will train and test. If train, it will yield (return a generator, iterable only once) with the X and Y, where X is a group of images If test, it will yield the X on which we will predict.","metadata":{}},{"cell_type":"code","source":"def custom_dcom_image_generator(batch_size, dataset, test=False, debug=False):\n    \n    fnames = dataset[['StudyInstanceUID', 'SeriesInstanceUID']]\n    \n    if not test:\n        Y = dataset[['rv_lv_ratio_lt_1']]\n        prefix = 'input/rsna-str-pulmonary-embolism-detection/train'\n        \n    else:\n        prefix = 'input/rsna-str-pulmonary-embolism-detection/train'\n    \n    X = np.array([])\n    batch = 0\n    for st, sr in fnames.values:\n        if debug:\n            print(f\"Current file: ../{prefix}/{st}/{sr}\")\n\n        temp = get_exam(f\"../{prefix}/{st}/{sr}\")\n        temp = temp.reshape((1, temp.shape[0], temp.shape[1], temp.shape[2], temp.shape[3]))\n        \n        if X.size==0:\n            X = temp\n        else:\n            X = np.concatenate((X, temp))\n            \n        \n        del st, sr\n        \n        #If we reached the end of the batch\n        if len(X) == batch_size:\n            if test:\n                #yield is used to save memory\n                yield X\n                del X\n            else:\n                yield X, Y[batch*batch_size:(batch+1)*batch_size].values\n                del X\n                \n            gc.collect()\n            X = np.array([])\n            batch += 1\n        \n    if test:\n        yield X\n    else:\n        yield X, Y[batch*batch_size:(batch+1)*batch_size].values\n        del Y\n    del X\n    gc.collect()\n    return","metadata":{"execution":{"iopub.status.busy":"2022-01-31T22:17:37.239001Z","iopub.execute_input":"2022-01-31T22:17:37.239514Z","iopub.status.idle":"2022-01-31T22:17:37.249607Z","shell.execute_reply.started":"2022-01-31T22:17:37.239471Z","shell.execute_reply":"2022-01-31T22:17:37.248763Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def train_model(model_path, train_data, train_size, batch_size, max_train_time, debug):\n    #Train loop\n    for n, (x, y) in enumerate(custom_dcom_image_generator(batch_size, train_data, False, debug)):\n\n        if len(x) < batch_size:#Tries to filter out empty or short data\n            model = load_model(model_path)\n            break\n        \n        clear_output(wait=True)\n        print(\"Training batch: %i - %i\" %(batch_size*n, batch_size*(n+1)))\n        model = load_model(model_path)\n        #print(\"check 1\")\n        hist = model.fit(\n            x[:train_size], \n            y[:train_size],\n\n            callbacks = checkpoint,\n\n            #validation_split=0.2,\n            epochs=50,\n            #batch_size=2,\n            verbose=debug\n        )\n        print(hist)\n        \n        try:\n            history = np.concatenate(history, hist.history)\n        except:\n            history = hist.history\n        #print(history)\n            \n\n        print(\"Metrics for batch validation:\")\n        model.evaluate(x[train_size:],\n                       y[train_size:]\n                      )\n\n        #To make sure that our model doesn't train overtime\n        if time.time() - start >= max_train_time:\n            print(\"Time's up!\")\n            break\n\n        model.save('3d_CNN_rvlv.h5')\n        del model, x, y, hist\n        #del x, y, hist\n        K.clear_session()\n        gc.collect()\n    \n    \n    return history, model\n    #return model","metadata":{"execution":{"iopub.status.busy":"2022-01-31T22:17:40.830765Z","iopub.execute_input":"2022-01-31T22:17:40.831026Z","iopub.status.idle":"2022-01-31T22:17:40.840215Z","shell.execute_reply.started":"2022-01-31T22:17:40.830996Z","shell.execute_reply":"2022-01-31T22:17:40.839225Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"history= np.array([])\nstart = time.time()\ntrain_data=train.sample(frac=1)\nbatch_size = 10\ndebug = 0\n\ntrain_size = int(batch_size*0.1)\n\nmax_train_time = 3600 * 4 #hours to seconds of training\n\ncheckpoint = MC(filepath='../working/3d_CNN_rvlv.h5', monitor='val_loss', save_best_only=True, verbose=1)\n#Train loop\n#history,trained_model = train_model('../working/3d_CNN_rvlv.h5', train_data, train_size, batch_size, max_train_time, debug)\nhistory, trained_model = train_model('../working/3d_CNN_rvlv.h5', train_data, train_size, batch_size, max_train_time, debug)\ntrained_model.save('3d_CNN_rvlv_trained.h5')","metadata":{"execution":{"iopub.status.busy":"2022-01-31T22:19:14.153807Z","iopub.execute_input":"2022-01-31T22:19:14.154064Z","iopub.status.idle":"2022-01-31T23:41:12.285729Z","shell.execute_reply.started":"2022-01-31T22:19:14.154035Z","shell.execute_reply":"2022-01-31T23:41:12.284885Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"history","metadata":{"execution":{"iopub.status.busy":"2022-01-31T23:41:12.287184Z","iopub.execute_input":"2022-01-31T23:41:12.287437Z","iopub.status.idle":"2022-01-31T23:41:12.298335Z","shell.execute_reply.started":"2022-01-31T23:41:12.287387Z","shell.execute_reply":"2022-01-31T23:41:12.297402Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"predictions = {}\nstopper = 3600 * 4 #4 hours limit for prediction\npred_start_time = time.time()\\\n\np, c = time.time(), time.time()\nbatch_size = 1\nl = 0\nn = test.shape[0]\n\nfor x in custom_dcom_image_generator(batch_size, test, True, False):\n    clear_output(wait=True)\n    #model = load_model(\"../input/modelss/3d_CNN_rvlv_trained.h5\")\n    model = load_model(\"../working/3d_CNN_rvlv.h5\")\n    print(x.shape)\n    preds = model.predict(x, verbose=1)\n    \n    try:\n        print(preds)\n        predictions += preds.tolist()\n    \n    except Exception as e:\n        print(e)\n        predictions = preds.tolist()\n            \n            \n    l = (l+batch_size)%n\n    p, c = c, time.time()\n    print(\"One batch time: %.2f seconds\" %(c-p))\n    print(\"ETA: %.2f\" %((n-l)*(c-p)/batch_size))\n    \n    if c - pred_start_time >= stopper:\n        print(\"Time's up!\")\n        break\n    \n    del model\n    K.clear_session()\n    \n    del x, preds\n    gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-01-31T20:43:03.573682Z","iopub.execute_input":"2022-01-31T20:43:03.574161Z","iopub.status.idle":"2022-01-31T20:43:18.394925Z","shell.execute_reply.started":"2022-01-31T20:43:03.574112Z","shell.execute_reply":"2022-01-31T20:43:18.393748Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"predictions","metadata":{"execution":{"iopub.status.busy":"2022-01-31T00:26:51.239414Z","iopub.execute_input":"2022-01-31T00:26:51.239705Z","iopub.status.idle":"2022-01-31T00:26:51.248353Z","shell.execute_reply.started":"2022-01-31T00:26:51.239673Z","shell.execute_reply":"2022-01-31T00:26:51.247572Z"},"trusted":true},"execution_count":null,"outputs":[]}]}