{"cells":[{"metadata":{"_uuid":"15a9d472fcd107b6c52619b6dcee42af867d31ef"},"cell_type":"markdown","source":"Based on Jonne's notebook (see fork link), using [YichengChen's evaluation metric function](https://www.kaggle.com/chenyc15/mean-average-precision-metric)\n"},{"metadata":{"_uuid":"d320f90f432c331afea02ef66fcddb59448ea200"},"cell_type":"markdown","source":"# Approach\n\n* Firstly a convolutional neural network is used to segment the image, using the bounding boxes directly as a mask. \n* Secondly connected components is used to separate multiple areas of predicted pneumonia.\n* Finally a bounding box is simply drawn around every connected component.\n\n# Network\n\n* The network consists of a number of residual blocks with convolutions and downsampling blocks with max pooling.\n* At the end of the network a single upsampling layer converts the output to the same shape as the input.\n\nAs the input to the network is 256 by 256 (instead of the original 1024 by 1024) and the network downsamples a number of times without any meaningful upsampling (the final upsampling is just to match in 256 by 256 mask) the final prediction is very crude. If the network downsamples 4 times the final bounding boxes can only change with at least 16 pixels."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"collapsed":true},"cell_type":"code","source":"import os\nimport csv\nimport random\nimport pydicom\nimport numpy as np\nimport pandas as pd\nfrom skimage import measure\nfrom skimage.transform import resize\n\nimport tensorflow as tf\nfrom tensorflow import keras\n\nfrom matplotlib import pyplot as plt","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8d58cefae21c951e077c02c1a989d020f18465cc"},"cell_type":"markdown","source":"# Load pneumonia locations\n\nTable contains [filename : pneumonia location] pairs per row. \n* If a filename contains multiple pneumonia, the table contains multiple rows with the same filename but different pneumonia locations. \n* If a filename contains no pneumonia it contains a single row with an empty pneumonia location.\n\nThe code below loads the table and transforms it into a dictionary. \n* The dictionary uses the filename as key and a list of pneumonia locations in that filename as value. \n* If a filename is not present in the dictionary it means that it contains no pneumonia."},{"metadata":{"trusted":true,"_uuid":"e08496a85ef9b0823595c3745d2677c6e84b6a3a","collapsed":true},"cell_type":"code","source":"# empty dictionary\npneumonia_locations = {}\n# load table\nwith open(os.path.join('../input/stage_1_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":"5a4fa6da9833fd7cbd476d19584ba35695b38ddb"},"cell_type":"markdown","source":"# Load filenames"},{"metadata":{"trusted":true,"_uuid":"ccd0b0d52cafd125558ed5560a9cc8fa15760bc5","collapsed":true},"cell_type":"code","source":"# load and shuffle filenames\nfolder = '../input/stage_1_train_images'\nfilenames = os.listdir(folder)\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":{"trusted":true,"collapsed":true,"_uuid":"276b80f59fa9acdfd9307d74cee5f705cd6aa5b6"},"cell_type":"markdown","source":" # Data generator\n\nThe dataset is too large to fit into memory, so we need to create a generator that loads data on the fly.\n\n* The generator takes in some filenames, batch_size and other parameters.\n\n* The generator outputs a random batch of numpy images and numpy masks.\n    "},{"metadata":{"trusted":true,"_uuid":"86b3f780a03cddda78c6adfde461d6ff8dad5672","collapsed":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 pneumonia_locations:\n            # loop through pneumonia\n            for location in 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":{"_uuid":"ebc622a4b406354cc4ef28801eab72346a724d8b"},"cell_type":"markdown","source":"# Network"},{"metadata":{"trusted":true,"_uuid":"9fc2b108689637a6037b48ebab3f7659b8704bf9","collapsed":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":{"_uuid":"bee3bc9363e9c1829eadf17da7a67fbd1a6a369e"},"cell_type":"markdown","source":"# Train network\n"},{"metadata":{"trusted":true,"_uuid":"4369be30f61440eb6858d57829fa541c4ee893bf","collapsed":true},"cell_type":"code","source":"# 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=keras.optimizers.Adam(lr=.02),\n              loss=keras.losses.binary_crossentropy,\n              metrics=['accuracy', mean_iou])\n\n# create train and validation generators\nfolder = '../input/stage_1_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, \n                              epochs=8, shuffle=True, verbose=2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3666ba4cac9ed2c3029b824af220404bfcc16f23","collapsed":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":{"_uuid":"df3c6c5d2402b3cfe51d3e78bf1d8ca2566f8c89"},"cell_type":"markdown","source":"# Evaluation metric functions, from Yicheng Chen's kernel"},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"6044a5636ebda68f5ef29511bec7a9f505ca7211"},"cell_type":"code","source":"# helper function to calculate IoU\ndef iou(box1, box2):\n    x11, y11, w1, h1 = box1\n    x21, y21, w2, h2 = box2\n    assert w1 * h1 > 0\n    assert w2 * h2 > 0\n    x12, y12 = x11 + w1, y11 + h1\n    x22, y22 = x21 + w2, y21 + h2\n\n    area1, area2 = w1 * h1, w2 * h2\n    xi1, yi1, xi2, yi2 = max([x11, x21]), max([y11, y21]), min([x12, x22]), min([y12, y22])\n    \n    if xi2 <= xi1 or yi2 <= yi1:\n        return 0\n    else:\n        intersect = (xi2-xi1) * (yi2-yi1)\n        union = area1 + area2 - intersect\n        return intersect / union\n    \ndef map_iou(boxes_true, boxes_pred, scores, thresholds = [0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75]):\n    \"\"\"\n    Mean average precision at differnet intersection over union (IoU) threshold\n    \n    input:\n        boxes_true: Mx4 numpy array of ground true bounding boxes of one image. \n                    bbox format: (x1, y1, w, h)\n        boxes_pred: Nx4 numpy array of predicted bounding boxes of one image. \n                    bbox format: (x1, y1, w, h)\n        scores:     length N numpy array of scores associated with predicted bboxes\n        thresholds: IoU shresholds to evaluate mean average precision on\n    output: \n        map: mean average precision of the image\n    \"\"\"\n    \n    # According to the introduction, images with no ground truth bboxes will not be \n    # included in the map score unless there is a false positive detection (?)\n        \n    # return None if both are empty, don't count the image in final evaluation (?)\n    if len(boxes_true) == 0 and len(boxes_pred) == 0:\n        return None\n    \n    assert boxes_true.shape[1] == 4 or boxes_pred.shape[1] == 4, \"boxes should be 2D arrays with shape[1]=4\"\n    if len(boxes_pred):\n        assert len(scores) == len(boxes_pred), \"boxes_pred and scores should be same length\"\n        # sort boxes_pred by scores in decreasing order\n        boxes_pred = boxes_pred[np.argsort(scores)[::-1], :]\n    \n    map_total = 0\n    \n    # loop over thresholds\n    for t in thresholds:\n        matched_bt = set()\n        tp, fn = 0, 0\n        for i, bt in enumerate(boxes_true):\n            matched = False\n            for j, bp in enumerate(boxes_pred):\n                miou = iou(bt, bp)\n                if miou >= t and not matched and j not in matched_bt:\n                    matched = True\n                    tp += 1 # bt is matched for the first time, count as TP\n                    matched_bt.add(j)\n            if not matched:\n                fn += 1 # bt has no match, count as FN\n                \n        fp = len(boxes_pred) - len(matched_bt) # FP is the bp that not matched to any bt\n        m = tp / (tp + fn + fp)\n        map_total += m\n    \n    return map_total / len(thresholds)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"25dc7ad848190bb37349a80f96ed2bdb5c3821b0"},"cell_type":"markdown","source":"# Predict validaiton images"},{"metadata":{"trusted":true,"_uuid":"2c9277e4ec9f12712dd690002c540b396278c504","scrolled":false,"collapsed":true},"cell_type":"code","source":"prob_thresholds = [0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5]  # Prediction thresholds to consider a pixel positive\nnthresh = len(prob_thresholds)\n\n# load and shuffle filenames\nfolder = '../input/stage_1_train_images'\ntest_filenames = valid_filenames\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# loop through validation set\ncount = 0\n\nns = nthresh*[0]     # Keep separate statistics for each threshold\nnfps = nthresh*[0]\nntps = nthresh*[0]\noverall_maps = nthresh*[0.]\n\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        count = count + 1\n        maxpred = np.max(pred)\n        # resize predicted mask\n        pred = resize(pred, (1024, 1024), mode='reflect')\n        boxes_preds = []  # Keep lists of boxes and scores for each candidate threshold\n        scoress = []\n        for thresh in prob_thresholds:\n            # threshold predicted mask\n            comp = pred[:, :, 0] > thresh\n            # apply connected components\n            comp = measure.label(comp)\n            # apply bounding boxes\n            boxes_pred = np.empty((0,4),int)\n            scores = np.empty((0))\n            for region in measure.regionprops(comp):\n                # retrieve x, y, height and width\n                y, x, y2, x2 = region.bbox\n                boxes_pred = np.append(boxes_pred, [[x, y, x2-x, y2-y]], axis=0)\n                # proxy for confidence score\n                conf = np.mean(pred[y:y2, x:x2])\n                scores = np.append( scores, conf )\n            boxes_preds = boxes_preds + [boxes_pred]\n            scoress = scoress + [scores]\n        boxes_true = np.empty((0,4),int)\n        fn = filename.split('.')[0]\n        # if image contains pneumonia\n        if fn in pneumonia_locations:\n            # loop through pneumonia\n            for location in pneumonia_locations[fn]:\n                x, y, w, h = location\n                boxes_true = np.append(boxes_true, [[x, y, w, h]], axis=0)\n        for i in range(nthresh):\n            if ( boxes_true.shape[0]==0 and boxes_preds[i].shape[0]>0 ):  # false positive\n                ns[i] = ns[i] + 1       # increment denominator but add nothing to numerator\n                nfps[i] = nfps[i] + 1   # track number of false positive cases, for curiosity\n            elif ( boxes_true.shape[0]>0 ):  # actual positive\n                ns[i] = ns[i] + 1       # increment denominator & add contribution to numerator\n                contrib = map_iou( boxes_true, boxes_preds[i], scoress[i] ) \n                overall_maps[i] = overall_maps[i] + contrib\n                if ( boxes_preds[i].shape[0]>0 ):  # true positive\n                    ntps[i] = ntps[i] + 1  # track number of true positive cases, for curiosity\n\n    # stop if we've got them all\n    if count >= len(test_filenames):\n        break\n\nfor i, thresh in enumerate(prob_thresholds):\n    print( \"\\nProbability threshold \", thresh )\n    overall_maps[i] = overall_maps[i] / (ns[i]+1e-7)\n    print( \"False positive cases:  \", nfps[i] )\n    print( \"True positive cases: \", ntps[i] )\n    print( \"Overall evaluation score: \", overall_maps[i] )\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"9586db3beaa6033e4813294aef5d2923c02d24cd"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}