{"cells":[{"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 nodules.\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 nodule locations\n\nTable contains [filename : nodule location] pairs per row. \n* If a filename contains multiple nodules, the table contains multiple rows with the same filename but different nodule locations. \n* If a filename contains no nodules it contains a single row with an empty nodule 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 nodule locations in that filename as value. \n* If a filename is not present in the dictionary it means that it contains no nodules."},{"metadata":{"trusted":true,"_uuid":"e08496a85ef9b0823595c3745d2677c6e84b6a3a","collapsed":true},"cell_type":"code","source":"# empty dictionary\nnodule_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        nodule = rows[5]\n        # if row contains a nodule add label to dictionary\n        # which contains a list of nodule locations per filename\n        if nodule == '1':\n            # convert string to float to int\n            location = [int(float(i)) for i in location]\n            # save nodule location in dictionary\n            if filename in nodule_locations:\n                nodule_locations[filename].append(location)\n            else:\n                nodule_locations[filename] = [location]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5a4fa6da9833fd7cbd476d19584ba35695b38ddb"},"cell_type":"markdown","source":"# Load filenames"},{"metadata":{"trusted":true,"_uuid":"ccd0b0d52cafd125558ed5560a9cc8fa15760bc5"},"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":{"_uuid":"bd0633867d9d32180eb776703205767e8e891c50"},"cell_type":"markdown","source":"# Exploration"},{"metadata":{"trusted":true,"_uuid":"daa7156380489227e510e8ef086b55c58b4b3ad8"},"cell_type":"code","source":"print('Total train images:',len(filenames))\nprint('Images with nodule:', len(nodule_locations))\n\nns = [len(value) for value in nodule_locations.values()]\nplt.figure()\nplt.hist(ns)\nplt.xlabel('Nodules per image')\nplt.xticks(range(1, np.max(ns)+1))\nplt.show()\n\nheatmap = np.zeros((1024, 1024))\nws = []\nhs = []\nfor values in nodule_locations.values():\n    for value in values:\n        x, y, w, h = value\n        heatmap[y:y+h, x:x+w] += 1\n        ws.append(w)\n        hs.append(h)\nplt.figure()\nplt.title('Nodule location heatmap')\nplt.imshow(heatmap)\nplt.figure()\nplt.title('Nodule height lengths')\nplt.hist(hs, bins=np.linspace(0,1000,50))\nplt.show()\nplt.figure()\nplt.title('Nodule width lengths')\nplt.hist(ws, bins=np.linspace(0,1000,50))\nplt.show()\nprint('Minimum nodule height:', np.min(hs))\nprint('Minimum nodule width: ', np.min(ws))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1851594939b259a95dd906639a1ef87fe8400f83"},"cell_type":"code","source":"from sklearn.cluster import KMeans\nkmeans = KMeans(n_clusters=2).fit(np.array([ws, hs]).T)\ncenters = kmeans.cluster_centers_\nplt.figure()\nplt.scatter(ws, hs, marker='.')\nplt.xlabel('width')\nplt.ylabel('height')\nfor center in centers:\n    print(center)\n    plt.scatter(center[0], center[1], c='red')\nplt.show()","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, nodule_locations=None, batch_size=32, image_size=512, shuffle=True, predict=False, augment = False):\n        self.folder = folder\n        self.filenames = filenames\n        self.nodule_locations = nodule_locations\n        self.batch_size = batch_size\n        self.image_size = image_size\n        self.augment = augment\n        self.shuffle = shuffle\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 nodules\n        if filename in nodule_locations:\n            # loop through nodules\n            for location in nodule_locations[filename]:\n                # add 1's at the location of the nodule\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\n## BatchNormormalization: Fine-Tune your Booster\n\n* Batch Normalization is work like normalize the sample data so we can train for better accuracy and result. It is same like we can do in Machine learning.\n* We can even apply this normalization to the input of neural networks. It fastens up training as in linear regression. But since the 2nd layer changes this distribution, the \n\n**Defination : Batch Normalization**\n* To reduce this **problem of internal covariate shift,  Batch Normalization adds Normalization “layer” between each layers.** \n* An important thing to note here is that **normalization has to be done separately for each dimension (input neuron), over the ‘mini-batches’**, and not altogether with all dimensions. Hence the name ‘batch’ normalization.\n![](https://cdn-images-1.medium.com/max/720/1*WRio7MD4JDeLww-CyrxEbg.png)    \n\n* Due to this normalization “layers” between each fully connected layers, the range of input distribution of each layer stays the same, no matter the changes in the previous layer. Given x inputs from k-th neuron.\n\n* Normalization brings all the inputs centered around 0. This way, there is not much change in each layer input. So, layers in the network can learn from the back-propagation simultaneously, without waiting for the previous layer to learn. This fastens up the training of networks.\n\n![](https://cdn-images-1.medium.com/max/720/1*c9mS1eNNJQqIFKvu8_Ji7g.png)\n\n* We apply Batch Normalization to the best-performing ImageNet classification network, and show that we can match its performance using only 7% of the training steps, and can further exceed its accuracy by a substantial margin. — [Original BatchNorm Paper](https://arxiv.org/pdf/1502.03167.pdf)\n\n* But, sometimes, it is difficult to keep track of all the mini-batch mean and variances. In such cases, exponential weighted “moving average” can be used to update population mean and variance:\n\n**Momentum**\n![](https://cdn-images-1.medium.com/max/720/1*nePTWq-SqDKOoTyl2ztR-A.png)\n\n* **Here, α is the momentum, γ = σ and β = mean**\n * Momentum is the importance given to the previous moving average, when calculating the population average for inference. If you can’t understand what momentum is, It’s nothing but the smoothing that we can adjust in Tensorboard. Momentum is the “lag” in learning mean and variance, so that noise due to mini-batch can be ignored.\n* **Actual(light) and lagged(bold) values with momentum 0.99 and 0.75**\n* **By default, momentum would be set a high value about 0.99**, meaning high lag and slow learning. When batch sizes are small, the no. of steps run will be more. So high momentum will result in slow but steady learning (more lag) of the moving mean. So, in this case, it is helpful.\n* But when the batch size is bigger, as I have used, i.e 5K images (out of 50K) in single step, the number of steps is less. Also, the statistics of mini-batch are mostly same as that of the population. At these times, momentum has to be less, so that the mean and variance are updated quickly. Hence a ground rule is that:\n>**Small batch size => High Momentum (0.9–0.99) **  \n>**Big batch size => Low Momentum (0.6–0.85)**\n\n**Calculation of moving mean and variance**\n* Here **α** is the **“momentum” given to previous moving statistic, around 0.9.** And those with **B subscript are mini-batch mean and mini-batch variance.** This is the implementation found in most libraries, where the momentum can be set manually.\n* An important thing to note here is that the **moving mean and moving variance** are calculated at training time, with training dataset and not at testing time."},{"metadata":{"trusted":true,"_uuid":"9fc2b108689637a6037b48ebab3f7659b8704bf9","collapsed":true},"cell_type":"code","source":"def create_downsample(channels, inputs):\n    x = keras.layers.BatchNormalization()(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()(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()(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=5):\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()(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","scrolled":false},"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=.01),loss=keras.losses.binary_crossentropy,metrics=['accuracy', mean_iou])\n\n# create train and validation generators\nfolder = '../input/stage_1_train_images'\ntrain_gen = generator(folder, train_filenames, nodule_locations, batch_size=32, image_size=256, shuffle=True, augment=True, predict=False)\nvalid_gen = generator(folder, valid_filenames, nodule_locations, batch_size=32, image_size=256, shuffle=False, predict=False)\n\nhistory = model.fit_generator(train_gen, validation_data=valid_gen, epochs=3, shuffle=True, verbose=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"257c31943ff5226d7d60135f07e3168d40bd4e85"},"cell_type":"code","source":"model.save(\"model.h5\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3666ba4cac9ed2c3029b824af220404bfcc16f23"},"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":"25dc7ad848190bb37349a80f96ed2bdb5c3821b0"},"cell_type":"markdown","source":"# Predict test images"},{"metadata":{"trusted":true,"_uuid":"2c9277e4ec9f12712dd690002c540b396278c504","scrolled":false},"cell_type":"code","source":"# load and shuffle filenames\nfolder = '../input/stage_1_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('submission1.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"a953a21bd4cb1f7a8d7d1bb875272f81cd080be9"},"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}