{"cells":[{"metadata":{"_uuid":"c7c6a863-09f9-4124-b7e8-ad3b44b00e0e","_cell_guid":"dbbfd3f3-1a25-4999-a264-eceb20a4931b","trusted":true},"cell_type":"code","source":"# -*- coding: utf-8 -*-\n\"\"\"final-submission.ipynb\n\nAutomatically generated by Colaboratory.\n\nOriginal file is located at\n    https://colab.research.google.com/drive/1cNndOCXSkssZtnELJgpebPyzV3J8dqy_\n\"\"\"\n\n# Installing packages required\n!conda install -c conda-forge gdcm -y\n!pip install git+https://github.com/titu1994/keras-efficientnets.git\n!pip install keras_applications==1.0.8 --no-deps\n!pip install keras==2.2.4\n# !pip install keras-preprocessing==1.2\n\n# Imports\nimport numpy as np \nimport pandas as pd\nimport os\nimport gc\nimport time\nfrom IPython.display import clear_output\nimport random\n\n\nfrom tensorflow.keras.models import load_model\nfrom tensorflow.keras.callbacks import ModelCheckpoint as MC\nfrom tensorflow.keras import backend as K\nimport tensorflow as tf\nimport keras\nfrom tensorflow import keras\nfrom keras import applications\nfrom tensorflow.keras.models import Model, Sequential\nfrom tensorflow.keras.layers import Input, Dense, Dropout, Conv2D\n\n# Image processing imports\nfrom scipy import ndimage\nimport matplotlib.pyplot as plt\nimport vtk\nfrom vtk.util import numpy_support\nimport cv2\n\n# Imports for working with DICOM Images\nimport pydicom\nimport scipy.ndimage\nimport gdcm\n\nfrom os import listdir, mkdir\n\n# Kaggle working directory\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)\n\n\"\"\"## Importing and Checking Datasets\n\nTraining Data\n\"\"\"\n\nprint('Reading train data...')\ntrain = pd.read_csv(\"../input/rsna-str-pulmonary-embolism-detection/train.csv\")\nprint(train.shape)\ntrain.head()\n\n\"\"\"Test Data\"\"\"\n\nprint('Reading test data...')\ntest = pd.read_csv(\"../input/rsna-str-pulmonary-embolism-detection/test.csv\")\nprint(test.shape)\ntest.head()\n\n\"\"\"Sample Submission File for Final Calculation\"\"\"\n\nprint('Reading sample data...')\nss = pd.read_csv(\"../input/rsna-str-pulmonary-embolism-detection/sample_submission.csv\")\nprint(ss.shape)\nss.head()\n\n\"\"\"Checking the IDs and Keys\"\"\"\n\nids = ss.id\ncounts = [1 for _ in range(10)]\ntypes = []\n\nfor i in ids:\n    n = '_'.join(i.split('_')[1:])\n    if n not in types:\n        types.append(n)\n    else:\n        counts[types.index(n)] += 1\n        \nfor x in range(len(counts)):\n    print(str(types[x]) + \" ------> \" + str(counts[x]) )\n\n# Helper Functions\n\n# Function to get DICOM Image Array","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"164b7ac7-92a5-4352-8651-b6d8229c85c3","_cell_guid":"fc41688a-657f-4abd-86c2-f67dad7069a5","trusted":true},"cell_type":"code","source":"\ndef load_scans(dcm_path):\n    '''\n    Reads the DICOM Image from a path and extracts a slice\n    '''   \n    slice_new = pydicom.dcmread(dcm_path)\n    return slice_new\n\ndef convert_to_rgb(array):\n    '''\n    Converts a dicom array to RGB values\n    '''\n    array = array.reshape((512, 512, 1))\n    \n    return np.stack([array, array, array], axis=2).reshape((512, 512, 3))\n\n# Image Pre-processing and Augmentation\n\ndef random_crop(img):\n    '''\n    Randmly crops image to get it to a desired dimension and returns the cropped image\n    '''\n    DIM = 320\n    \n    new_image = tf.image.random_crop(img, size=[DIM, DIM])\n    arr = new_image.numpy()\n    \n    return np.stack([arr, arr, arr], axis=2).reshape((DIM, DIM, 3))\n\n\ndef random_rotation(img, angle=8):\n    '''\n    Rotates the image randomly in positive or negative direction depending on specified angle range\n    '''\n    degree = random.uniform(0, 1)*angle\n    orientation = random.choice([-1, 1])\n    new_image = ndimage.rotate(img, degree*orientation, reshape=False)\n    \n    return new_image\n\n\ndef random_flip(img, threshold=0.5):\n    '''\n    Based on a uniform probability distribution and specified threshold, \n    horizontally flips the image \n    '''\n    \n    chance = random.uniform(0,1)\n    \n    if chance > threshold:\n        new_img = tf.image.flip_left_right(img)\n    else:\n        new_img = img\n    \n    return new_img\n\ndef change_contrast(img, factor=2):\n    '''\n    Change contrast of the image\n    '''\n    \n    new_img = tf.image.adjust_contrast(img, factor)\n    \n    return new_img\n\n# Windowing\n\n# Set the threshold of the pixels in image\ndef set_outside_scanner_to_air(raw_pixelarrays, neg_thresh=-1000):\n    '''\n    In OSIC we find outside-scanner-regions with raw-values of -2000. \n    This function thresholds between air (0) and this default (-2000) using -1000\n    ''' \n    raw_pixelarrays[raw_pixelarrays <= neg_thresh] = 0\n    return raw_pixelarrays\n\ndef transform_to_hu(slices):\n    '''\n    Convert the Image to HU and return the 16 bit image\n    '''\n    \n    image = slices.pixel_array.astype(np.int16)\n    \n    image = set_outside_scanner_to_air(image)\n    \n    b = slices.RescaleIntercept\n    m = slices.RescaleSlope\n    \n    if b != 1:\n            image = m * image.astype(np.float64)\n            image = image.astype(np.int16)\n    \n    return np.array(image, dtype=np.int16)\n\n\ndef set_manual_window(hu_image, custom_center, custom_width):\n    '''\n    Window the image based on manually defined center and width points\n    '''\n    \n    window_image = hu_image.copy()\n    \n    min_value = custom_center - (custom_width/2)\n    max_value = custom_center + (custom_width/2)\n    \n    window_image[window_image < min_value] = min_value\n    window_image[window_image > max_value] = max_value\n    \n    return window_image\n    \n\ndef three_channel_stacker(hu_scan, r_window=[40,400], g_window=[100, 700], b_window=[-600,1500], DIM=512 ):\n    '''\n    Create an RGB 3-channel Image from the HU filtered image  using windowing for each channel\n    '''\n    \n    r_channel = set_manual_window(hu_scan, r_window[0], r_window[1] )\n    g_channel = set_manual_window(hu_scan, g_window[0], g_window[1])\n    b_channel = set_manual_window(hu_scan, b_window[0], b_window[1])\n    \n    r_channel = r_channel.reshape((DIM, DIM, 1))\n    g_channel = g_channel.reshape((DIM, DIM, 1))\n    b_channel = b_channel.reshape((DIM, DIM, 1))\n    \n    return np.stack([r_channel, g_channel, b_channel], axis=2).reshape((DIM, DIM, 3))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"992fe12c-aca5-4e83-8aa7-48c9a73408bd","_cell_guid":"4ba2008a-20f2-4a50-abdb-73c75943cb2f","trusted":true},"cell_type":"code","source":"\"\"\"Testing out Pre-Processing\"\"\"\n\nexample = \"../input/rsna-str-pulmonary-embolism-detection/train/6897fa9de148/2bfbb7fd2e8b/822dd7790999.dcm\"\n\nds = load_scans(example)\nscan = load_scans(example)\nhu_scan = transform_to_hu(scan)\nfinal_image = three_channel_stacker(hu_scan)\n\nf, ax = plt.subplots(1,4, figsize=(30,30))\n# ax[0].imshow(ds)\nax[0].imshow(hu_scan)\nax[1].imshow(final_image)\nax[2].imshow(random_flip(final_image))\nax[3].imshow(change_contrast(random_flip(random_rotation(final_image)), 2))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1d0ef2b7-ac9c-4e45-8970-0decdaf78259","_cell_guid":"5a04d04d-7ace-486c-939c-d9c62af8f426","trusted":true},"cell_type":"code","source":"\n\"\"\"## Model Creation\"\"\"\n\ninputs = Input((512, 512, 3))\n\nbase_model = keras.applications.Xception(include_top=False, weights=\"imagenet\")\nbase_model.trainable = False\n\noutputs = base_model(inputs, training=False)\noutputs = keras.layers.GlobalAveragePooling2D()(outputs)\noutputs = Dropout(0.25)(outputs)\noutputs = Dense(1024, activation='relu')(outputs)\noutputs = Dense(256, activation='relu')(outputs)\noutputs = Dense(64, activation='relu')(outputs)\n\npe_present_on_image = Dense(1, activation='sigmoid', name='pe_present_on_image')(outputs)\nrv_lv_ratio_gte_1 = Dense(1, activation='sigmoid', name='rv_lv_ratio_gte_1')(outputs)\nrv_lv_ratio_lt_1 = Dense(1, activation='sigmoid', name='rv_lv_ratio_lt_1')(outputs) \nleftsided_pe = Dense(1, activation='sigmoid', name='leftsided_pe')(outputs)\nchronic_pe = Dense(1, activation='sigmoid', name='chronic_pe')(outputs)\nrightsided_pe = Dense(1, activation='sigmoid', name='rightsided_pe')(outputs)\nacute_and_chronic_pe = Dense(1, activation='sigmoid', name='acute_and_chronic_pe')(outputs)\ncentral_pe = Dense(1, activation='sigmoid', name='central_pe')(outputs)\nindeterminate = Dense(1, activation='sigmoid', name='indeterminate')(outputs)\n\nmodel = Model(inputs=inputs, outputs={'pe_present_on_image':pe_present_on_image,\n                                      'rv_lv_ratio_gte_1':rv_lv_ratio_gte_1,\n                                      'rv_lv_ratio_lt_1':rv_lv_ratio_lt_1,\n                                      'leftsided_pe':leftsided_pe,\n                                      'chronic_pe':chronic_pe,\n                                      'rightsided_pe':rightsided_pe,\n                                      'acute_and_chronic_pe':acute_and_chronic_pe,\n                                      'central_pe':central_pe,\n                                      'indeterminate':indeterminate\n                                     })\n\nopt = keras.optimizers.Adam(lr=0.001)\n\nmodel.compile(optimizer=opt, loss='binary_crossentropy', metrics=['accuracy'])\n\nmodel.summary()\nmodel.save('improved_model.h5')\ndel model\nK.clear_session()\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e696bb71-7172-4f5b-a7d9-d6c84c4a99ac","_cell_guid":"07c08e92-0dfb-4c7d-9548-46bc13cc845b","trusted":true},"cell_type":"code","source":"\n\n\n\"\"\"## Image Generator\"\"\"\n\ndef data_preprocessor(batch_size, dataset, test=False, debug=True):\n    \n    scan_ids = dataset[['StudyInstanceUID', 'SeriesInstanceUID', 'SOPInstanceUID']]\n\n    if not test:\n        Y = dataset[['pe_present_on_image', \n                     'rv_lv_ratio_gte_1', \n                     'rv_lv_ratio_lt_1', \n                     'leftsided_pe',\n                     'chronic_pe', \n                     'rightsided_pe', \n                     'acute_and_chronic_pe', \n                     'central_pe', \n                     'indeterminate'\n                    ]]\n        location = 'input/rsna-str-pulmonary-embolism-detection/train'\n        \n    else:\n        location = 'input/rsna-str-pulmonary-embolism-detection/test'\n    \n    X = []\n    batch = 0\n    for study, series, sliceOfPatient in scan_ids.values:\n\n        # Extract the location of the images\n        if debug:\n            print(f\"Current file: ../{location}/{study}/{series}/{sliceOfPatient}.dcm\")\n        \n        \n        # Extract the HU frames and window the image\n        scan = load_scans(f\"../{location}/{study}/{series}/{sliceOfPatient}.dcm\")\n        image = transform_to_hu(scan)\n        image = three_channel_stacker(image)\n        \n        # Apply the augmentations\n        image = random_flip(image)\n        image = random_rotation(image)\n        image = change_contrast(image)\n\n        \n        X.append(image)\n        del study, series, sliceOfPatient\n        \n        if len(X) == batch_size:\n            if test:\n                yield np.array(X)\n                del X\n            else:\n                yield np.array(X), Y[batch*batch_size:(batch+1)*batch_size].values\n                del X\n                \n            gc.collect()\n            X = []\n            batch += 1\n        \n    if test:\n        yield np.array(X)\n    else:\n        yield np.array(X), Y[batch*batch_size:(batch+1)*batch_size].values\n        del Y\n    del X\n    gc.collect()\n    return","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"36da1fa7-0e65-411f-82ac-adc04458a9e8","_cell_guid":"2a59d068-98a8-4c0a-9212-1043c508dd1b","trusted":true},"cell_type":"code","source":"\n\"\"\"## Training\"\"\"\n\nhistory = {}\nstart = time.time()\ndebug = 0\nbatch_size = 1000\ntrain_size = int(batch_size*0.9)    # Use 90% of the batch for training\n\n# Set the maxout of training time \nmax_training_hours = 4.5  \nmax_train_time = 3600 * max_training_hours            \n\ncheckpoint = MC(filepath='../working/improved_model.h5', monitor='val_loss', save_best_only=True, verbose=1)\n\n#Train loop\nfor n, (x, y) in enumerate(data_preprocessor(batch_size, train.sample(frac=1), test=False, debug=False)):\n    \n    if len(x) < 10: #Tries to filter out empty or short data\n        break\n        \n    clear_output(wait=True)\n    \n    print(\"Training batch: %i - %i\" %(batch_size*n, batch_size*(n+1)))\n    \n    model = load_model('../working/improved_model.h5')\n    hist = model.fit(\n        x[:train_size], #Y values are in a dict as there's more than one target for training output\n        {\n         'pe_present_on_image':y[:train_size, 0],\n         'rv_lv_ratio_gte_1':y[:train_size, 1],\n         'rv_lv_ratio_lt_1':y[:train_size, 2],\n         'leftsided_pe':y[:train_size, 3],\n         'chronic_pe':y[:train_size, 4],\n         'rightsided_pe':y[:train_size, 5],\n         'acute_and_chronic_pe':y[:train_size, 6],\n         'central_pe':y[:train_size, 7],\n         'indeterminate':y[:train_size, 8]\n        },\n\n        callbacks = checkpoint,\n\n        validation_split=0.2,\n        epochs=3,\n        batch_size=8,\n        verbose=debug\n    )\n    \n    print(\"Metrics for batch validation:\")\n    model.evaluate(x[train_size:], {\n                                    'pe_present_on_image':y[train_size:, 0],\n                                    'rv_lv_ratio_gte_1':y[train_size:, 1],\n                                    'rv_lv_ratio_lt_1':y[train_size:, 2],\n                                    'leftsided_pe':y[train_size:, 3],\n                                    'chronic_pe':y[train_size:, 4],\n                                    'rightsided_pe':y[train_size:, 5],\n                                    'acute_and_chronic_pe':y[train_size:, 6],\n                                    'central_pe':y[train_size:, 7],\n                                    'indeterminate':y[train_size:, 8]\n                                   }\n                   \n                  )\n    \n    try:\n        for key in hist.history.keys():\n            history[key] = np.concatenate([history[key], hist.history[key]], axis=0)\n    except:\n        for key in hist.history.keys():\n            history[key] = hist.history[key]\n            \n    #To make sure that our model don't train overtime\n    if time.time() - start >= max_train_time:\n        print(\"Max Training time Reached!\")\n        break\n        \n    model.save('improved_model.h5')\n    \n    del model, x, y, hist\n    \n    K.clear_session()\n    gc.collect()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e4e6d182-a5c8-415c-9a85-fa1594c45e5e","_cell_guid":"36bc0755-b9b7-4cd4-9977-6eb5c892ee33","trusted":true},"cell_type":"code","source":"\n\"\"\"Plot the Losses and Accuracies\"\"\"\n\nfor key in history.keys():\n    if key.startswith('val'):\n        continue\n    else:\n        epoch = range(len(history[key]))\n        plt.plot(epoch, history[key]) #X=epoch, Y=value\n        plt.plot(epoch, history['val_'+key])\n        plt.title(key)\n        if 'accuracy' in key:\n            plt.axis([0, len(history[key]), -0.1, 1.1]) #Xmin, Xmax, Ymin, Ymax\n        plt.legend(['train', 'validation'], loc='upper right')\n        plt.show()","execution_count":null,"outputs":[]}],"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":4,"nbformat_minor":4}