{"cells":[{"metadata":{},"cell_type":"markdown","source":"#### 1.Importing Libraries and Loading Data"},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport pydicom\nimport pandas as pd\nfrom glob import glob\nimport os\nimport csv\nimport random\nimport pydicom\nfrom skimage import io\nfrom skimage import measure\nfrom skimage.transform import resize\n\nimport matplotlib.pyplot as plt\nfrom matplotlib.patches import Rectangle\nimport matplotlib.patches as patches\n%matplotlib inline\nimport tensorflow as tf\nfrom tensorflow import keras\n\ndet_class_path = '../input/rsna-pneumonia-detection-challenge/stage_2_detailed_class_info.csv'\nbbox_path = '../input/rsna-pneumonia-detection-challenge/stage_2_train_labels.csv'\ndicom_dir = '../input/rsna-pneumonia-detection-challenge/stage_2_train_images/'\n\n# load and shuffle filenames\ndicom_dir = '../input/rsna-pneumonia-detection-challenge/stage_2_train_images/'\nfilenames = os.listdir(dicom_dir)\nrandom.shuffle(filenames)\n# split into train and validation filenames\nn_valid_samples = 2560\ntrain_filenames = filenames[n_valid_samples:]\nvalid_filenames = filenames[:n_valid_samples]\nprint('n train samples', len(train_filenames))\nprint('n valid samples', len(valid_filenames))\nn_train_samples = len(filenames) - n_valid_samples","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**The most interesting group here is the No Lung Opacity / Not Normal since they are cases that look like opacity but are not. \nSo the first step might be to divide the test images into clear groups and then only perform the bounding box prediction on the suspicious images.**"},{"metadata":{"trusted":true},"cell_type":"code","source":"# empty dictionary\npneumonia_locations = {}\n# load table\nwith open(os.path.join('../input/rsna-pneumonia-detection-challenge/stage_2_train_labels.csv'), mode='r') as infile:\n    # open reader\n    reader = csv.reader(infile)\n    # skip header\n    next(reader, None)\n    # loop through rows\n    for rows in reader:\n        # retrieve information\n        filename = rows[0]\n        location = rows[1:5]\n        pneumonia = rows[5]\n        # if row contains pneumonia add label to dictionary\n        # which contains a list of pneumonia locations per filename\n        if pneumonia == '1':\n            # convert string to float to int\n            location = [int(float(i)) for i in location]\n            # save pneumonia location in dictionary\n            if filename in pneumonia_locations:\n                pneumonia_locations[filename].append(location)\n            else:\n                pneumonia_locations[filename] = [location]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"det_class_df = pd.read_csv(det_class_path)\nprint(det_class_df.shape[0], 'class infos loaded')\nprint(det_class_df['patientId'].value_counts().shape[0], 'patient cases')\ndet_class_df.groupby('class').size().plot.bar()\ndet_class_df.sample(3)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Load the Bounding Box Data\n**Here we show the bounding boxes**"},{"metadata":{"trusted":true},"cell_type":"code","source":"bbox_df = pd.read_csv(bbox_path)\nprint(bbox_df.shape[0], 'boxes loaded')\nprint(bbox_df['patientId'].value_counts().shape[0], 'patient cases')\nbbox_df.sample(3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# we first try a join and see that it doesn't work (we end up with too many boxes)\ncomb_bbox_df = pd.merge(bbox_df, det_class_df, how='inner', on='patientId')\nprint(comb_bbox_df.shape[0], 'combined cases')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Concatenate\n\n**We have to concatenate the two datasets and then we get class and target information on each region**"},{"metadata":{"trusted":true},"cell_type":"code","source":"comb_bbox_df = pd.concat([bbox_df, \n                        det_class_df.drop('patientId',1)], 1)\nprint(comb_bbox_df.shape[0], 'combined cases')\ncomb_bbox_df.sample(3)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Distribution of Boxes and Labels\nThe values below show the number of boxes and the patients that have that number."},{"metadata":{"trusted":true},"cell_type":"code","source":"box_df = comb_bbox_df.groupby('patientId').\\\n    size().\\\n    reset_index(name='boxes')\ncomb_box_df = pd.merge(comb_bbox_df, box_df, on='patientId')\nbox_df.\\\n    groupby('boxes').\\\n    size().\\\n    reset_index(name='patients')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"How are class and target related?\nI assume that all the Target=1 values fall in the Lung Opacity class, but it doesn't hurt to check."},{"metadata":{"trusted":true},"cell_type":"code","source":"comb_bbox_df.groupby(['class', 'Target']).size().reset_index(name='Patient Count')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Images\nNow that we have the boxes and labels loaded we can examine a few images."},{"metadata":{"trusted":true},"cell_type":"code","source":"image_df = pd.DataFrame({'path': glob(os.path.join(dicom_dir, '*.dcm'))})\nimage_df['patientId'] = image_df['path'].map(lambda x: os.path.splitext(os.path.basename(x))[0])\nprint(image_df.shape[0], 'images found')\nimg_pat_ids = set(image_df['patientId'].values.tolist())\nbox_pat_ids = set(comb_box_df['patientId'].values.tolist())\n# check to make sure there is no funny business\nassert img_pat_ids.union(box_pat_ids)==img_pat_ids, \"Patient IDs should be the same\"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image_bbox_df = pd.merge(comb_box_df, \n                         image_df, \n                         on='patientId',\n                        how='left').sort_values('patientId')\nprint(image_bbox_df.shape[0], 'image bounding boxes')\nimage_bbox_df.head(5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Enrich the image fields\nWe have quite a bit of additional data in the DICOM header we can easily extract to help learn more about the patient like their age, view position and gender which can make the model much more precise"},{"metadata":{"trusted":true},"cell_type":"code","source":"DCM_TAG_LIST = ['PatientAge', 'BodyPartExamined', 'ViewPosition', 'PatientSex']\ndef get_tags(in_path):\n    c_dicom = pydicom.read_file(in_path, stop_before_pixels=False)\n    tag_dict = {c_tag: getattr(c_dicom, c_tag, '') \n         for c_tag in DCM_TAG_LIST}\n    tag_dict['path'] = in_path\n    return pd.Series(tag_dict)\nimage_meta_df = image_df.apply(lambda x: get_tags(x['path']), 1)\n# show the summary\nimage_meta_df['PatientAge'] = image_meta_df['PatientAge'].map(int)\nimage_meta_df['PatientAge'].hist()\nimage_meta_df.drop('path',1).describe(exclude=np.number)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image_full_df = pd.merge(image_df,\n                         image_meta_df,\n                         on='path')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Create Sample Data Set\nWe create a sample dataset covering different cases, and number of boxes"},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_df = image_bbox_df.\\\n    groupby(['Target','class', 'boxes']).\\\n    apply(lambda x: x[x['patientId']==x.sample(1)['patientId'].values[0]]).\\\n    reset_index(drop=True)\nsample_df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Show the position and bounding box\nHere we can see the position (point) and the bounding box for each of the different image types"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, m_axs = plt.subplots(2, 3, figsize = (20, 10))\nfor c_ax, (c_path, c_rows) in zip(m_axs.flatten(),\n                    sample_df.groupby(['path'])):\n    c_dicom = pydicom.read_file(c_path)\n    c_ax.imshow(c_dicom.pixel_array, cmap='bone')\n    c_ax.set_title('{class}'.format(**c_rows.iloc[0,:]))\n    for i, (_, c_row) in enumerate(c_rows.dropna().iterrows()):\n        c_ax.plot(c_row['x'], c_row['y'], 's', label='{class}'.format(**c_row))\n        c_ax.add_patch(Rectangle(xy=(c_row['x'], c_row['y']),\n                                width=c_row['width'],\n                                height=c_row['height'], \n                                 alpha = 0.5))\n        if i==0: c_ax.legend()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Bounding Box Distribution\nHere we just look at the bounding box distribution to get a better idea how this looks over the whole dataset"},{"metadata":{"trusted":true},"cell_type":"code","source":"pos_bbox = image_bbox_df.query('Target==1')\npos_bbox.plot.scatter(x='x', y='y')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax1 = plt.subplots(1, 1, figsize = (10, 10))\nax1.set_xlim(0, 1024)\nax1.set_ylim(0, 1024)\nfor _, c_row in pos_bbox.sample(1000).iterrows():\n    ax1.add_patch(Rectangle(xy=(c_row['x'], c_row['y']),\n                 width=c_row['width'],\n                 height=c_row['height'],\n                           alpha=5e-3))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Show the boxes as segmentation\n\n**By showing them as segmentations we can get a better probability map for where the opacity regions are most likely to occur**"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Show the boxes themselves\nX_STEPS, Y_STEPS = 1024, 1024\nxx, yy = np.meshgrid(np.linspace(0, 1024, X_STEPS),\n           np.linspace(0, 1024, Y_STEPS), \n           indexing='xy')\nprob_image = np.zeros_like(xx)\nfor _, c_row in pos_bbox.sample(5000).iterrows():\n    c_mask = (xx>=c_row['x']) & (xx<=(c_row['x']+c_row['width']))\n    c_mask &= (yy>=c_row['y']) & (yy<=c_row['y']+c_row['height'])\n    prob_image += c_mask\nfig, ax1 = plt.subplots(1, 1, figsize = (10, 10))\nax1.imshow(prob_image, cmap='hot')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Overlay the Probability on a few images\nDoes the probability we calculate seem to make sense? or have we flipped something somewhere?"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, m_axs = plt.subplots(2, 3, figsize = (20, 10))\nfor c_ax, (c_path, c_rows) in zip(m_axs.flatten(),\n                    sample_df.groupby(['path'])):\n    c_img_arr = pydicom.read_file(c_path).pixel_array\n    # overlay\n    c_img = plt.cm.gray(c_img_arr)\n    c_img += 0.25*plt.cm.hot(prob_image/prob_image.max())\n    c_img = np.clip(c_img, 0, 1)\n    c_ax.imshow(c_img)\n    \n    c_ax.set_title('{class}'.format(**c_rows.iloc[0,:]))\n    for i, (_, c_row) in enumerate(c_rows.dropna().iterrows()):\n        c_ax.plot(c_row['x'], c_row['y'], 's', label='{class}'.format(**c_row))\n        c_ax.add_patch(Rectangle(xy=(c_row['x'], c_row['y']),\n                                width=c_row['width'],\n                                height=c_row['height'], \n                                 alpha = 0.5,\n                                fill=False))\n        if i==0: c_ax.legend()\nfig.savefig('overview.png', figdpi = 600)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Save the preprocessed results\n\nWe can use the preprocessed results with the appropriate DICOM tags to make model training step easier"},{"metadata":{"trusted":true},"cell_type":"code","source":"image_bbox_df.to_csv('image_bbox_full.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class generator(keras.utils.Sequence):\n    \n    def __init__(self, folder, filenames, pneumonia_locations=None, batch_size=32, image_size=256, shuffle=True, augment=False, predict=False):\n        self.folder = folder\n        self.filenames = filenames\n        self.pneumonia_locations = pneumonia_locations\n        self.batch_size = batch_size\n        self.image_size = image_size\n        self.shuffle = shuffle\n        self.augment = augment\n        self.predict = predict\n        self.on_epoch_end()\n        \n    def __load__(self, filename):\n        # load dicom file as numpy array\n        img = pydicom.dcmread(os.path.join(self.folder, filename)).pixel_array\n        # create empty mask\n        msk = np.zeros(img.shape)\n        # get filename without extension\n        filename = filename.split('.')[0]\n        # if image contains pneumonia\n        if filename in self.pneumonia_locations:\n            # loop through pneumonia\n            for location in self.pneumonia_locations[filename]:\n                # add 1's at the location of the pneumonia\n                x, y, w, h = location\n                msk[y:y+h, x:x+w] = 1\n        # resize both image and mask\n        img = resize(img, (self.image_size, self.image_size), mode='reflect')\n        msk = resize(msk, (self.image_size, self.image_size), mode='reflect') > 0.5\n        # if augment then horizontal flip half the time\n        if self.augment and random.random() > 0.5:\n            img = np.fliplr(img)\n            msk = np.fliplr(msk)\n        # add trailing channel dimension\n        img = np.expand_dims(img, -1)\n        msk = np.expand_dims(msk, -1)\n        return img, msk\n    \n    def __loadpredict__(self, filename):\n        # load dicom file as numpy array\n        img = pydicom.dcmread(os.path.join(self.folder, filename)).pixel_array\n        # resize image\n        img = resize(img, (self.image_size, self.image_size), mode='reflect')\n        # add trailing channel dimension\n        img = np.expand_dims(img, -1)\n        return img\n        \n    def __getitem__(self, index):\n        # select batch\n        filenames = self.filenames[index*self.batch_size:(index+1)*self.batch_size]\n        # predict mode: return images and filenames\n        if self.predict:\n            # load files\n            imgs = [self.__loadpredict__(filename) for filename in filenames]\n            # create numpy batch\n            imgs = np.array(imgs)\n            return imgs, filenames\n        # train mode: return images and masks\n        else:\n            # load files\n            items = [self.__load__(filename) for filename in filenames]\n            # unzip images and masks\n            imgs, msks = zip(*items)\n            # create numpy batch\n            imgs = np.array(imgs)\n            msks = np.array(msks)\n            return imgs, msks\n        \n    def on_epoch_end(self):\n        if self.shuffle:\n            random.shuffle(self.filenames)\n        \n    def __len__(self):\n        if self.predict:\n            # return everything\n            return int(np.ceil(len(self.filenames) / self.batch_size))\n        else:\n            # return full batches only\n            return int(len(self.filenames) / self.batch_size)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def create_downsample(channels, inputs):\n    x = keras.layers.BatchNormalization(momentum=0.9)(inputs)\n    x = keras.layers.LeakyReLU(0)(x)\n    x = keras.layers.Conv2D(channels, 1, padding='same', use_bias=False)(x)\n    x = keras.layers.MaxPool2D(2)(x)\n    return x\n\ndef create_resblock(channels, inputs):\n    x = keras.layers.BatchNormalization(momentum=0.9)(inputs)\n    x = keras.layers.LeakyReLU(0)(x)\n    x = keras.layers.Conv2D(channels, 3, padding='same', use_bias=False)(x)\n    x = keras.layers.BatchNormalization(momentum=0.9)(x)\n    x = keras.layers.LeakyReLU(0)(x)\n    x = keras.layers.Conv2D(channels, 3, padding='same', use_bias=False)(x)\n    return keras.layers.add([x, inputs])\n\ndef create_network(input_size, channels, n_blocks=2, depth=4):\n    # input\n    inputs = keras.Input(shape=(input_size, input_size, 1))\n    x = keras.layers.Conv2D(channels, 3, padding='same', use_bias=False)(inputs)\n    # residual blocks\n    for d in range(depth):\n        channels = channels * 2\n        x = create_downsample(channels, x)\n        for b in range(n_blocks):\n            x = create_resblock(channels, x)\n    # output\n    x = keras.layers.BatchNormalization(momentum=0.9)(x)\n    x = keras.layers.LeakyReLU(0)(x)\n    x = keras.layers.Conv2D(1, 1, activation='sigmoid')(x)\n    outputs = keras.layers.UpSampling2D(2**depth)(x)\n    model = keras.Model(inputs=inputs, outputs=outputs)\n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# define iou or jaccard loss function\ndef iou_loss(y_true, y_pred):\n    y_true = tf.reshape(y_true, [-1])\n    y_pred = tf.reshape(y_pred, [-1])\n    intersection = tf.reduce_sum(y_true * y_pred)\n    score = (intersection + 1.) / (tf.reduce_sum(y_true) + tf.reduce_sum(y_pred) - intersection + 1.)\n    return 1 - score\n\n# combine bce loss and iou loss\ndef iou_bce_loss(y_true, y_pred):\n    return 0.5 * keras.losses.binary_crossentropy(y_true, y_pred) + 0.5 * iou_loss(y_true, y_pred)\n\n# mean iou as a metric\ndef mean_iou(y_true, y_pred):\n    y_pred = tf.round(y_pred)\n    intersect = tf.reduce_sum(y_true * y_pred, axis=[1, 2, 3])\n    union = tf.reduce_sum(y_true, axis=[1, 2, 3]) + tf.reduce_sum(y_pred, axis=[1, 2, 3])\n    smooth = tf.ones(tf.shape(intersect))\n    return tf.reduce_mean((intersect + smooth) / (union - intersect + smooth))\n\n# create network and compiler\nmodel = create_network(input_size=256, channels=32, n_blocks=2, depth=4)\nmodel.compile(optimizer='adam',\n              loss=iou_bce_loss,\n              metrics=['accuracy', mean_iou])\n\n# cosine learning rate annealing\ndef cosine_annealing(x):\n    lr = 0.001\n    epochs = 10\n    return lr*(np.cos(np.pi*x/epochs)+1.)/2\nlearning_rate = tf.keras.callbacks.LearningRateScheduler(cosine_annealing)\n\n# create train and validation generators\nfolder = '../input/rsna-pneumonia-detection-challenge/stage_2_train_images'\ntrain_gen = generator(folder, train_filenames, pneumonia_locations, batch_size=32, image_size=256, shuffle=True, augment=True, predict=False)\nvalid_gen = generator(folder, valid_filenames, pneumonia_locations, batch_size=32, image_size=256, shuffle=False, predict=False)\n\nhistory = model.fit_generator(train_gen, validation_data=valid_gen, callbacks=[learning_rate], epochs=10, workers=4, use_multiprocessing=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(12,4))\nplt.subplot(131)\nplt.plot(history.epoch, history.history[\"loss\"], label=\"Train loss\")\nplt.plot(history.epoch, history.history[\"val_loss\"], label=\"Valid loss\")\nplt.legend()\nplt.subplot(132)\nplt.plot(history.epoch, history.history[\"acc\"], label=\"Train accuracy\")\nplt.plot(history.epoch, history.history[\"val_acc\"], label=\"Valid accuracy\")\nplt.legend()\nplt.subplot(133)\nplt.plot(history.epoch, history.history[\"mean_iou\"], label=\"Train iou\")\nplt.plot(history.epoch, history.history[\"val_mean_iou\"], label=\"Valid iou\")\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for imgs, msks in valid_gen:\n    # predict batch of images\n    preds = model.predict(imgs)\n    # create figure\n    f, axarr = plt.subplots(4, 8, figsize=(20,15))\n    axarr = axarr.ravel()\n    axidx = 0\n    # loop through batch\n    for img, msk, pred in zip(imgs, msks, preds):\n        # plot image\n        axarr[axidx].imshow(img[:, :, 0])\n        # threshold true mask\n        comp = msk[:, :, 0] > 0.5\n        # apply connected components\n        comp = measure.label(comp)\n        # apply bounding boxes\n        predictionString = ''\n        for region in measure.regionprops(comp):\n            # retrieve x, y, height and width\n            y, x, y2, x2 = region.bbox\n            height = y2 - y\n            width = x2 - x\n            axarr[axidx].add_patch(patches.Rectangle((x,y),width,height,linewidth=2,edgecolor='b',facecolor='none'))\n        # threshold predicted mask\n        comp = pred[:, :, 0] > 0.5\n        # apply connected components\n        comp = measure.label(comp)\n        # apply bounding boxes\n        predictionString = ''\n        for region in measure.regionprops(comp):\n            # retrieve x, y, height and width\n            y, x, y2, x2 = region.bbox\n            height = y2 - y\n            width = x2 - x\n            axarr[axidx].add_patch(patches.Rectangle((x,y),width,height,linewidth=2,edgecolor='r',facecolor='none'))\n        axidx += 1\n    plt.show()\n    # only plot one batch\n    break","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# load and shuffle filenames\nfolder = '../input/rsna-pneumonia-detection-challenge/stage_2_test_images'\ntest_filenames = os.listdir(folder)\nprint('n test samples:', len(test_filenames))\n\n# create test generator with predict flag set to True\ntest_gen = generator(folder, test_filenames, None, batch_size=25, image_size=256, shuffle=False, predict=True)\n\n# create submission dictionary\nsubmission_dict = {}\n# loop through testset\nfor imgs, filenames in test_gen:\n    # predict batch of images\n    preds = model.predict(imgs)\n    # loop through batch\n    for pred, filename in zip(preds, filenames):\n        # resize predicted mask\n        pred = resize(pred, (1024, 1024), mode='reflect')\n        # threshold predicted mask\n        comp = pred[:, :, 0] > 0.5\n        # apply connected components\n        comp = measure.label(comp)\n        # apply bounding boxes\n        predictionString = ''\n        for region in measure.regionprops(comp):\n            # retrieve x, y, height and width\n            y, x, y2, x2 = region.bbox\n            height = y2 - y\n            width = x2 - x\n            # proxy for confidence score\n            conf = np.mean(pred[y:y+height, x:x+width])\n            # add to predictionString\n            predictionString += str(conf) + ' ' + str(x) + ' ' + str(y) + ' ' + str(width) + ' ' + str(height) + ' '\n        # add filename and predictionString to dictionary\n        filename = filename.split('.')[0]\n        submission_dict[filename] = predictionString\n    # stop if we've got them all\n    if len(submission_dict) >= len(test_filenames):\n        break\n\n# save dictionary as csv file\nsub = pd.DataFrame.from_dict(submission_dict,orient='index')\nsub.index.names = ['patientId']\nsub.columns = ['PredictionString']\nsub.to_csv('submission.csv')","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}