{"cells":[{"metadata":{"_cell_guid":"0edf1de3-15a9-4f64-a560-e24cf44eefab","_uuid":"daca9bb62139ee9d5d864e3f1f26d87ba01adff5"},"cell_type":"markdown","source":"# Overview\nThis dataset contains a large number of segmented nuclei images. The images were acquired under a variety of conditions and vary in the cell type, magnification, and imaging modality (brightfield/darkfield/fluorescence). The dataset is designed to challenge an algorithm's ability to **generalize across these variations**.\n\nEach image is represented by an associated ImageId. Files belonging to an image are contained in a folder with this ImageId. Within this folder are two subfolders:\n\n- `images` contains the image file.\n- `masks` contains the segmented masks of each nucleus. This folder is only included in the training set. Each mask contains one nucleus. Masks are not allowed to overlap (no pixel belongs to two masks).\n\n**The second stage dataset will contain images from unseen experimental conditions.** To deter hand labeling, it will also contain images that are ignored in scoring. The metric used to score this competition requires that your submissions are in *run-length encoded* format. Please see the [evaluation page](https://www.kaggle.com/c/data-science-bowl-2018#evaluation) for details.\n\nAs with any human-annotated dataset, you may find various forms of errors in the data. You may manually correct errors you find in the training set. The dataset will not be updated/re-released unless it is determined that there are a large number of systematic errors. The masks of the stage 1 test set will be released with the release of the stage 2 test set.\n\nFile descriptions\n- /stage1_train/* - training set images (images and annotated masks)\n- /stage1_test/* - stage 1 test set images (images only, you are predicting the masks)\n- /stage2_test/* (released later) - stage 2 test set images (images only, you are predicting the masks)\n- stage1_sample_submission.csv - a submission file containing the ImageIds for which you must predict during stage 1\n- stage2_sample_submission.csv (released later) - a submission file containing the ImageIds for which you must predict during stage 2\n- stage1_train_labels.csv - a file showing the run-length encoded representation of the training images. This is provided as a convenience and is redundant with the mask image files."},{"metadata":{"_cell_guid":"5919cb7f-780f-4ce7-b3b9-54776aee7b84","_uuid":"f84b37e2a891d7d861655277656f4111c7b0004e"},"cell_type":"markdown","source":"**Notes:**\n\nAlpha channels\n> Simply to make things one step more confusing, some image files (notably, like PNG) contain more image information in an extra channel. This channel is essentially no different than an RGB or CMYK channel, with 256 shades of gray representing the image areas. However, Alpha channels have on crucial difference: they denote transparency, not color information.\n\n[Source](https://www.howtogeek.com/howto/42393/rgb-cmyk-alpha-what-are-image-channels-and-what-do-they-mean/)\n\nWe can simply drop this channel from our input. For example, a 256x256x4 image would become a 256x256x3 image."},{"metadata":{"collapsed":true,"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":false},"cell_type":"code","source":"import os\nimport warnings\nfrom pathlib import Path\n\nimport numpy as np\nimport pandas as pd\n\nfrom skimage.io import imread, imshow, imread_collection, concatenate_images\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\nTRAIN_DIR = Path('../input/data-science-bowl-2018/stage1_train')\nTEST_DIR = Path('../input/data-science-bowl-2018/stage1_test')\n\nIMG_TYPE = '.png'         # Image type\nIMG_CHANNELS = 3          # Default number of channels\nIMG_DIR_NAME = 'images'   # Folder name including the image\nMASK_DIR_NAME = 'masks'   # Folder name including the masks\n\nwarnings.filterwarnings('ignore', category=UserWarning, module='skimage')\nwarnings.filterwarnings('ignore', category=FutureWarning, module='skimage')","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"b57513d6-2e1e-4342-864b-9f13ad7072f4","_uuid":"96a0187f37e416b6de1e349487eaaf09f4630661"},"cell_type":"markdown","source":"# Utils\nBefore we get started, I'm going to define a few helper functions that we can use throughout this notebook. "},{"metadata":{"collapsed":true,"_cell_guid":"47afe701-6e22-404d-bc2e-a9940ffb17c5","_uuid":"2053460bb89eb27538316821bbcbbaf364a78aef","trusted":true},"cell_type":"code","source":"from scipy import ndimage as ndi\nfrom skimage.morphology import erosion, square\n\ndef read_image(observation_id, directory):\n    return imread(sorted((directory / observation_id).glob('images/*.png'))[0])\n\ndef read_masks(observation_id, directory):\n    return imread_collection(sorted((directory / observation_id).glob('masks/*.png')))\n\ndef segment_mask(masks):\n    '''Combine a list of masks into a single image.'''\n    mask = np.sum(masks, axis=0)\n    return np.clip(mask, 0, 1).astype(np.uint8)\n\ndef segment_soft_mask(masks):\n    '''\n    EXPERIMENTAL\n    Try a soft encoding for masks (as opposed to a 0/1 hard encoding) where the probability of a \n    mask is a function of the distance from the center of the nearest nucleus.\n    \n    RESULT\n    This didn't end up working out too well, the masks were way too small and required tuning of\n    the cutoff parameter. \n    '''\n    final_mask = np.zeros(masks[0].shape) # pixel locations with a value of 0 denote the background\n    for i, mask in enumerate(masks):\n        distance = ndi.distance_transform_edt(mask)\n        final_mask = np.maximum(final_mask, distance)\n    return final_mask / np.max(final_mask)\n\ndef segment_eroded_mask(masks, size=2):\n    '''Remove pixels at the boundary of a mask. Useful for ensuring that no two masks are touching.'''\n    masks = [erosion(mask, square(size)) for mask in masks]\n    mask = np.sum(masks, axis=0)\n    return np.clip(mask, 0, 1).astype(np.uint8)\n\ndef instance_mask(masks):\n    '''Returns an overlay where each instance location is labeled by an integer starting at 1 and incresasing.'''\n    all_labels = np.zeros(masks[0].shape) # pixel locations with a value of 0 denote the background\n    for i, mask in enumerate(masks):\n        mask = mask > 0\n        label = (mask)*(i+1) # pixel locations with a value of i denote the ith mask\n        all_labels = np.maximum(all_labels, label) # for overlapping masks, use the higher value - this shouldn't ever happen for this dataset\n    return all_labels.astype(np.uint8)\n\ndef separate_instances(label_image):\n    '''\n    Input: Labeled pixel map where each integer corresponds with one nucleus. \n    Returns: A list of masks where each mask shows the complete pixel mapping for one nucleus.\n    '''\n    all_masks = []\n    for i in range(1, np.max(label_image)+1):\n        mask = (label_image == i).astype(np.uint8)\n        all_masks.append(mask)\n    return all_masks","execution_count":1,"outputs":[]},{"metadata":{"_cell_guid":"99ba76af-9e2d-4ee0-8faf-0a44e2ddf248","_uuid":"c44a3859768f109d7916ab2b9db3390f5413e52a"},"cell_type":"markdown","source":"# Data exploration"},{"metadata":{"_cell_guid":"28a1b98a-e775-4f94-b8c5-1dfb472caf12","_uuid":"d4aa4bc17d69421bcc51391dfff5e6d4552ecfe4"},"cell_type":"markdown","source":"Experimental conditions were acquired from [Allen Goodman](https://www.kaggle.com/c/data-science-bowl-2018/discussion/48130#273531). We can use this information to diagnose which experimental conditions are most troubling for the model we develop. \n\n(In the notebook sidebar, click \"Add Dataset\" to use.)"},{"metadata":{"collapsed":true,"_cell_guid":"a3af7485-20fa-4216-8569-f372476cda69","_uuid":"12c1f2204d1af848b21fd58565d521ef78b4b91c","trusted":false},"cell_type":"code","source":"experimental_conditions = pd.read_csv('../input/dsb-observation-types/classes.csv')\nexperimental_conditions.head()","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"e8acddb7-1fdd-4816-becd-1693a32cabba","_uuid":"5055a3e15a0a37e8464d28aaf421d436d172ee26","trusted":false},"cell_type":"code","source":"experimental_conditions.groupby(['foreground', 'background'])['background'].agg('count')","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"0f9163eb-83db-4988-ba15-64c3f82c2034","_uuid":"c027a5f280f8af1979350d0131439e40f51cafba"},"cell_type":"markdown","source":"---"},{"metadata":{"collapsed":true,"_cell_guid":"a0996cb5-3518-4fa5-896a-1c85621cd7e1","_uuid":"d5b64fda700c25dc8bdb530ded0f3386a0049510","trusted":false},"cell_type":"code","source":"observations = os.listdir(TRAIN_DIR)\nprint(f'{len(observations)} training examples were found in {TRAIN_DIR}')\nprint(f'{len(os.listdir(TEST_DIR))} training examples were found in {TEST_DIR}')","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"0964f59b-c1bf-4fd2-a7f0-f661a7c6c922","_uuid":"b6957f573b479a12233bec4962f700d17ed0de4c","trusted":false},"cell_type":"code","source":"train_observations = observations[:-60]\nval_observations = observations[-60:]","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"39fe8966-4a45-4630-851b-6e8510cd1063","_uuid":"d00a8e3bd8460e6b9981bba6f7eb22c8647527ba"},"cell_type":"markdown","source":"Let's see what files are listed in the `images/` and `masks/` directories. We expect to see one image and a variable number of masks for each observation. "},{"metadata":{"collapsed":true,"_cell_guid":"60624fab-d5dc-4636-9992-eb5aefe473e8","_uuid":"da717200e59699644ce7b597b89584ad9507df45","trusted":false},"cell_type":"code","source":"sample = train_observations[0]\n\n# for each observation, images/ contains one photo\nimage_files = sorted((TRAIN_DIR / sample).glob('images/*.png'))\nprint(f\"Files found in 'images/': \\n{image_files}\")\n\n# insepecting an example of an image\nimage = imread(image_files[0])\nimshow(image)\n\nprint(f'\\nImage dimensions: {image.shape}')","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"6d59faf0-00dc-4ff4-aadf-fc535ed50de7","_uuid":"fe5399392d40d1ef9ca4a28696b5a254700a9a18"},"cell_type":"markdown","source":"We can combine the masks into a single image to view the labels as a semantic segmentation. "},{"metadata":{"collapsed":true,"_cell_guid":"de59efc8-e839-437d-b263-83039eccf4b0","_uuid":"ee49ab7bb11507da62f616fbc641fe741a45b4da","trusted":false},"cell_type":"code","source":"# for each observation, masks/ contains n masks where n is the number of nuclei identified in the image\nmask_files = sorted((TRAIN_DIR / sample).glob('masks/*.png'))\nprint(f\"Files found in 'masks/': \\n{mask_files}\")\n\n# insepcting an example of the masks\nmasks = imread_collection(sorted((TRAIN_DIR / observations[0]).glob('masks/*.png')))\nimshow(np.sum(masks, axis=0).astype(np.uint8)) # combine masks \n\nprint(f'\\nNumber of masks: {len(masks)}')","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"3b9b3ed3-ae3d-441e-a001-1fbfb7f70c15","_uuid":"ed2f83e9639fb9f23e3415e93a46695f7b1fd30b"},"cell_type":"markdown","source":"We're going to be building a model which outputs a segmentation map such as the one you see above. \n\nHowever, the competition would like us to output a* single mask for each nucleus* detected in the image. Using `skimage.morphology.label`, we can separate out  disconnected segments into distinct instances. This assigns an integer value for each ***object*** in the image. A pixel value of 0 indicates \"background\", meaning no nucleus was detected at that location."},{"metadata":{"collapsed":true,"_cell_guid":"06603c5c-3b97-487d-939b-8247327938d8","_uuid":"f0963de602f5eeb32665495cbba5404b5c886ece","trusted":false},"cell_type":"code","source":"from skimage.morphology import label\nfrom skimage.color import label2rgb\n\n# treat masks as a segmentation problem, then use skimage.morphology.label to identify the instances\nlabel_image = label(np.sum(masks, axis=0))\n\nimage_label_overlay = label2rgb(label_image, image=image, bg_label=0)\n\nfig, ax = plt.subplots(figsize=(10, 5))\nax.imshow(image_label_overlay)","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"499a8bff-f82a-4ce3-9dad-1538d361c06f","_uuid":"7fe61b79ec33b0bb75f9520d1f774100aa05dcd4","trusted":false},"cell_type":"code","source":"masks = separate_instances(label_image)\n\n# Show the first 8 masks in a labeled image\nfig, ax = plt.subplots(ncols=4, nrows=2, figsize=(14, 8))\nfig.suptitle(\"Masks for image\", fontsize=16)\nfor i in range(min(len(masks), 8)):\n    ax[i // 4, i % 4].imshow(masks[i])\nfor x in ax.ravel(): x.axis(\"off\")","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"fa67801e-14c5-42af-aca3-1d79ced4a597","_uuid":"2311f3ecd588c473408782a176b0fd9ce9dbdd32","trusted":false},"cell_type":"code","source":"from skimage.measure import regionprops\nimport matplotlib.patches as mpatches\nimport random\n\ndef show_random_observations(observations, n=3, bboxes=False):\n    fig, ax = plt.subplots(ncols=2, nrows=n, figsize=(20, n*6))\n\n    for i in range(n):\n        # Load an example image\n        sample = random.choice(observations)\n        image = read_image(sample, TRAIN_DIR)\n        masks = instance_mask(read_masks(sample, TRAIN_DIR))\n        image_label_overlay = label2rgb(masks, image=image, bg_label=0)\n\n        ax[i, 0].imshow(image)\n        ax[i, 0].set_title('input image')\n\n        ax[i, 1].imshow(image_label_overlay)\n        ax[i, 1].set_title('target labels')\n\n        # also show bounding boxes just for fun\n        if bboxes:\n            for region in regionprops(masks):\n                # take regions with large enough areas\n                if region.area >= 1:\n                    # draw rectangle around segmented coins\n                    minr, minc, maxr, maxc = region.bbox\n                    rect = mpatches.Rectangle((minc, minr), maxc - minc, maxr - minr,\n                                              fill=False, edgecolor='red', linewidth=2)\n                    ax[i, 1].add_patch(rect)\n\n        print(sample) # in case I want to make a note of a specific observation","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"143c6f89-5e54-4d67-a88a-18c1a1343668","_uuid":"8c499662e801989e7118d0b9096184651cbcb89a","trusted":false},"cell_type":"code","source":"show_random_observations(train_observations, n=4)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"e4ce7159-7ed5-4068-b545-0daa64c4b20b","_uuid":"6855be543909d8365d7515aeabe08e292ff24621"},"cell_type":"markdown","source":"Let's get some basic stats about our training data. "},{"metadata":{"collapsed":true,"_cell_guid":"25c07a20-dadd-4755-a667-530a35af85ad","_uuid":"77c23f5e3a14657ec07855be23c7ba68c0bac44e","trusted":false},"cell_type":"code","source":"from collections import namedtuple\n\n# let's get some basic stats about our dataset \nSummary = namedtuple('Summary', ['observation_id', 'image_size', 'n_masks'])\ntrain_data_summary = []\n\n# create data frame of observation_id, image_size, n_masks\nfor observation_id in train_observations:\n    train_data_summary.append(Summary(observation_id=observation_id,\n                                      image_size=read_image(observation_id, TRAIN_DIR).shape,\n                                      n_masks=len(read_masks(observation_id, TRAIN_DIR))\n                                     )\n                             )\n    \ndf = pd.DataFrame(train_data_summary, columns=Summary._fields)\ndf.head()","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"6bfd6a23-cc64-4b78-9a3f-c3de27ab1702","_uuid":"81d03d8fe56b71067748b6e6bd8b5e7824efaf01","trusted":false},"cell_type":"code","source":"df['image_size'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"280c5b8d-efe1-4ff2-b141-87af2789014b","_uuid":"b1e14033b32890b6008bcfe97cb32f1eecf73dd7","trusted":false},"cell_type":"code","source":"df['n_masks'].plot(kind='hist')","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"a7bcf926-e0f2-43bb-945d-eb6f281a83f0","_uuid":"a3545b056cd0b0055e2c85a33db5aeb5e2fa9e11"},"cell_type":"markdown","source":"The majority of images have between 0 and 50 nuceli, although some have up to 350 nuclei in a single image!"},{"metadata":{"_cell_guid":"3466c281-291d-4619-a2fc-cd05c89f24f3","_uuid":"edb423e42838e7d54d2d98d5907280e37f06dd96"},"cell_type":"markdown","source":"# Strategy\nThere's two approaches we could take here. \n\n1. We could treat this as a semantic segmentation problem first, and then attempt to separate out the semantic mask into individual nuclei masks. \n2. We could directly treat this as an instance segmentation problem, both detecting the nuclei in the image and calculating the corresponding mask for each nucleus in an end to end fashion. \n\nThe [U-Net](https://lmb.informatik.uni-freiburg.de/people/ronneber/u-net/) architecture is popular for image segmentation tasks. This outputs a pixel map of all of the regions of the image where a nucleus is predicted, but it does not distinguish between nuclei.  \n\n![](https://i.imgur.com/GKBDAzD.png)\n\nIt's been reported to perform [quite well](https://www.kaggle.com/c/data-science-bowl-2018/discussion/54426) for this competition. **We'll build a U-Net model in this notebook. **\n\n[Mask R-CNN](https://research.fb.com/publications/mask-r-cnn/) is a popular instance segmentation model. You can see an implementation of this model discussed [here](https://www.kaggle.com/c/data-science-bowl-2018/discussion/54089).\n\n![](https://cdn-images-1.medium.com/max/800/1*IWWOPIYLqqF9i_gXPmBk3g.png)\n\nIf we wanted to compare our deep learning method against a baseline, we could consider the [watershed approach](http://scikit-image.org/docs/dev/auto_examples/segmentation/plot_watershed.html#sphx-glr-auto-examples-segmentation-plot-watershed-py) as a traditional CV method to compare our model against."},{"metadata":{"_cell_guid":"c0ab36c9-c6c4-47ac-8e41-a360b07c5897","_uuid":"e00eb0644ef0bc3d3244199e23de3ba477fa17dc"},"cell_type":"markdown","source":"# Preprocessing techniques\n\nIn the [U-Net paper](https://arxiv.org/abs/1505.04597), the authors discuss a loss weighting scheme which encourages their model to pay close attention to the segmentation map at the boundary of cells. \n\n![](https://i.imgur.com/CUiEfOU.png)\n\nKeras doesn't currently have native support for weighting a pixel-wise loss function in this manner. However, we can cleverly sneak this information into our target as suggested [here](https://github.com/keras-team/keras/issues/6629#issuecomment-302756549). \n\n![](https://i.imgur.com/gu4UG1i.png)\n\nThis may seem a bit awkward, but we've successfully developed a manner to encode our loss weighting scheme into the target such that when it comes time to calculate the loss, we can disentangle the weights from the mask, calculate the loss using the mask, and then weight the loss accordingly. \n\nIf you're building a model in Tensorflow or PyTorch you can weight your loss function accordingly without needing this \"hack\"."},{"metadata":{"collapsed":true,"_cell_guid":"139931c9-8ec5-4bac-9e10-fb13df076c0a","_uuid":"daff9bef9dca776b2e5e285c80e9fdd139872ac3","trusted":false},"cell_type":"code","source":"def encode_target(masks, w0=5, sigma=2):\n    # ref : https://www.kaggle.com/piotrczapla/tensorflow-u-net-starter-lb-0-34/notebook\n    masks = [erosion(mask, square(2)) for mask in masks]\n    \n    merged_mask = segment_mask(masks)\n    weight = np.zeros(merged_mask.shape)\n    # calculate weight for important pixels\n    distances = np.array([ndi.distance_transform_edt(m==0) for m in masks])\n    shortest_dist = np.sort(distances, axis=0)\n    # distance to the border of the nearest cell \n    d1 = shortest_dist[0]\n    # distance to the border of the second nearest cell\n    d2 = shortest_dist[1] if len(shortest_dist) > 1 else np.zeros(d1.shape)\n\n    weight = w0 * np.exp(-(d1+d2)**2/(2*sigma**2)).astype(np.float32)\n    weight = 1 + (merged_mask == 0) * weight\n    return merged_mask - weight\n\ndef decode_target(encoding):\n    target_mask = np.array(encoding == 0, dtype=np.uint8)\n    weights = (-1 * encoding) + target_mask\n    \n    return target_mask, weights","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"cbd5cd7b-6d85-4519-b4ab-04211e344103","_uuid":"1f8c106a74bfadaf01cdf8385fd0826e04d4ced4","trusted":false},"cell_type":"code","source":"encoding = encode_target(read_masks(train_observations[12], TRAIN_DIR))\nmask, weight = decode_target(encoding)\n\nfig, ax = plt.subplots(ncols=3, figsize=(16,8))\n\nax[0].imshow(mask) \nax[0].set_title('Semantic segmentation')\n\nax[1].imshow(weight) \nax[1].set_title('Separation weights')\n\nax[2].imshow(mask + weight) \nax[2].set_title('Weights imposed on segmentation mask')","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"0fa27073-9e6d-4e96-9732-0f22bc0b6535","_uuid":"c772a50e4e6bb3479a8635e940136d722d59a5b9","trusted":false},"cell_type":"code","source":"from keras.losses import binary_crossentropy\nimport keras.backend as K\n\ndef weighted_binary_crossentropy(y_true, y_pred):\n    '''\n    Calculates the weighted pixel-wise binary cross entropy. Expects target to be encoded as `(mask - weights)`. \n    '''\n    # mask <- where value==0\n    target_mask = K.cast(K.equal(y_true, 0), 'float32') \n    \n    # weights calculated as described above\n    weights = (-1 * y_true) + target_mask\n    \n    cce = binary_crossentropy(target_mask, y_pred)  \n    wcce = cce * K.squeeze(weights, axis=-1)\n    return K.mean(wcce, axis=-1)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"c5680a8c-a656-4753-9d0b-e542d5462440","_uuid":"13318d21c8d41f2d2a23c700bb9ab4903e7e0c04"},"cell_type":"markdown","source":"# Data preparation"},{"metadata":{"collapsed":true,"_cell_guid":"3a301eb0-f50f-4203-aa99-3d1e0d44b240","_uuid":"abd2a997fdd6da3e9f6e9a1f0e2cad4a27a32d39","trusted":false},"cell_type":"code","source":"def prepare_target(observation):\n    masks = read_masks(observation, TRAIN_DIR)\n    masks_resized = [cv2.resize(mask, (IMG_HEIGHT, IMG_WIDTH), interpolation=cv2.INTER_AREA) for mask in masks]\n    encoding = encode_target(masks_resized)\n    return encoding","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"5e5c27a3-7830-4ca5-8768-dc96141bcefc","_uuid":"a77ad9955d811fdd5542308f74fc695df67172fd","trusted":false},"cell_type":"code","source":"import cv2\n\nIMG_WIDTH = 256\nIMG_HEIGHT = 256\nIMG_CHANNELS = 3\n\n# images\nx_train = [read_image(observation, TRAIN_DIR) for observation in train_observations]\nx_train = [cv2.resize(image, (IMG_HEIGHT, IMG_WIDTH), interpolation=cv2.INTER_AREA) for image in x_train]\nx_train = np.array(x_train, dtype=np.uint8)\nx_train = x_train[:,:,:,:IMG_CHANNELS]\n\nx_val = [read_image(observation, TRAIN_DIR) for observation in val_observations]\nx_val = [cv2.resize(image, (IMG_HEIGHT, IMG_WIDTH), interpolation=cv2.INTER_AREA) for image in x_val]\nx_val = np.array(x_val, dtype=np.uint8)\nx_val = x_val[:,:,:,:IMG_CHANNELS]\n\n# targets\ny_train = [prepare_target(observation) for observation in train_observations]\ny_train = np.array(y_train)\ny_train = np.expand_dims(y_train, axis=-1)\n\ny_val = [prepare_target(observation) for observation in val_observations]\ny_val = np.array(y_val)\ny_val = np.expand_dims(y_val, axis=-1)","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"bdfad6e9-a34b-4844-9729-bd0e29b2e508","_uuid":"aa4c1bf4e4b18c2abd48b94b2275ce20e4c74c46","trusted":false},"cell_type":"code","source":"ix = 12\n\nencoding = y_train[ix]\nmask, weights = decode_target(encoding)\nfig, ax = plt.subplots(ncols=3, figsize=(16,8))\nax[0].imshow(x_train[ix])\nax[0].set_title('Image')\nax[1].imshow(np.squeeze(mask))\nax[1].set_title('Mask')\nax[2].imshow(np.squeeze(weights))\nax[2].set_title('Weights')","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"eea84808-a5ea-4085-b8dd-31706a917597","_uuid":"a445789cdbea2a131ea104b4c19b3ecd55a1f24a"},"cell_type":"markdown","source":"We'll create ***generator*** objects for randomly flipping our inputs in order to augment our dataset size. "},{"metadata":{"collapsed":true,"_cell_guid":"0de45e41-8050-4f11-8b61-0da56491ed5c","_uuid":"35162d8ff3530efb4722dd906d5ce78de806a40a","trusted":false},"cell_type":"code","source":"from keras.preprocessing.image import ImageDataGenerator\n\ndata_gen_args = dict(horizontal_flip=True, \n                     vertical_flip=True)\n\nimage_datagen = ImageDataGenerator(rescale=1./255, **data_gen_args)\nmask_datagen = ImageDataGenerator(**data_gen_args)\n\n\n# Provide the same seed and keyword arguments to the flow methods\nseed = 1\nbatch_size = 4\n\n# ------ training data ------\ntrain_image_generator = image_datagen.flow(x_train, batch_size=batch_size, seed=seed)\ntrain_mask_generator = mask_datagen.flow(y_train, batch_size=batch_size, seed=seed)\n\n# combine generators into one which yields image and masks\ntrain_generator = zip(train_image_generator, train_mask_generator)\ntrain_steps = np.ceil(len(x_train) / batch_size)\n\n# ------ validation data ------\nval_image_generator = image_datagen.flow(x_val, batch_size=batch_size, seed=seed)\nval_mask_generator = mask_datagen.flow(y_val, batch_size=batch_size, seed=seed)\n\n# combine generators into one which yields image and masks\nval_generator = zip(val_image_generator, val_mask_generator)\nval_steps = np.ceil(len(x_val) / batch_size)","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"a9f5c648-859f-49f0-a91f-80fd332194db","_uuid":"e7c89da325957052910ce20c34aca02fbeba0f75","trusted":false},"cell_type":"code","source":"image, encoding = next(train_generator)\nmask, weights = decode_target(encoding)\n\nfig, ax = plt.subplots(ncols=3, figsize=(16, 8))\nax[0].imshow(image[0])\nax[1].imshow(np.squeeze(mask[0]))\nax[2].imshow(np.squeeze(weights[0]))","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"666a609d-51b1-4557-8172-cfff10762f0c","_uuid":"7c0351a7835d695c8d6020822643b72245ce05f3"},"cell_type":"markdown","source":"# Standard U-Net Model"},{"metadata":{"collapsed":true,"_cell_guid":"69615cee-2099-4334-95ee-0b17a610b556","_uuid":"ade25d3006eb61949e0dc83aa8ab6040e765bb78","trusted":false},"cell_type":"code","source":"# define the u-net\nfrom keras.models import Model, load_model\nfrom keras.layers import Input, BatchNormalization, Activation\nfrom keras.layers.core import Lambda, SpatialDropout2D\nfrom keras.layers.convolutional import Conv2D, Conv2DTranspose\nfrom keras.layers.pooling import MaxPooling2D\nfrom keras.layers.merge import concatenate, add\nfrom keras.callbacks import EarlyStopping, ModelCheckpoint\nfrom keras import regularizers\nfrom keras import backend as K\n\nimport tensorflow as tf","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"aa9d21fe-44be-4a5b-8081-62d8d06f750b","_uuid":"ba2e8c2930b7a0337c520e21420939f05cb56d2e","trusted":false},"cell_type":"code","source":"def conv_block(inputs, filters, filter_size=3, drop_prob=0.2, regularizer=regularizers.l2(0.0001)):\n    x = Conv2D(filters, filter_size, padding='same', kernel_regularizer=regularizer)(inputs)\n    x = BatchNormalization()(x)\n    x = Activation('relu')(x)\n    x = Conv2D(filters, filter_size, padding='same', kernel_regularizer=regularizer)(x)\n    x = BatchNormalization()(x)\n    x = Activation('relu')(x)\n    x = SpatialDropout2D(drop_prob)(x)\n    return x\n\ndef downsample(block):\n    x = MaxPooling2D(pool_size=(2, 2)) (block)\n    return x\n\ndef upsample(block, skip_connection, filters, regularizer=regularizers.l2(0.0001)):\n    x = Conv2DTranspose(filters, (3, 3), strides=(2, 2), padding='same', kernel_regularizer=regularizer)(block)\n    stack = concatenate([skip_connection, x])\n    return stack","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"7f3dc05c-cebd-49ab-831c-6de3411f6f3b","_uuid":"ce0c7f2d73e19fba941e401cfc12403ec351f874","trusted":false},"cell_type":"code","source":"from keras.optimizers import SGD\n\ndef build_unet(IMG_WIDTH=256, IMG_HEIGHT=256, IMG_CHANNELS=3, drop_prob=0.2):\n    \n    regularizer=regularizers.l2(0.0001)\n    \n    # ---- Model ----\n    inputs = Input((IMG_HEIGHT, IMG_WIDTH, IMG_CHANNELS))\n    \n    # Downsample\n    encode_1 = conv_block(inputs, 16, regularizer=regularizer)\n    down_1 = downsample(encode_1)\n    \n    encode_2 = conv_block(down_1, 32, regularizer=regularizer)\n    down_2 = downsample(encode_2)\n    \n    encode_3 = conv_block(down_2, 64, regularizer=regularizer)\n    down_3 = downsample(encode_3)\n    \n    encode_4 = conv_block(down_3, 128, regularizer=regularizer)\n    down_4 = downsample(encode_4)\n    \n    bridge = conv_block(down_4, 256, regularizer=regularizer)\n    \n    up_4 = upsample(bridge, encode_4, 128)\n    decode_4 = conv_block(up_4, 128, regularizer=regularizer)\n\n    up_3 = upsample(decode_4, encode_3, 64)\n    decode_3 = conv_block(up_3, 64, regularizer=regularizer)\n\n    up_2 = upsample(decode_3, encode_2, 32)\n    decode_2 = conv_block(up_2, 32, regularizer=regularizer)\n\n    up_1 = upsample(decode_2, encode_1, 16)\n    decode_1 = conv_block(up_1, 16, regularizer=regularizer)\n\n    outputs = Conv2D(1, (1, 1), activation='sigmoid', kernel_regularizer=regularizer)(decode_1)\n\n    model = Model(inputs=[inputs], outputs=[outputs])\n    model.compile(optimizer='adam', loss=weighted_binary_crossentropy)\n    return model\n\nmodel = build_unet()","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"c0a40772-df35-4cd4-8042-f6fac5ece4d9","_uuid":"d7e91d54a8c7495f9073a1d31d6b96afd00c7586","scrolled":false,"trusted":false},"cell_type":"code","source":"model.summary()","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"f9809672-3c7d-4a58-bbd2-87d53ef33b0b","_uuid":"8d5318e2cdca8273d4b407a932622c398fa9b0ae","trusted":false},"cell_type":"code","source":"from keras.callbacks import LearningRateScheduler\n\ndef step_decay_schedule(initial_lr=1e-3, decay_factor=0.75, step_size=10):\n    '''\n    Wrapper function to create a LearningRateScheduler with step decay schedule.\n    '''\n    def schedule(epoch):\n        return initial_lr * (decay_factor ** np.floor(epoch/step_size))\n    \n    return LearningRateScheduler(schedule)\n\nlr_sched = step_decay_schedule(initial_lr=1e-3, decay_factor=0.90, step_size=10)","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"325dc36b-3e0f-464a-b2fb-4a37d6665229","_uuid":"2e81a00b6e2a72348324bc0ee3b2dba105ba2445","trusted":false},"cell_type":"code","source":"checkpointer = ModelCheckpoint('unet_best.h5', verbose=1, save_best_only=True)\n\nresults = model.fit_generator(train_generator, steps_per_epoch=train_steps, epochs=75, \n                              validation_data=val_generator, validation_steps=val_steps,\n                              callbacks=[checkpointer, lr_sched])","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"a71fbfbb-bda2-4f49-ba56-6323c54f74c9","_uuid":"94e8e1143051ac920ea72e578bc218e08f8d9725","trusted":false},"cell_type":"code","source":"model.load_weights('unet_best.h5')","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"4552c857-d873-450a-9e2f-6525b223e211","_uuid":"53e5783938b19494db82d07e08ab62ae1d24861a","trusted":false},"cell_type":"code","source":"plt.plot(results.history['loss'])\nplt.plot(results.history['val_loss'])\nplt.title('model loss')\nplt.ylabel('loss')\nplt.xlabel('epoch')\nplt.legend(['train', 'val'], loc='upper left')","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"fb3c57d5-e95d-4873-b879-73f17a59dcdd","_uuid":"394e42b96cdb6411096e285c2222c7c21eb92005"},"cell_type":"markdown","source":"### Quickly inspect model performance"},{"metadata":{"collapsed":true,"_cell_guid":"7baa5cac-9a23-4436-be9d-ae9ed81fa647","_uuid":"f1cccc28cf7a256f5c40da90096b00663ba91c2b","trusted":false},"cell_type":"code","source":"def show_results(ix, model):\n    '''Quick helper function to display predictions.''' \n    image = x_val[ix]\n    key = val_observations[ix]\n\n    masks = read_masks(key, TRAIN_DIR)\n    masks_resized = [cv2.resize(mask, (IMG_HEIGHT, IMG_WIDTH), interpolation=cv2.INTER_AREA) for mask in masks]\n    target = label2rgb(instance_mask(masks_resized), image=image, bg_label=0)\n\n\n    pred_mask = model.predict((x_val[ix]/255)[None])[0]\n    pred_mask = np.squeeze((pred_mask > 0.5).astype(np.uint8))\n    label_image = label(pred_mask)\n    image_label_overlay = label2rgb(np.squeeze(label_image), image=image, bg_label=0)\n\n    mask, weights = decode_target(np.squeeze(y_val[ix]))\n\n    fig, ax = plt.subplots(2, 3, figsize=(20, 12))\n\n    ax[0,0].imshow(image)\n    ax[0,0].set_title('input')\n\n    ax[1,0].imshow(weights)\n    ax[1,0].set_title('loss weights')\n\n    ax[0,1].imshow(mask)\n    ax[0,1].set_title('semantic target')\n\n    ax[0,2].imshow(target)\n    ax[0,2].set_title('instance target')\n\n    ax[1,1].imshow(pred_mask*1.0)\n    ax[1,1].set_title('prediction: segmented mask')\n\n    ax[1,2].imshow(image_label_overlay)\n    ax[1,2].set_title('prediction: instance mask')\n\n    for x in ax.ravel(): x.axis(\"off\")","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"b177752d-e719-4dc2-b6af-eadbe4e3bd58","_uuid":"96ef40cf3cbf8a62f64baa26ca58f0f285cc83d0","trusted":false},"cell_type":"code","source":"show_results(48, model)","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"c3755537-2456-4c24-89b3-8710264686e2","_uuid":"5ab0d6a4b4e22ad4b7ba06b794987aeaf49ab934","trusted":false},"cell_type":"code","source":"show_results(54, model)","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"6233ea15-5c94-43cc-b6fd-3d2025bcf466","_uuid":"797291cd852b71a706d7ee1970244e8512706ce4","trusted":false},"cell_type":"code","source":"show_results(26, model)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"7d176a3d-088a-4952-bc53-0fa824abba94","_uuid":"7850602219b780cf10fba1ec3d3562b55a956ad8"},"cell_type":"markdown","source":"The simple U Net works well for easy images, but struggles on some of the more challenging ones. This was also observed by the [competition winners](https://www.kaggle.com/c/data-science-bowl-2018/discussion/54741) and their solution was to go deeper. Unfortunately, we're limited in our capacity to build deeper models due to Kaggle's memory constraints on kernels. "},{"metadata":{"collapsed":true,"_cell_guid":"e7087a27-6a69-4243-a4e1-ebd0720b9230","_uuid":"1249b0466d5c43ed98315d2c8fe2d01b4a01b66f","trusted":false},"cell_type":"code","source":"# Kaggle Kernels have limited memory, so we'll remove this model from memory\nimport gc\ndel model \ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"7951bb7f-ddff-448d-9403-84e9b55cc82b","_uuid":"e33837a4bd91b986ba14e3018d29cf049a658c9f"},"cell_type":"markdown","source":"# Fully Convolutional DenseNet model\n\nNext, we'll explore a  more advanced architecture composed of dense convolution blocks. These dense blocks allow for any layer to access all of the feature maps from previous layers within a block; this reuse allows us to perform convolutions at much lower channel depths and simply copy the earlier feature maps.\n\n[The One Hundred Layers Tiramisu: Fully Convolutional DenseNets for Semantic Segmentation](https://arxiv.org/abs/1611.09326)\n\n![](https://i.imgur.com/g0p6LPY.png)"},{"metadata":{"collapsed":true,"_cell_guid":"74090ff8-94b5-4750-8073-ce54607e9da2","_uuid":"a6f8ecd29a793eb8794020749d8850c5d045c1ec","trusted":false},"cell_type":"code","source":"def dense_block(stack, n_layers, growth_rate, filter_size=3, drop_prob=0.2):\n    \n    for layer in range(n_layers):\n        x = BatchNormalization()(stack)\n        x = Activation('relu')(x)\n        x = Conv2D(growth_rate, filter_size, padding='same', kernel_regularizer=regularizers.l2(0.0001)) (x)\n        x = SpatialDropout2D(drop_prob)(x)\n        stack = concatenate([stack, x])\n        \n    return stack\n\ndef downsample(block, n_filters, drop_prob=0.2):\n    x = BatchNormalization()(block)\n    x = Activation('relu')(x)\n    x = Conv2D(n_filters, (1, 1), padding='same', kernel_regularizer=regularizers.l2(0.0001)) (x)\n    x = SpatialDropout2D(drop_prob)(x)\n    x = MaxPooling2D(pool_size=(2, 2)) (x)\n    return x\n\ndef upsample(block, skip_connection, n_filters):\n    x = Conv2DTranspose(n_filters, (3, 3), strides=(2, 2), padding='same', kernel_regularizer=regularizers.l2(0.0001)) (block)\n    stack = concatenate([skip_connection, x])\n    return stack","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"e21fbce2-b9df-4d64-8fe3-5fffb59df279","_uuid":"eecf47a65ddf6a7264861aa7f2a22ccc6fe7435e","trusted":false},"cell_type":"code","source":"def build_fcdensenet(m, IMG_WIDTH=256, IMG_HEIGHT=256, IMG_CHANNELS=3):\n    '''\n    Keras Implementation of the Fully Convolutional DenseNet. \n    m: List containing the number of feature maps for each dense block. \n    '''\n    regularizer = regularizers.l2(0.0001)\n    \n    # ---- Model ----\n    inputs = Input((IMG_HEIGHT, IMG_WIDTH, IMG_CHANNELS))\n    conv = Conv2D(m[0], (3,3), padding='same', kernel_regularizer=regularizer) (inputs)\n\n    encode_1 = dense_block(conv, 4, 12)\n    down_1 = downsample(encode_1, m[1])\n\n    encode_2 = dense_block(down_1, 4, 12)\n    down_2 = downsample(encode_2, m[2])\n\n    encode_3 = dense_block(down_2, 4, 12)\n    down_3 = downsample(encode_3, m[3])\n\n    encode_4 = dense_block(down_3, 4, 12)\n    down_4 = downsample(encode_4, m[4])\n\n    bridge = dense_block(down_4, 4, 12)\n\n    up_4 = upsample(bridge, encode_4, m[6])\n    decode_4 = dense_block(up_4, 4, 12)\n\n    up_3 = upsample(decode_4, encode_3, m[7])\n    decode_3 = dense_block(up_3, 4, 12)\n\n    up_2 = upsample(decode_3, encode_2, m[8])\n    decode_2 = dense_block(up_2, 4, 12)\n\n    up_1 = upsample(decode_2, encode_1, m[9])\n    decode_1 = dense_block(up_1, 4, 12)\n\n    outputs = Conv2D(1, (1, 1), activation='sigmoid', kernel_regularizer=regularizer) (decode_1)\n\n    model = Model(inputs=[inputs], outputs=[outputs])\n    model.compile(optimizer='adam', loss=weighted_binary_crossentropy)\n    return model","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"01749161-42b0-4a0a-9b69-05a006623079","_uuid":"d6f1e1d9222d2ad24bfde120b7d4f7c723ef065f","trusted":false},"cell_type":"code","source":"# determine number of feature maps for each layer, following guidance from the original paper\nm = [32]\nencoder_blocks = 5 # 4 dense + bridge\ndecoder_blocks = 4\nfor i in range(encoder_blocks):\n    m.append(m[i] + 4*12)\n    \nfor i in range(decoder_blocks):\n    m.append(m[4-i] + 4*12 + 4*12) # skip connection + feature maps from the upsampled block + feature maps in the new block\n\ndense_model = build_fcdensenet(m)","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"484cec99-2b4d-4608-a37d-ddafd2505434","_uuid":"f1328bd625b603ca201e50f4175f9845881a74f5","trusted":false},"cell_type":"code","source":"dense_model.summary()","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"e8623355-ed21-48b6-b982-965ec3b78a7d","_uuid":"dbb5487b27fda286a10c9b8e2dad45328300581a","trusted":false},"cell_type":"code","source":"lr_sched = step_decay_schedule(initial_lr=1e-3, decay_factor=0.90, step_size=10)\ncheckpointer = ModelCheckpoint('fcdense_best.h5', verbose=1, save_best_only=True)\n\nresults = dense_model.fit_generator(train_generator, steps_per_epoch=train_steps, epochs=75, \n                                   validation_data=val_generator, validation_steps=val_steps,\n                                   callbacks=[checkpointer, lr_sched])","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"2a12e6c5-8ba7-4caa-8061-c5779bffe686","_uuid":"8ed063f0d24e96edafebbdff9663354344d0b5f2","trusted":false},"cell_type":"code","source":"dense_model.load_weights('fcdense_best.h5')","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"29ff2017-b5d5-4f29-b009-b219b8e4581e","_uuid":"2d1bdfd304a57345794051d7e8eeeef30e6c142e","trusted":false},"cell_type":"code","source":"plt.plot(results.history['loss'])\nplt.plot(results.history['val_loss'])\nplt.title('model loss')\nplt.ylabel('loss')\nplt.xlabel('epoch')\nplt.legend(['train', 'val'], loc='upper left')","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"71e30f4d-7157-4ade-ac40-495e2a04bee8","_uuid":"c4a09e5294696f0d5c8ffc7930c1e2ccae1f841c","trusted":false},"cell_type":"code","source":"show_results(48, dense_model)","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"1ae64a2d-951e-473b-8fd4-83ae97d6c1e2","_uuid":"04247b99cc1c3ee8ab1c1f5a42e59e82cb5e92ce","trusted":false},"cell_type":"code","source":"show_results(54, dense_model)","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"13b9fda5-eb40-4f8f-ade4-445934f5d0b2","_uuid":"141443c1e9cba3677fd204c85e974a9f802451ee","trusted":false},"cell_type":"code","source":"show_results(26, dense_model)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"d8518375-4f13-4b18-bfc0-aa582aadd973","_uuid":"be9679f05c82f3ef0ed164a7f04e8f603c6a057e"},"cell_type":"markdown","source":"# Evaluation\n\nThis competition evaluates submissions according to the mean average precision. \n\nI'll discuss the relevant concepts for calculating this value below. \n\n\n### Intersection over union\nThe IoU metric is a measure for how close our prediction is to the true label. We'll use this value to determine whether or not our mask prediction was \"successful\", as defined by having an IoU score above some specified threshold. \n![](https://i.imgur.com/cgd1ZNX.png)\n\n### Calculating precision \n![](https://i.imgur.com/09FRSWM.png)\n\nThe precision can be calculated as: \n$$precision = \\frac{TP}{(TP + FP + FN)}$$"},{"metadata":{"_cell_guid":"f8dfce03-8efe-495a-b110-d3a06205c532","_uuid":"20602e3fcd93f15be286df8ff654f3e5cb4b54f4"},"cell_type":"markdown","source":"We'll create a matrix where we measure the IoU of each predicted mask against each target mask. \n\n![](https://i.imgur.com/Oswiisg.png)\n\n### Average precision\n\nFor this competition, the precision is calculated by averaging over a range of IoU thresholds. \n\n![](https://i.imgur.com/H4D5nEI.png)\n"},{"metadata":{"collapsed":true,"_cell_guid":"7dfff8e4-14f5-446a-b3ac-d38d67e59221","_uuid":"442b715dcf934542ee4a5b078297b4dc69e2641f","trusted":false},"cell_type":"code","source":"def iou_at_thresholds(target_mask, pred_mask, thresholds=np.arange(0.5,1,0.05)):\n    '''Returns True if IoU is greater than the thresholds.'''\n    intersection = np.logical_and(target_mask, pred_mask)\n    union = np.logical_or(target_mask, pred_mask)\n    iou = np.sum(intersection > 0) / np.sum(union > 0)\n    return iou > thresholds\n\ndef calculate_iou_tensor(target_masks, pred_masks, thresholds=np.arange(0.5,1,0.05)):\n    iou_tensor = np.zeros([len(thresholds), len(pred_masks), len(target_masks)])\n\n    # TODO: Use tiling to make this faster\n    for i, p_mask in enumerate(pred_masks):\n        for j, t_mask in enumerate(target_masks):\n            iou_tensor[:, i, j] = iou_at_thresholds(t_mask, p_mask, thresholds)\n\n    return iou_tensor\n\ndef calculate_average_precision(target_masks, pred_masks, thresholds=np.arange(0.5,1,0.05)):\n    '''Calculates the average precision over a range of thresholds for one observation (with a single class).'''\n    iou_tensor = calculate_iou_tensor(target_masks, pred_masks, thresholds=thresholds)\n    \n    TP = np.sum((np.sum(iou_tensor, axis=2) == 1), axis=1)\n    FP = np.sum((np.sum(iou_tensor, axis=1) == 0), axis=1)\n    FN = np.sum((np.sum(iou_tensor, axis=2) == 0), axis=1)\n\n    precision = TP / (TP + FP + FN)\n\n    return np.mean(precision)\n\ndef calculate_mean_average_precision(y_true, y_pred):\n    '''\n    # Arguments\n        y_true: A list of lists each containing the target masks for a given observation.\n        y_pred: A list of lists each containing the predicted masks for a given observation.\n    '''\n    average_precision = []\n    thresholds=np.arange(0.5,1,0.05)\n    \n    for target, prediction in zip(y_true, y_pred):\n        ap = calculate_average_precision(target, prediction, thresholds=thresholds)\n        average_precision.append(ap)\n        \n    return average_precision","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"653d500d-24b6-4e91-9732-3dcb8087216e","_uuid":"faef342c1124e659224517894a70e3e483e90ed1","trusted":false},"cell_type":"code","source":"# resize individual masks for final evaluation\ndef resize_masks(masks):\n    return [cv2.resize(mask, (IMG_HEIGHT, IMG_WIDTH), interpolation=cv2.INTER_AREA) for mask in masks]\n\ny_val_masks = [read_masks(observation, TRAIN_DIR) for observation in val_observations]\ny_val_masks = [resize_masks(masks) for masks in y_val_masks]","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"573d22eb-f8cc-4399-9da8-7596f0b57e4b","_uuid":"de4e18a4e6eed86b061f53c0e44db60ad9d2526d","trusted":false},"cell_type":"code","source":"from skimage.morphology import dilation\n\ndef separate_instances_with_tricks(label_image):\n    pred_masks = []\n    for i in range(1, np.max(label_image)+1):\n        mask = (label_image == i).astype(np.uint8)\n        if np.sum(mask) > 5:\n            dilated_mask = dilation(mask, square(3))\n            pred_masks.append(mask)\n    return pred_masks\n\npreds = dense_model.predict(x_val/255)\npreds = np.squeeze((preds > 0.5).astype(np.uint8))\npreds_masks = [label(pred) for pred in preds]\npreds_masks = [separate_instances_with_tricks(label) for label in preds_masks]","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"48bb56b8-84de-407a-b789-80c7e4d1a925","_uuid":"e2018444d7073932a45753b59bfed5e8a0b7421f","trusted":false},"cell_type":"code","source":"avg_prec = calculate_mean_average_precision(y_val_masks, preds_masks)","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"4243564f-b718-4a98-b76b-9be4156a67c4","_uuid":"33fe2bd10cfacd512899f005c15f1bf30d640730","trusted":false},"cell_type":"code","source":"import seaborn as sns\nsns.distplot(avg_prec)\nprint(f'Mean average precision: {np.mean(avg_prec)}')","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"78c5157f-37d0-4993-94d7-8c313f63cbfe","_uuid":"258f2729b80e15b02265519a01b293176c977fac"},"cell_type":"markdown","source":"# Generate test predictions"},{"metadata":{"collapsed":true,"_cell_guid":"193db631-693a-4981-816b-dc61faeea8c2","_uuid":"5d671e6656643d128f85bf896cacffc0075a1fac","trusted":false},"cell_type":"code","source":"test_data = {}\n\n# read in the data\nfor observation in os.listdir(TEST_DIR):\n    test_data[observation] = {'image': read_image(observation, TEST_DIR),\n                              'size': read_image(observation, TEST_DIR).shape[:2]}","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"1c931c05-9553-435a-b5b5-cfb67e84afe8","_uuid":"cff329afdad4a339cd10a84726528e2fd205d94e","trusted":false},"cell_type":"code","source":"x_test = [observation['image'] for observation in test_data.values()]\nx_test = [cv2.resize(image, (IMG_HEIGHT, IMG_WIDTH), interpolation=cv2.INTER_AREA)[:,:,:IMG_CHANNELS] for image in x_test]\nx_test = np.array(x_test)","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"ddd5391b-c42f-4eab-a291-b95ba20047b3","_uuid":"f0ea8ea44e1a9a8dd0851959166c37cfff0551a7","trusted":false},"cell_type":"code","source":"sizes = [data['size'] for data in test_data.values()]\n\npreds = dense_model.predict(x_test/255)\npreds = [cv2.resize(image, size, interpolation=cv2.INTER_CUBIC) for image, size in zip(preds, sizes)]\npreds_masks = [label(pred > 0.5) for pred in preds]\npreds_masks = [separate_instances(label) for label in preds_masks]","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"f7cf59bf-fb3c-4ebf-84a9-b467ad195ef1","_uuid":"394e59db6c299607dd9c2c71940ea88dc5af439a","trusted":false},"cell_type":"code","source":"def 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","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"e335376b-6848-4a08-aa64-8617eeee7205","_uuid":"c181b1ef4c7d48fa2960d15767aaa168d5a35b7d","trusted":false},"cell_type":"code","source":"from collections import namedtuple\n\nMask = namedtuple('Mask', ['observation_id', 'rle'])\nrle_preds = []\n\nfor _id, preds in zip(test_data.keys(), preds_masks):\n    for pred in preds:\n        if np.sum(pred) > 10:\n            rle_preds.append(Mask(observation_id=_id, \n                                  rle=rle_encoding(pred)\n                                 )\n                            )\n\nrle_df = pd.DataFrame(rle_preds, columns=['ImageId', 'EncodedPixels'])\nrle_df.head()","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"754a725e-60cb-4bb8-bbfd-15454081184e","_uuid":"467700081afb7c16017a32957592b02ba008fc56","trusted":false},"cell_type":"code","source":"rle_df.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"facfeb04-0021-4dcc-a680-a26b4f9e625d","_uuid":"440149905aa2f709a04aa847c28ddca8183b6542"},"cell_type":"markdown","source":"# Other ideas to explore\n\n* Creating a second output channel for predicting cell contours. We can use this information to separate cells. \n* Experiment with different loss functions such as Dice and Focal loss. \n* Graduate from a Kaggle Kernel to a large machine to explore deeper models, larger batch sizes, etc. "}],"metadata":{"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"},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"}},"nbformat":4,"nbformat_minor":1}