{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Peek\n\nFirst, let's see the structure of our files."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"!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\nimport 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\n\nimport vtk\nfrom vtk.util import numpy_support\nimport cv2\n\n\nimport pydicom\nimport scipy.ndimage\nimport gdcm\n\nfrom os import listdir, mkdir\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)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Load Data\n\nNext, we load all '.csv' files into memory and peek into their makeup."},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"print('Reading train data...')\ntrain = pd.read_csv(\"../input/rsna-str-pulmonary-embolism-detection/train.csv\")\nprint(train.shape)\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Reading test data...')\ntest = pd.read_csv(\"../input/rsna-str-pulmonary-embolism-detection/test.csv\")\nprint(test.shape)\ntest.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Reading sample data...')\nss = pd.read_csv(\"../input/rsna-str-pulmonary-embolism-detection/sample_submission.csv\")\nprint(ss.shape)\nss.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train['weight'] = train.groupby('StudyInstanceUID')\\\n                       .pe_present_on_image.transform('mean')\nX_train = train.loc[train.weight>0]\nprint( X_train.shape[0] / train.shape[0] )","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Check Targets and Input Image\n\nTo make sure that we will be making the correct training targets, we check all ID from the sample submission."},{"metadata":{"trusted":true},"cell_type":"code","source":"ids = ss.id\ncounter = [1 for _ in range(10)]\nmapper = []\nfor i in ids:\n    n = '_'.join(i.split('_')[1:])\n    if n not in mapper:\n        mapper.append(n)\n    else:\n        counter[mapper.index(n)] += 1\nprint(\"List of keys:\")\nprint(mapper, sep='\\n')\nprint()\nprint(\"Count of items per key:\")\nprint(counter)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"After checking the keys, we will now work on the fuction to get the image array from a DICOM image. We will be using a code snippet by [eladwar](https://www.kaggle.com/eladwar) from his notebook [here](https://www.kaggle.com/eladwar/20-seconds-or-less)."},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_scans(dcm_path):\n    # otherwise we sort by ImagePositionPatient (z-coordinate) or by SliceLocation\n#     slices = [pydicom.dcmread(dcm_path + \"/\" + file) for file in listdir(dcm_path)]\n#     slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))\n#     return slices\n    slice_new = pydicom.dcmread(dcm_path)\n    return slice_new\n# for file in listdir(\"../input/rsna-str-pulmonary-embolism-detection/train/6897fa9de148/2bfbb7fd2e8b\"):\n#     print(file)\nexample = \"../input/rsna-str-pulmonary-embolism-detection/train/6897fa9de148/2bfbb7fd2e8b/822dd7790999.dcm\"\nscans = load_scans(example)\nprint(scans.pixel_array)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def random_crop(img):\n    new_image = tf.image.random_crop(img, size=[320, 320])\n    return new_image.numpy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy import ndimage\nimport random\ndef random_rotation(img):\n#     new_image = tf.keras.preprocessing.image.random_rotation(img, 0.16)\n#     return new_image\n    degree = random.uniform(0, 1)*8\n    orientation = random.choice([-1, 1])\n    new_image = ndimage.rotate(img, degree*orientation, reshape=False)\n    return new_image\n\n    \n        \n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import tensorflow as tf\n\ndef convert_to_rgb(array):\n    array = array.reshape((512, 512, 1))\n    return np.stack([array, array, array], axis=0).reshape((512, 512, 3))\n\n#test read a dcom file and view it\nfpath = \"../input/rsna-str-pulmonary-embolism-detection/train/0003b3d648eb/d2b2960c2bbf/00ac73cfc372.dcm\"\n\nimport matplotlib.pyplot as plt\n# img = scans[10].pixel_array\nimg = scans.pixel_array\nnew_image = random_crop(img)\nnew_image = convert_to_rgb(img)\n# new_image = random_rotation(new_image)\n\n\nf, ax = plt.subplots(1,2, figsize=(16,20))\nax[0].imshow(img)\nax[1].imshow(new_image)\nprint(img.shape)\nprint(new_image.shape)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"After defining our image reader, we will test it with a sample DICOM image to load."},{"metadata":{},"cell_type":"markdown","source":"# Model Creation"},{"metadata":{},"cell_type":"markdown","source":"I decided to try and use a pre-trained Xception model as a feature-extractor. To simplify my coding, I did a multi-output model with each output having one node activated by a sigmoid. Also, since we are dealing with numbers between 1 and 0, I decided to use binary_crossentropy as the loss function."},{"metadata":{"trusted":true},"cell_type":"code","source":"# !pip install keras_applications==1.0.8 --no-deps\n# !pip install keras==2.2.1\n!pip install efficientnet\n\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\nimport efficientnet.tfkeras as efn \n\n# inputs = Input((512, 512, 3))\n# #x = Conv2D(3, (1, 1), activation='relu')(inputs)\n# base_model = keras.applications.Xception(\n#     include_top=False,\n#     weights=\"imagenet\"\n# )\ninp = keras.Input(shape=(320, 320, 1)) # INPUT IS UINT8\nx = keras.layers.Concatenate()([inp/255., inp/255., inp/255.])\nbase_model = efn.EfficientNetB4(weights='imagenet', include_top=False)\n\nbase_model.trainable = False\n\noutputs = base_model(x, training=False)\noutputs = keras.layers.GlobalAveragePooling2D()(outputs)\noutputs = Dropout(0.5)(outputs)\noutputs = Dense(1024, activation='relu')(outputs)\noutputs = Dense(256, activation='relu')(outputs)\noutputs = Dense(64, activation='relu')(outputs)\nppoi = Dense(1, activation='sigmoid', name='pe_present_on_image')(outputs)\nrlrg1 = Dense(1, activation='sigmoid', name='rv_lv_ratio_gte_1')(outputs)\nrlrl1 = Dense(1, activation='sigmoid', name='rv_lv_ratio_lt_1')(outputs) \nlspe = Dense(1, activation='sigmoid', name='leftsided_pe')(outputs)\ncpe = Dense(1, activation='sigmoid', name='chronic_pe')(outputs)\nrspe = Dense(1, activation='sigmoid', name='rightsided_pe')(outputs)\naacpe = Dense(1, activation='sigmoid', name='acute_and_chronic_pe')(outputs)\ncnpe = Dense(1, activation='sigmoid', name='central_pe')(outputs)\nindt = Dense(1, activation='sigmoid', name='indeterminate')(outputs)\n\nmodel = Model(inputs=inp, outputs={'pe_present_on_image':ppoi,\n                                      'rv_lv_ratio_gte_1':rlrg1,\n                                      'rv_lv_ratio_lt_1':rlrl1,\n                                      'leftsided_pe':lspe,\n                                      'chronic_pe':cpe,\n                                      'rightsided_pe':rspe,\n                                      'acute_and_chronic_pe':aacpe,\n                                      'central_pe':cnpe,\n                                      'indeterminate':indt})\n\nopt = keras.optimizers.Adam(lr=0.0001)\n\nmodel.compile(optimizer=opt,\n              loss='binary_crossentropy',\n              metrics=['accuracy'])\n\nmodel.summary()\nmodel.save('pe_detection_model.h5')\ndel model\nK.clear_session()\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Training Model"},{"metadata":{},"cell_type":"markdown","source":"Before we can train our model, we would be needing an image generator. This is so that our training code would look much cleaner and nicer."},{"metadata":{"trusted":true},"cell_type":"code","source":"from vtk.util import numpy_support\nimport cv2\n\nreader = vtk.vtkDICOMImageReader()\ndef get_img(path):\n    reader.SetFileName(path)\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    return ArrayDicom","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def custom_dcom_image_generator(batch_size, dataset, test=False, debug=True):\n    \n    fnames = dataset[['StudyInstanceUID', 'SeriesInstanceUID', 'SOPInstanceUID']]\n    \n    if not test:\n        Y = dataset[['pe_present_on_image', 'rv_lv_ratio_gte_1', 'rv_lv_ratio_lt_1', 'leftsided_pe',\n                     'chronic_pe', 'rightsided_pe', 'acute_and_chronic_pe', 'central_pe', 'indeterminate'\n                    ]]\n        prefix = 'input/rsna-str-pulmonary-embolism-detection/train'\n        \n    else:\n        prefix = 'input/rsna-str-pulmonary-embolism-detection/test'\n    \n    X = []\n    batch = 0\n    for st, sr, so in fnames.values:\n        if debug:\n            print(f\"Current file: ../{prefix}/{st}/{sr}/{so}.dcm\")\n\n        dicom = load_scans(f\"../{prefix}/{st}/{sr}/{so}.dcm\")\n        image = random_crop(dicom.pixel_array)\n        image = random_rotation(image)\n        X.append(image)\n        \n        del st, sr, so\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":{"trusted":true},"cell_type":"code","source":"history = {}\nstart = time.time()\ndebug = 0\nbatch_size = 1000\ntrain_size = int(batch_size*0.9)\nsample_wght = train.groupby('StudyInstanceUID').pe_present_on_image.transform('mean') * 5.6222\nmax_train_time = 3600 * 3 #hours to seconds of training\n\ncheckpoint = MC(filepath='../working/pe_detection_model.h5', monitor='val_loss', save_best_only=True, verbose=1)\n#Train loop\nfor n, (x, y) in enumerate(custom_dcom_image_generator(batch_size, X_train.sample(frac=1), False, debug)):\n    \n    if len(x) < 10: #Tries to filter out empty or short data\n        break\n        \n    clear_output(wait=True)\n    print(\"Training batch: %i - %i\" %(batch_size*n, batch_size*(n+1)))\n    \n    model = load_model('../working/pe_detection_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        {'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        callbacks = checkpoint,\n        sample_weight = train.groupby('StudyInstanceUID').pe_present_on_image.transform('mean') * 5.6222,\n        validation_split=0.2,\n        epochs=5,\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    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(\"Time's up!\")\n        break\n        \n    model.save('pe_detection_model.h5')\n    del model, x, y, hist\n    K.clear_session()\n    gc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We will now look at the history of the training for our data. Since the data is trained across several different batches, there should be some spikes among the values reflected here."},{"metadata":{"trusted":true},"cell_type":"code","source":"for key in history.keys():\n    if key.startswith('val'):\n        continue\n    else:\n        epoch = range(len(history[key]))\n        plt.figure(num=None, figsize=(8, 6), dpi=80, facecolor='w', edgecolor='k')\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":{},"cell_type":"markdown","source":"# End of Part 1"},{"metadata":{},"cell_type":"markdown","source":"# Start of Part 2"},{"metadata":{},"cell_type":"markdown","source":"# Prediction\n\nWe will now proceed to predict our data. Since we can't submit this straight to the competition, due to using internet for downloading the Xception's ImageNet weights, this will be like a trial version for the prediction function and the csv output production."},{"metadata":{"trusted":true},"cell_type":"code","source":"# predictions = {}\n# stopper = 3600 * 4 #4 hours limit for prediction\n# pred_start_time = time.time()\n\n# p, c = time.time(), time.time()\n# batch_size = 500\n    \n# l = 0\n# n = test.shape[0]\n\n# for x in custom_dcom_image_generator(batch_size, test, True, False):\n#     clear_output(wait=True)\n#     model = load_model(\"../working/pe_detection_model.h5\")\n#     preds = model.predict(x, batch_size=8, verbose=1)\n    \n#     try:\n#         for key in preds.keys():\n#             predictions[key] += preds[key].flatten().tolist()\n            \n#     except Exception as e:\n#         print(e)\n#         for key in preds.keys():\n#             predictions[key] = preds[key].flatten().tolist()\n            \n#     l = (l+batch_size)%n\n#     print('Total predicted:', len(predictions['indeterminate']),'/', 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()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We will now check the shape of each key in predictions to make sure we predicted everything."},{"metadata":{"trusted":true},"cell_type":"code","source":"# for key in predictions.keys():\n#     print(key, np.array(predictions[key]).shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Next, we will convert predictions into a dataframe based on the test dataframe. We will also copy all unique **StudyInstanceUID** for predicting later."},{"metadata":{"trusted":true},"cell_type":"code","source":"# test_ids = []\n# for v in test.StudyInstanceUID:\n#     if v not in test_ids:\n#         test_ids.append(v)\n        \n# test_preds = test.copy()\n# test_preds = pd.concat([test_preds, pd.DataFrame(predictions)], axis=1)\n# test_preds.to_csv('test_predictions.csv', index=False)\n# test_preds","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here, since we need to predict for each label in each Study, we will aggregate the Mean for each group of the same label. We will also apply the label hierarchy for our final submission. *(Note: This is just a sample run for the submission)*"},{"metadata":{"trusted":true},"cell_type":"code","source":"# from scipy.special import softmax\n\n# label_agg = {key:[] for key in \n#              ['id', 'negative_exam_for_pe', 'rv_lv_ratio_gte_1',\n#               'rv_lv_ratio_lt_1', 'leftsided_pe', 'chronic_pe',\n#               'rightsided_pe', 'acute_and_chronic_pe',\n#               'central_pe', 'indeterminate']\n#             }\n\n# for uid in test_ids:\n#     temp = test_preds.loc[test_preds.StudyInstanceUID ==uid]\n#     label_agg['id'].append(uid)\n    \n#     n = temp.shape[0]\n#     #Check for any image level presence of PE of high confidence\n#     positive = any(temp.pe_present_on_image >= 0.5) #50% threshhold\n    \n#     #Only one from positive, negative and indeterminate should have value>0.5\n#     #per exam\n#     if positive: \n#         label_agg['indeterminate'].append(temp.indeterminate.min()/2)\n#         label_agg['negative_exam_for_pe'].append(0)\n#     else:\n#         if any(temp.indeterminate >= 0.5):\n#             label_agg['indeterminate'].append(temp.indeterminate.max())\n#             label_agg['negative_exam_for_pe'].append(1)\n#         else:\n#             label_agg['indeterminate'].append(temp.indeterminate.min()/2)\n#             label_agg['negative_exam_for_pe'].append(1)\n    \n#     #I decided that the total ratio should be equal to 1, so I used softmax\n#     #We modify the weights by multiplying the bigger by 2 and dividing the smaller by 2\n#     a, b = temp[['rv_lv_ratio_gte_1', 'rv_lv_ratio_lt_1']].mean().values\n#     if a > b:\n#         a, b = a*2, b/2\n#     elif a < b:\n#         a, b = a/2, b*2\n#     a, b = softmax([a, b])\n#     if positive:\n#         label_agg['rv_lv_ratio_gte_1'].append(a)\n#         label_agg['rv_lv_ratio_lt_1'].append(b)\n#     else:\n#         label_agg['rv_lv_ratio_gte_1'].append(a/2)\n#         label_agg['rv_lv_ratio_lt_1'].append(b/2)\n    \n#     #Next is for Chronic (C), Acute-Chronic (AC) and Acute (A) PE\n#     #We need to see if we got a high confidence value from either C or AC\n#     #If there is, we add it to a 50% based score for high confidence\n#     #and half weight for low confidence score\n#     if any(temp['acute_and_chronic_pe'] > 0.5): #50% confidence level\n#         label_agg['acute_and_chronic_pe'].append(0.5 + temp['acute_and_chronic_pe'].mean()/2)\n#         label_agg['chronic_pe'].append(temp['chronic_pe'].mean()/2)\n        \n#     elif any(temp['chronic_pe'] > 0.5):\n#         label_agg['acute_and_chronic_pe'].append(temp['acute_and_chronic_pe'].mean()/2)\n#         label_agg['chronic_pe'].append(0.5 + temp['chronic_pe'].mean()/2)\n        \n#     else: #Else, we set both to half values, as we declare the A as the value\n#         label_agg['acute_and_chronic_pe'].append(temp['acute_and_chronic_pe'].mean()/2)\n#         label_agg['chronic_pe'].append(temp['chronic_pe'].mean()/2)\n    \n#     #for right, left, central, we use the same metric above\n#     for key in ['leftsided_pe', 'rightsided_pe', 'central_pe']:\n#         if positive:\n#             label_agg[key].append(0.5 + temp[key].mean()/2)\n#         else:\n#             label_agg[key].append(temp[key].mean()/2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Lastly, we need to add the result of \"pe_present_on_image\" that is 1-to-1 for each label. We will ensemble a new DataFrame from the exam ids and labels that we have."},{"metadata":{"trusted":true},"cell_type":"code","source":"# uid = []\n# labels = []\n# df = pd.DataFrame(label_agg)\n# for key in ['negative_exam_for_pe', 'rv_lv_ratio_gte_1', 'rv_lv_ratio_lt_1', 'leftsided_pe', 'chronic_pe',\n#             'rightsided_pe', 'acute_and_chronic_pe', 'central_pe', 'indeterminate']:\n#     for i in df.id:\n#         uid.append('_'.join([i, key]))\n#         labels.append(df.loc[df.id==i][key].values[0])\n# del df\n# gc.collect()\n\n# uid += test_preds.SOPInstanceUID.tolist()\n# labels += test_preds['pe_present_on_image'].tolist()\n\n# sub = pd.DataFrame({\"id\":uid, 'label':labels})\n# sub","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Finally, we will save our submission into a file. But first, to make sure that we fill up any unpredicted variables (so that there will be no NaN in our values), we will fill up all missing values with 0.2 as a placeholder."},{"metadata":{"trusted":true},"cell_type":"code","source":"# sub.fillna(0.2, inplace=True)\nsub.to_csv('submission.csv', index=False)","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}