{"cells":[{"metadata":{"_cell_guid":"b110dd56-4ccf-4cf8-8d0b-663763d7ac3e","_uuid":"e3cfa2c10efaa058c7b990faa1faf5b3702cc8e0"},"cell_type":"markdown","source":"*March 2018*\n\n**Introduction**\n\n   Problem case:\n\nTL;DR: In a nutshell, create an algorithm to automate nucleus detection in divergent images  to advance medical discovery and unlock faster cures. [(Read more here)](https://datasciencebowl.com)\n\n\n   Solution:\n\nThe solution is to build a CNN in the shape of a U, known as a U-Net. The U-Net, according to the [authors at the University of Freiburg](https://lmb.informatik.uni-freiburg.de/people/ronneber/u-net/), is best suited for convolutional neural networks for Biomedical Image Segmentation.\n\nThe CNN will be built on the training data and applied to the test data.\n\n**U-Net architecture flow:**\n\n![U-Net architecture flow:](https://lmb.informatik.uni-freiburg.de/Publications/2015/RFB15a/u-net-architecture.png)\n\nThe main idea of this 'fully connected network' is to supplement a usual contracting network by successive layers, where pooling operators are replaced by upsampling operators. Hence, these layers (the right side of the 'U') **increase** the resolution of the output. \n\n*Excerpt* from [original paper:](https://arxiv.org/pdf/1505.04597.pdf) \n\n**Network Architecture**\n\nThe network architecture consists of a contracting path (left side) and an expansive path (right side). The contracting path follows the typical architecture of a convolutional network. It consists of the repeated application of two 3x3 convolutions (unpadded convolutions), each followed by a rectified linear unit (ReLU) and a 2x2 max pooling operation with stride 2 for downsampling. At each downsampling step we double the number of feature channels. Every step in the expansive path consists of an upsampling of the feature map followed by a 2x2 convolution (“up-convolution”) that halves the number of feature channels, a concatenation with the correspondingly cropped feature map from the contracting path, and two 3x3 convolutions, each followed by a ReLU. The cropping is necessary due to the loss of border pixels in every convolution. At the final layer a 1x1 convolution is used to map each 64- component feature vector to the desired number of classes. In total the network has 23 convolutional layers.\n\n*End excerpt.*\n\n"},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","collapsed":true,"trusted":false},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport cv2\nimport os\nimport numpy as np\n\n\ndef make_df(train_path, test_path, img_size):\n    train_ids = next(os.walk(train_path))[1]\n    test_ids = next(os.walk(test_path))[1]\n    X_train = np.zeros((len(train_ids), img_size, img_size, 3), dtype=np.uint8)\n    Y_train = np.zeros((len(train_ids), img_size, img_size, 1), dtype=np.bool)\n    for i, id_ in enumerate(train_ids):\n        path = train_path + id_\n        img = cv2.imread(path + '/images/' + id_ + '.png')\n        img = cv2.resize(img, (img_size, img_size))\n        X_train[i] = img\n        mask = np.zeros((img_size, img_size, 1), dtype=np.bool)\n        for mask_file in next(os.walk(path + '/masks/'))[2]:\n            mask_ = cv2.imread(path + '/masks/' + mask_file, 0)\n            mask_ = cv2.resize(mask_, (img_size, img_size))\n            mask_ = mask_[:, :, np.newaxis]\n            mask = np.maximum(mask, mask_)\n        Y_train[i] = mask\n    X_test = np.zeros((len(test_ids), img_size, img_size, 3), dtype=np.uint8)\n    sizes_test = []\n    for i, id_ in enumerate(test_ids):\n        path = test_path + id_\n        img = cv2.imread(path + '/images/' + id_ + '.png')\n        sizes_test.append([img.shape[0], img.shape[1]])\n        img = cv2.resize(img, (img_size, img_size))\n        X_test[i] = img\n\n    return X_train, Y_train, X_test, sizes_test","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","collapsed":true},"cell_type":"markdown","source":"Define the U-Net model\n\nIn this next section, we define the U-Net model, 9 parts in all. Each part in the architecture consists of 4 operations:\n\n    1. Convolve\n    2. Dropout\n    3. Convolve\n    4. Maxpooling\n \n If you do not what there are about please check out Stanfords fabled [Deep Learning specialization](https://www.coursera.org/specializations/deep-learning) on Coursera, run by professor Andrew Ng, the unicorn of AI. \n Now when it comes to the activation function, there are a whole host that you can test for increased accuracy, like elu, relu, selu and tanh."},{"metadata":{"_cell_guid":"fb1b3149-eac2-4ef7-a25f-0487e89ccf13","_uuid":"c87f786d0f898f372283cebb7ba18fd4866a81b3","collapsed":true,"trusted":false},"cell_type":"code","source":"from keras.layers.convolutional import Conv2D, Conv2DTranspose\nfrom keras.layers.pooling import MaxPooling2D\nfrom keras.layers.merge import concatenate\nfrom keras.models import Model\nfrom keras.layers import Input\nfrom keras.layers.core import Dropout, Lambda\n\n\ndef Unet(img_size):\n    inputs = Input((img_size, img_size, 3))\n    s = Lambda(lambda x: x / 255)(inputs)\n\n    c1 = Conv2D(16, (3, 3), activation='tanh', kernel_initializer='he_normal', padding='valid')(s)\n    c1 = Dropout(0.1)(c1)\n    c1 = Conv2D(16, (3, 3), activation='tanh', kernel_initializer='he_normal', padding='valid')(c1)\n    p1 = MaxPooling2D((2, 2))(c1)\n    print('c1')\n    c2 = Conv2D(32, (3, 3), activation='tanh', kernel_initializer='he_normal', padding='valid')(p1)\n    c2 = Dropout(0.1)(c2)\n    c2 = Conv2D(32, (3, 3), activation='tanh', kernel_initializer='he_normal', padding='valid')(c2)\n    p2 = MaxPooling2D((2, 2))(c2)\n    print('c2')\n    c3 = Conv2D(64, (3, 3), activation='tanh', kernel_initializer='he_normal', padding='valid')(p2)\n    c3 = Dropout(0.2)(c3)\n    c3 = Conv2D(64, (3, 3), activation='tanh', kernel_initializer='he_normal', padding='valid')(c3)\n    p3 = MaxPooling2D((2, 2))(c3)\n    print('c3')\n    c4 = Conv2D(128, (3, 3), activation='tanh', kernel_initializer='he_normal', padding='valid')(p3)\n    c4 = Dropout(0.2)(c4)\n    c4 = Conv2D(128, (3, 3), activation='tanh', kernel_initializer='he_normal', padding='valid')(c4)\n    p4 = MaxPooling2D(pool_size=(2, 2))(c4)\n    print('c4')\n    c5 = Conv2D(256, (3, 3), activation='tanh', kernel_initializer='he_normal', padding='valid')(p4)\n    c5 = Dropout(0.3)(c5)\n    c5 = Conv2D(256, (3, 3), activation='tanh', kernel_initializer='he_normal', padding='valid')(c5)\n    print('c5')\n    u6 = Conv2DTranspose(128, (2, 2), strides=(2, 2), padding='same')(c5)\n    u6 = concatenate([u6, c4])\n    c6 = Conv2D(128, (3, 3), activation='tanh', kernel_initializer='he_normal', padding='valid')(u6)\n    c6 = Dropout(0.2)(c6)\n    c6 = Conv2D(128, (3, 3), activation='tanh', kernel_initializer='he_normal', padding='valid')(c6)\n    print('c6')\n    u7 = Conv2DTranspose(64, (2, 2), strides=(2, 2), padding='same')(c6)\n    u7 = concatenate([u7, c3])\n    c7 = Conv2D(64, (3, 3), activation='tanh', kernel_initializer='he_normal', padding='valid')(u7)\n    c7 = Dropout(0.2)(c7)\n    c7 = Conv2D(64, (3, 3), activation='tanh', kernel_initializer='he_normal', padding='valid')(c7)\n    print('c7')\n    u8 = Conv2DTranspose(32, (2, 2), strides=(2, 2), padding='same')(c7)\n    u8 = concatenate([u8, c2])\n    c8 = Conv2D(32, (3, 3), activation='tanh', kernel_initializer='he_normal', padding='valid')(u8)\n    c8 = Dropout(0.1)(c8)\n    c8 = Conv2D(32, (3, 3), activation='tanh', kernel_initializer='he_normal', padding='valid')(c8)\n    print('c8')\n    u9 = Conv2DTranspose(16, (2, 2), strides=(2, 2), padding='same')(c8)\n    u9 = concatenate([u9, c1], axis=3)\n    c9 = Conv2D(16, (3, 3), activation='tanh', kernel_initializer='he_normal', padding='valid')(u9)\n    c9 = Dropout(0.1)(c9)\n    c9 = Conv2D(16, (3, 3), activation='tanh', kernel_initializer='he_normal', padding='valid')(c9)\n    print('c9')\n    outputs = Conv2D(1, (1, 1), activation='sigmoid')(c9)\n\n    model = Model(inputs=[inputs], outputs=[outputs])\n\n    return model","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"f0ae7cd5-c63b-4f19-b252-74ef99e92f05","_uuid":"5b7f550ee669c2137122bb948d5a0116a457cf14"},"cell_type":"markdown","source":"Define your generator\n\nIn this section we are going to generate images by performing some operations on the images, like horizontally flipping them, and vertically too."},{"metadata":{"_cell_guid":"0e359665-6741-4c74-a99d-991a4ecc3ef1","_uuid":"6e71cefbd0e26eace77e5c4474a03c3ea67fbb31","collapsed":true,"trusted":false},"cell_type":"code","source":"from keras.preprocessing.image import ImageDataGenerator\n\n\ndef generator(xtr, xval, ytr, yval, batch_size):\n    data_gen_args = dict(horizontal_flip=True,\n                         vertical_flip=True,\n                         rotation_range=90.,\n                         width_shift_range=0.1,\n                         height_shift_range=0.1,\n                         zoom_range=0.1)\n    image_datagen = ImageDataGenerator(**data_gen_args)\n    mask_datagen = ImageDataGenerator(**data_gen_args)\n    image_datagen.fit(xtr, seed=10)\n    mask_datagen.fit(ytr, seed=10)\n    image_generator = image_datagen.flow(xtr, batch_size=batch_size, seed=10)\n    mask_generator = mask_datagen.flow(ytr, batch_size=batch_size, seed=10)\n    train_generator = zip(image_generator, mask_generator)\n\n    val_gen_args = dict()\n    image_datagen_val = ImageDataGenerator(**val_gen_args)\n    mask_datagen_val = ImageDataGenerator(**val_gen_args)\n    image_datagen_val.fit(xval, seed=10)\n    mask_datagen_val.fit(yval, seed=10)\n    image_generator_val = image_datagen_val.flow(xval, batch_size=batch_size, seed=10)\n    mask_generator_val = mask_datagen_val.flow(yval, batch_size=batch_size, seed=10)\n    val_generator = zip(image_generator_val, mask_generator_val)\n\n    return train_generator, val_generator","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"3271c001-65c8-4b99-8266-6a8b47b63e85","_uuid":"e10f08e5f9d661d14b30c3075207ca9d737ba702"},"cell_type":"markdown","source":"The IoU metric is an evaluation metric that tells how well your detected bounding boxes overlap the actual object on your image. Here is the basic intuition of it:\n![IoU](https://i.stack.imgur.com/JlHnn.jpg)\n\nSo later when the model is compiled, this evaluation metric will be used to determine the best model to be compiled according to the model training data. For more on this metric and it's use, check out this [paper on Fully Convolutional Networks for Semantic Segmentation.](https://people.eecs.berkeley.edu/~jonlong/long_shelhamer_fcn.pdf) \n\nA good loss funtion is based on the dice coefficient. Here is an excerpt from a [really good paper](https://arxiv.org/pdf/1711.03954.pdf) that will give you a good intuition about the dice coefficient and how it is used in a loss funtion:\n\n*begin excerpt:*\n\nThe dice coefficient is a popular and largely used cost function in segmentation problems. Considering the predicted region P and the groundtruth region G, and by denoting |P| and |G| the sum of elements in each area, the dice coefficient is twice the ratio of the intersection over the sum of areas:\nDiceCoef(P, G) = 2|P ∩ G| / ( |P| + |G| )\nA perfect segmentation result is given by a dice coefficient of 1, while a dice coefficient of 0 refers to a completely mistaken segmentation.\n\n*end excerpt.*\n\nAnother loss function is based on cross entropy. Cross-entropy loss, or log loss, measures the performance of a classification model whose output is a probability value between 0 and 1. Cross-entropy loss increases as the predicted probability diverges from the actual label. Read more about it [here.](http://ml-cheatsheet.readthedocs.io/en/latest/loss_functions.html)\n\nI am going to use the dice coefficient in a custom loss function that uses binary cross entropy too. It is aptly called the bce_dice_loss.\n\n"},{"metadata":{"_cell_guid":"1030d317-ab1d-4a3f-b1fb-fec9d907dafd","_uuid":"b4d0d5d584036a13557177f5491c60093370677a","collapsed":true,"trusted":false},"cell_type":"code","source":"import tensorflow as tf\nimport numpy as np\nfrom keras import backend as K\nfrom keras.losses import binary_crossentropy\n\n# Define IoU metric\ndef mean_iou(y_true, y_pred):\n    prec = []\n    for t in np.arange(0.5, 1.0, 0.05):\n        y_pred_ = tf.to_int32(y_pred > t)\n        score, up_opt = tf.metrics.mean_iou(y_true, y_pred_, 2, y_true)\n        K.get_session().run(tf.local_variables_initializer())\n        with tf.control_dependencies([up_opt]):\n            score = tf.identity(score)\n        prec.append(score)\n    return K.mean(K.stack(prec), axis=0)\n\ndef dice_coef(y_true, y_pred):\n    smooth = 1.\n    y_true_f = K.flatten(y_true)\n    y_pred_f = K.flatten(y_pred)\n    intersection = K.sum(y_true_f * y_pred_f)\n    return (2. * intersection + smooth) / (K.sum(y_true_f) + K.sum(y_pred_f) + smooth)\n\ndef bce_dice_loss(y_true, y_pred):\n    return 0.5 * binary_crossentropy(y_true, y_pred) - dice_coef(y_true, y_pred)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"3a7d0619-9335-4e56-9d2c-f75e16108735","_uuid":"278f6afea048862f149465d774205a699c731331"},"cell_type":"markdown","source":"RLE encoding for submission."},{"metadata":{"_cell_guid":"14e8c3ac-57bd-4212-b9ea-9f3fe1b791f8","_uuid":"841b23fa00c47d343071ddfd8257d3c7dc4e3178","collapsed":true,"trusted":false},"cell_type":"code","source":"from skimage.morphology import label\n\ndef rle_encoding(x):\n    dots = np.where(x.T.flatten() == 1)[0]\n    run_lengths = []\n    prev = -2\n    for b in dots:\n        if (b>prev+1): run_lengths.extend((b + 1, 0))\n        run_lengths[-1] += 1\n        prev = b\n    return run_lengths\n\ndef prob_to_rles(x, cutoff=0.5):\n    lab_img = label(x > cutoff)\n    for i in range(1, lab_img.max() + 1):\n        yield rle_encoding(lab_img == i)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"df819662-8dec-427d-a950-18da3d589021","_uuid":"d6cd30966f9f286274f4bc391d10201db8f72caa"},"cell_type":"markdown","source":"Now we are ready to run our model on the training data and apply it to our test data. \nFor my optimizer I am going to test using Adam, which converges *faster* than Stochastic Gradient Descent(SGD) but can generalize [*worse*.](https://arxiv.org/pdf/1705.08292.pdf)\nI'll also test SGD. Hopefully it finishes running before hardware allowance cutoff here on kaggle!"},{"metadata":{"_cell_guid":"01641fee-49eb-4fd7-9103-17fa899f17d6","_uuid":"ed21d489d4a7d3c1ccc41ed9129b0ab27d1e5787","collapsed":true,"trusted":false},"cell_type":"code","source":"import cv2\nimport os\nimport pandas as pd\nfrom sklearn.model_selection import train_test_split\n\n\nif __name__ == \"__main__\":\n    img_size = 256\n    batch_size = 32\n    train_path = '../input/stage1_train/'\n    test_path = '../input/stage2_test_final/'\n    \n    X_train, Y_train, X_test, sizes_test = make_df(train_path, test_path, img_size)\n    xtr, xval, ytr, yval = train_test_split(X_train, Y_train, test_size=0.1, random_state=7)\n    train_generator, val_generator = generator(xtr, xval, ytr, yval, batch_size)\n    print('imported')\n    print(xtr)\n    model = Unet(img_size)\n    print('Model summary:')\n    model.summary()\n    #model.compile(optimizer='adam', loss=bce_dice_loss, metrics=[mean_iou])\n    model.compile(optimizer='sgd', loss=bce_dice_loss, metrics=[mean_iou])\n    print('going to model fit generator now')\n    #model.fit_generator(train_generator, steps_per_epoch=len(xtr)/6, epochs=15,\n    model.fit_generator(train_generator, steps_per_epoch=25, epochs=,\n                        validation_data=val_generator, validation_steps=len(xval)/batch_size)\n    \n    preds_test = model.predict(X_test, verbose=1)\n\n    preds_test_upsampled = []\n    for i in range(len(preds_test)):\n        preds_test_upsampled.append(cv2.resize(preds_test[i], \n                                           (sizes_test[i][1], sizes_test[i][0])))\n        \n    test_ids = next(os.walk(test_path))[1]\n    new_test_ids = []\n    rles = []\n    for n, id_ in enumerate(test_ids):\n        rle = list(prob_to_rles(preds_test_upsampled[n]))\n        rles.extend(rle)\n        new_test_ids.extend([id_] * len(rle))\n        \n        \n    sub = pd.DataFrame()\n    sub['ImageId'] = new_test_ids\n    sub['EncodedPixels'] = pd.Series(rles).apply(lambda x: ' '.join(str(y) for y in x))\n    sub.to_csv('submission_tanh.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"b1ef6d8a-d810-4701-9c35-2c21fc0b7e10","_uuid":"cfcfc38799a4af19884103a4f6c4cc0d5a7c5341"},"cell_type":"markdown","source":"OK so thats about that, feel free to fork the notebook, play around with the data genarator parameters, and also the loss function, also the hyper parameters like epochs etc...\nHappy training! and if you liked the notebook please UPVOTE! :-)"}],"metadata":{"language_info":{"version":"3.6.5","nbconvert_exporter":"python","file_extension":".py","pygments_lexer":"ipython3","mimetype":"text/x-python","codemirror_mode":{"version":3,"name":"ipython"},"name":"python"},"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"}},"nbformat":4,"nbformat_minor":4}