{"cells":[{"metadata":{"_cell_guid":"ebba0809-6591-4386-8e4a-ba8d1935dd57","_uuid":"af1935f19e6ffeca47d2933c39bfcaf5a816535a"},"cell_type":"markdown","source":"This kernel will implement classical image techniques and will hopefully serve as a useful primer to people who have never worked with image data before. Ultimately, we will develop a simple pipeline using `scipy` and `numpy` (and a little bit of `scikit-image`) that we can apply to the test images -- in fact, we won't even use the training images except to optimize parameters.\n\nI'll keep updating this notebook to try and improve it - user Gabro Vecsei takes a similar approach and scores 0.22 in [this kernel]((https://www.kaggle.com/gaborvecsei/basic-pure-computer-vision-segmentation-lb-0-229). My main intention here, though, is to help out people who are new to analyses, not to score highly. "},{"metadata":{"trusted":true,"collapsed":true},"cell_type":"code","source":"import os\nimport zipfile\nimport pathlib\nimport imageio\nimport numpy as np\n\n\nzfs = {file.stem: zipfile.ZipFile(file) for file in pathlib.Path('../input/').glob('*.zip')}\nfor stem, zf in zfs.items():\n    if not stem.endswith('.csv'):\n        zf.extractall('./data/' + stem)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"2b168289-7746-4255-b933-356bd04914b7","_uuid":"25fa830192ea11807891ec2f71487dbeb18364e3","trusted":true,"collapsed":true},"cell_type":"code","source":"import pathlib\nimport imageio\nimport numpy as np\n\n# Glob the training data and load a single image path\ntraining_paths = pathlib.Path('./data/stage1_train').glob('*/images/*.png')\ntraining_sorted = sorted([x for x in training_paths])\ntraining_sorted\nim_path = training_sorted[45]\nim = imageio.imread(str(im_path))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"training_paths","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"457b5a62-435d-44c5-a7d8-30367c5212f2","_uuid":"c6854798a17802a06c246173b84848011d266fda"},"cell_type":"markdown","source":"# Dealing with color\n\nThe images in this dataset can be in RGB, RGBA and grayscale format, based on the \"modality\" in which they are acquired. For color images, there is a third dimension which encodes the \"channel\" (e.g. Red, Green, Blue). To make things simpler for this first pass, we can coerce all these images into grayscale using the `rgb2gray` function from `scikit-image`."},{"metadata":{"_cell_guid":"3f389854-b251-4723-babf-9fd0d71e3805","_uuid":"3f03d9fd9cb23494b936e2ebd4535878efc2ccfd","trusted":true},"cell_type":"code","source":"# Print the image dimensions\nprint('Original image shape: {}'.format(im.shape))\n\n# Coerce the image into grayscale format (if not already)\nfrom skimage.color import rgb2gray\nim_gray = rgb2gray(im)\nprint('New image shape: {}'.format(im_gray.shape))","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"c298bcd0-2c4b-42aa-b0d1-f5b829515841","_kg_hide-input":true,"_uuid":"b1e6a69a5ae9121ff4f2f1a4fa96698f532a9375","trusted":true},"cell_type":"code","source":"# Now, let's plot the data\nimport matplotlib.pyplot as plt\n\nplt.figure(figsize=(10,4))\n\nplt.subplot(1,2,1)\nplt.imshow(im)\nplt.axis('off')\nplt.title('Original Image')\n\nplt.subplot(1,2,2)\nplt.imshow(im_gray, cmap='gray')\nplt.axis('off')\nplt.title('Grayscale Image')\n\nplt.tight_layout()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"5f5383ed-7c97-4c27-8ea4-ac4ae35bbd25","_uuid":"0a93026712b2a2ef54412c43873ac58197636917"},"cell_type":"markdown","source":"# Removing background\nPerhaps the simplest approach for this problem is to assume that there are two classes in the image: objects of interest and the background. Under this assumption, we would expect the data to fall into a bimodal distribution of intensities. If we found the best separation value, we could \"mask\" out the background data, then simply count the objects we're left with.\n\nThe \"dumbest\" way we could find the threshold value would be to use a simple descriptive statistic, such as the mean or median. But there are other methods: the \"Otsu\" method is useful because it models the image as a bimodal distribution and finds the optimal separation value. "},{"metadata":{"_cell_guid":"f145ff54-8e4d-4546-9913-e695bb7d9254","_uuid":"b76b12a4b26dc1a351d717a047bd90ea6aa3f7dd","collapsed":true,"trusted":true},"cell_type":"code","source":"from skimage.filters import threshold_otsu\nthresh_val = threshold_otsu(im_gray)\nmask = np.where(im_gray > thresh_val, 1, 0)\n\n# Make sure the larger portion of the mask is considered background\nif np.sum(mask==0) < np.sum(mask==1):\n    mask = np.where(mask, 0, 1)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"65033fba-1681-46fd-aa6d-c3e3fec09c1f","_kg_hide-input":true,"_uuid":"14792479e2019e59f6cd0388b40798debb9e34dc","trusted":true,"collapsed":true},"cell_type":"code","source":"plt.figure(figsize=(10,4))\n\nplt.subplot(1,2,1)\nim_pixels = im_gray.flatten()\nplt.hist(im_pixels,bins=50)\nplt.vlines(thresh_val, 0, 100000, linestyle='--')\nplt.ylim([0,50000])\nplt.title('Grayscale Histogram')\n\nplt.subplot(1,2,2)\nmask_for_display = np.where(mask, mask, np.nan)\nplt.imshow(im_gray, cmap='gray')\nplt.imshow(mask_for_display, cmap='rainbow', alpha=0.5)\nplt.axis('off')\nplt.title('Image w/ Mask')\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"6b6afa37-f217-4c0b-a9d4-8e244d5384b9","_uuid":"47cc36785c751676f70dec7b25d2a4c56fe4cf0c"},"cell_type":"markdown","source":"# Deriving individual masks for each object\n\nFor this contest, we need to get a separate mask for each nucleus. One way we can do this is by looking for all objects in the mask that are connected, and assign each of them a number using `ndimage.label`.  Then, we can loop through each `label_id` and add it to an iterable, such as a list."},{"metadata":{"_cell_guid":"5c8b868f-1b2d-4939-914f-836d0aa7b50d","_uuid":"2475893dbeabff7c165283970eebea2fd2de9bd7","trusted":true,"collapsed":true},"cell_type":"code","source":"from scipy import ndimage\nlabels, nlabels = ndimage.label(mask)\n\nlabel_arrays = []\nfor label_num in range(1, nlabels+1):\n    label_mask = np.where(labels == label_num, 1, 0)\n    label_arrays.append(label_mask)\n\nprint('There are {} separate components / objects detected.'.format(nlabels))","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"14e0c916-d8c9-42df-acd0-3a89f4d0d072","_kg_hide-input":true,"_uuid":"53cb03a6457b538ad803ffe4fc1159b085f2a72e","trusted":true,"collapsed":true},"cell_type":"code","source":"# Create a random colormap\nfrom matplotlib.colors import ListedColormap\nrand_cmap = ListedColormap(np.random.rand(256,3))\n\nlabels_for_display = np.where(labels > 0, labels, np.nan)\nplt.imshow(im_gray, cmap='gray')\nplt.imshow(labels_for_display, cmap=rand_cmap)\nplt.axis('off')\nplt.title('Labeled Cells ({} Nuclei)'.format(nlabels))\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a510167c5bc06a6843a4255d8e3505df23cb95f1"},"cell_type":"markdown","source":"A quick glance reveals two problems (in this very simple image): \n\n- There are a few individual pixels that stand alone (e.g. top-right)\n- Some cells are combined into a single mask (e.g., top-middle)\n    \nUsing `ndimage.find_objects`, we can iterate through our masks, zooming in on the individual nuclei found to apply additional processing steps.  `find_objects` returns a list of the coordinate range for each labeled object in your image."},{"metadata":{"trusted":true,"_uuid":"34ab91ec5564752381fe300645aa4e4732288e14","collapsed":true},"cell_type":"code","source":"for label_ind, label_coords in enumerate(ndimage.find_objects(labels)):\n    cell = im_gray[label_coords]\n    \n    # Check if the label size is too small\n    if np.product(cell.shape) < 10: \n        print('Label {} is too small! Setting to 0.'.format(label_ind))\n        mask = np.where(labels==label_ind+1, 0, mask)\n\n# Regenerate the labels\nlabels, nlabels = ndimage.label(mask)\nprint('There are now {} separate components / objects detected.'.format(nlabels))","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"28cc9fda0ca6b5e0060591ca9ed7f41d52919a20","collapsed":true},"cell_type":"code","source":"fig, axes = plt.subplots(1,6, figsize=(10,6))\n\nfor ii, obj_indices in enumerate(ndimage.find_objects(labels)[0:6]):\n    cell = im_gray[obj_indices]\n    axes[ii].imshow(cell, cmap='gray')\n    axes[ii].axis('off')\n    axes[ii].set_title('Label #{}\\nSize: {}'.format(ii+1, cell.shape))\n\nplt.tight_layout()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c44a9f56a96eb0b65464f1b20a3738dc4051826d"},"cell_type":"markdown","source":"Label #2 has the \"adjacent cell\" problem: the two cells are being considered part of the same object. One thing we can do here is to see whether we can shrink the mask to \"open up\" the differences between the cells. This is called mask erosion. We can then re-dilate it to to recover the original proportions. "},{"metadata":{"trusted":true,"_uuid":"aea9549db2f19d7551fb2a021bd877572c8555fb","collapsed":true},"cell_type":"code","source":"# Get the object indices, and perform a binary opening procedure\ntwo_cell_indices = ndimage.find_objects(labels)[1]\ncell_mask = mask[two_cell_indices]\ncell_mask_opened = ndimage.binary_opening(cell_mask, iterations=8)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"cc2b1241e08673952bc6da9f4d05b95608a2c1e8","collapsed":true},"cell_type":"code","source":"fig, axes = plt.subplots(1,4, figsize=(12,4))\n\naxes[0].imshow(im_gray[two_cell_indices], cmap='gray')\naxes[0].set_title('Original object')\naxes[1].imshow(mask[two_cell_indices], cmap='gray')\naxes[1].set_title('Original mask')\naxes[2].imshow(cell_mask_opened, cmap='gray')\naxes[2].set_title('Opened mask')\naxes[3].imshow(im_gray[two_cell_indices]*cell_mask_opened, cmap='gray')\naxes[3].set_title('Opened object')\n\n\nfor ax in axes:\n    ax.axis('off')\nplt.tight_layout()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"ad2b865b-678c-41de-9db1-133462f20377","_uuid":"9f417451b0234f5f714be676cd5a217873b8fe98"},"cell_type":"markdown","source":"# Convert each labeled object to Run Line Encoding\nFinally, we need to encode each `label_mask` into a \"run line encoded\" string. Basically, we walk through the array, and when we find a pixel that is part of the mask, we index it and count how many subsequent pixels are also part of the mask. We repeat this each time we see new pixel start point.\n\nI found a nice function to do RLE from [Kaggle user Rakhlin's kernel](https://www.kaggle.com/rakhlin/fast-run-length-encoding-python), which I've copied here."},{"metadata":{"_cell_guid":"29ea6b76-766c-487d-8f88-d6b0e288e94f","_uuid":"4d64b21aa0f26d177863b3cae7fb91e80bc39180","trusted":true,"collapsed":true},"cell_type":"code","source":"def rle_encoding(x):\n    '''\n    x: numpy array of shape (height, width), 1 - mask, 0 - background\n    Returns run length as list\n    '''\n    dots = np.where(x.T.flatten()==1)[0] # .T sets Fortran order down-then-right\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 \" \".join([str(i) for i in run_lengths])\n\nprint('RLE Encoding for the current mask is: {}'.format(rle_encoding(label_mask)))","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"42577f67-935c-4269-b6ef-d7c166723779","_uuid":"dfcd51eaaa5fdf144ce7e9de6397e8abea0b50b9","collapsed":true},"cell_type":"markdown","source":"# Combine it into a single function\nNow that we've seen the basic steps to processing an image in a \"dumb\" way, we can combine it all into a single function. This function will take an image path, perform the processes outlined above, and spit out a dataframe with the RLE strings for each mask found. \n\nWe also create a wrapper function that will spit out a single DataFrame for all images in the dataset."},{"metadata":{"_cell_guid":"0d32b900-8f40-4eff-8465-0f715f24faf6","_uuid":"944f216004fa635201470b8eaa28a446d9679516","collapsed":true,"trusted":true},"cell_type":"code","source":"import pandas as pd\n\ndef analyze_image(im_path):\n    '''\n    Take an image_path (pathlib.Path object), preprocess and label it, extract the RLE strings \n    and dump it into a Pandas DataFrame.\n    '''\n    # Read in data and convert to grayscale\n    im_id = im_path.parts[-3]\n    im = imageio.imread(str(im_path))\n    im_gray = rgb2gray(im)\n    \n    # Mask out background and extract connected objects\n    thresh_val = threshold_otsu(im_gray)\n    mask = np.where(im_gray > thresh_val, 1, 0)\n    if np.sum(mask==0) < np.sum(mask==1):\n        mask = np.where(mask, 0, 1)    \n        labels, nlabels = ndimage.label(mask)\n    labels, nlabels = ndimage.label(mask)\n    \n    # Loop through labels and add each to a DataFrame\n    im_df = pd.DataFrame()\n    for label_num in range(1, nlabels+1):\n        label_mask = np.where(labels == label_num, 1, 0)\n        if label_mask.flatten().sum() > 10:\n            rle = rle_encoding(label_mask)\n            s = pd.Series({'ImageId': im_id, 'EncodedPixels': rle})\n            im_df = im_df.append(s, ignore_index=True)\n    \n    return im_df\n\n\ndef analyze_list_of_images(im_path_list):\n    '''\n    Takes a list of image paths (pathlib.Path objects), analyzes each,\n    and returns a submission-ready DataFrame.'''\n    all_df = pd.DataFrame()\n    for im_path in im_path_list:\n        im_df = analyze_image(im_path)\n        all_df = all_df.append(im_df, ignore_index=True)\n    \n    return all_df","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"58b07c91-855f-4e6d-a9d8-457750d9ad4a","_uuid":"b7e21d3b95d7936904bfed088a880a5c78f0fa18","collapsed":true,"trusted":true},"cell_type":"code","source":"testing = pathlib.Path('../input/stage1_test/').glob('*/images/*.png')\ndf = analyze_list_of_images(list(testing))\ndf.to_csv('submission.csv', index=None)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"51fdaa6a207574dc67d497eca81eb8d757021fe4"},"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}