getting_started/Image-Processing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Image Processing in Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Earlier, we gave a brief overview of Python and three powerful Python libraries (NumPy, MatPlotLib, and SciPy) we use when working with DeepCell. Today, we will become more familiar with SciPy and introduce a fourth package, a powerful suite of image processing tools called `scikit-image`. <i> NB: scikit-image has a description of tag `skimage` that can be used interchageably with the offical name.</i>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To install `scitkit-image`, execute one of the following command on the command line:\n",
    "\n",
    "`pip install scikit-image`\n",
    "\n",
    "`conda install scikit-image` (for Anaconda users)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial, we will cover:\n",
    "\n",
    "* Loading & Handling Image Data: Input/output, Data Types, and Colorspaces\n",
    "* Preprocessing: Contrast Adjustment, Background Subtraction, and Filtering \n",
    "* Foreground Detection: Thresholding and Morphological Operations\n",
    "* Object Detection and Segmentation: Labeling, Seeding and Expansion\n",
    "* Postprocessing: Affine Transformations\n",
    "* Writing Output to Files: Images, JSON, and NPZ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 1: Loading and Handling Image Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`skimage` has a number of useful functions we can import to load and manipulate images, as well as save the resulting images. Remember, though, this library builds on NumPy, so we will need to import that library as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from skimage.io import imread, imsave, imshow\n",
    "from skimage.transform import resize as imresize\n",
    "\n",
    "# Read an JPEG image into a numpy array\n",
    "#img = imread('resources/img_000000000_Phase_000.png')\n",
    "img = imread('resources/cat.jpg')\n",
    "print(img.dtype, img.shape)  # Prints \"uint8 (400, 248, 3)\"\n",
    "\n",
    "# We can tint the image by scaling each of the color channels\n",
    "# by a different scalar constant. The image has shape (400, 248, 3);\n",
    "# we multiply it by the array [1, 0.95, 0.9] of shape (3,);\n",
    "# numpy broadcasting means that this leaves the red channel unchanged,\n",
    "# and multiplies the green and blue channels by 0.95 and 0.9\n",
    "# respectively.\n",
    "img_tinted = img * [1, 0.95, 0.9]\n",
    "\n",
    "# Resize the tinted image to be 300 by 300 pixels.\n",
    "img_tinted = imresize(img_tinted, (300, 300))\n",
    "\n",
    "# Write the tinted image back to disk\n",
    "imsave('resources/cat_tinted.jpg', np.uint8(img_tinted))\n",
    "imshow(img_tinted.astype('uint8'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rather than \"hardcoding\" the file path, a more robust way of handling these operations is to designate a section of code to specify the directory path and file name as variables. \n",
    "\n",
    "If the file is not in the current working directory, you must also have a way of specifying the path to the directory where the file is stored. In our case, the example images are stored in the directory called 'resources' in the same folder as this notebook. Note that you can use either the full path - something like r\"/home/user/bootcamp/intro-to-deepcell/resources/example_cells_1.tif\" or the relative path, starting from the current working directory.\n",
    "\n",
    "NB: Paths and filenames can contain slashes, empty spaces and other special symbols, which can cause trouble for programming languages under certain circumstances. To circumvent such trouble, add the letter r before your string definition to create a so-called 'raw string', which is not affected by these problems (e.g. `my_raw_string = r\"some string with funny symbols: \\\\\\!/~***!\"`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a string variable with the name of the file you'd like to load (here: 'example_cells_1.tif').\n",
    "# Suggested name for the variable: filename\n",
    "filename = r'img_000000000_FITC_001.png'\n",
    "\n",
    "# Create a string variable with the path to the directory that contains the file you'd like to load.\n",
    "# Suggested name for the variable: dirpath\n",
    "dirpath = r'resources'  # Relative path\n",
    "#dirpath = r'/home/user/bootcamp/intro-to-deepcell/resources/img_000000000_FITC_001.png'  # Absolute path "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now combine the directory path and file name into one variable, the file path\n",
    "\n",
    "# Import the function 'join' from the module 'os.path'\n",
    "# This function automatically takes care of the slashes that need to be added when combining two paths.\n",
    "from os.path import join\n",
    "\n",
    "# Print the result to see that everything is correct\n",
    "# Suggested name for the variable: filepath\n",
    "filepath = join(dirpath, filename)\n",
    "print(filepath)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the file path established, we can load the image (using the `imread` function we imported earlier), make sure the load was successful, and display the image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load 'img_000000000_Far-red_001.png' and store it in a variable.\n",
    "# Suggested name for the variable: img\n",
    "img = imread(filepath)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check that the load went as expected\n",
    "\n",
    "# Check that 'img' is a variable of type 'ndarray' - use Python's built-in function 'type'.\n",
    "print(\"Loaded array is of type:\", type(img).__name__)\n",
    "\n",
    "# Print the shape of the array using the numpy-function 'shape'. \n",
    "print(\"Loaded array has shape:\", img.shape)\n",
    "\n",
    "# Check the datatype of the individual numbers in the array. You can use the array attribute 'dtype' to do so.\n",
    "print(\"Loaded values are of type:\", img.dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NB: The dtype should be 'uint16', because these are unsigned 16-bit integer images. Another common dtype for images is uint8. You can read more about the differences [here](https://www.mathworks.com/help/matlab/creating_plots/working-with-8-bit-and-16-bit-images.html) and [here](https://printaura.com/8-bit-vs-16-bit-images-whats-the-difference-which-to-use/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we're ready to take a look at the images. To plot the array as an image, use pyplot's functions `plt.imshow` followed by `plt.show`. \n",
    "\n",
    "You can check the documentation for `plt.imshow` and note the parameters that can be specified, such as colormap (cmap)\n",
    "and interpolation. Since we are working with scientific data, interpolation is unwelcome, so you should set it to \"none\". The most common cmap for grayscale images is naturally \"gray\". You may also want to adjust the size of the figure. You can do this by preparing the figure canvas with the function `plt.figure` before calling `plt.imshow`. The canvas size is adjusted using the keyword argument 'figsize' when calling `plt.figure`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.figure(figsize=(7, 7))\n",
    "plt.imshow(img, interpolation='none', cmap='gray')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Don't worry if it's dark, we'll fix that in a minute. In the meantime, for our peace of mind, here is a side-by-side example using our cat image from ealier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img2 = imread('resources/cat.jpg')\n",
    "img_tinted2 = img2 * [1, 0.95, 0.9]\n",
    "\n",
    "# Show the original image\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.imshow(img2)\n",
    "\n",
    "# Show the tinted image\n",
    "plt.subplot(1, 2, 2)\n",
    "\n",
    "# A slight gotcha with imshow is that it might give strange results if presented with data that is not uint8. \n",
    "# To work around this, we explicitly cast the image to uint8 before displaying it.\n",
    "plt.imshow(np.uint8(img_tinted2))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 2: Preprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Raw microscope images, like the ones we will be using to train our neural nets, often start out looking like garbage. Viz.:\n",
    "![alt text][phase000]\n",
    "\n",
    "[phase000]: ./resources/img_000000000_Phase_000.png \"A Bunch of Cells?\"\n",
    "\n",
    "Before we can really use this image for anything, we need to improve it. This is the (pre)process of preprocessing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Contrast Adjustment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Contrast adjustment is the process of remapping pixel intensity values in an image such that, generally, there is a broader range of lights and darks, making objects in the image more discernible.\n",
    "\n",
    "To make the process more quantitative (and, thus, easier to analyze), we often refer to the histogram of pixel intensities in an image. We will often overlay our histograms with cumulative distribution functions (CDF), which provide another quantitative way of visualizing our pixel distributions. Don't worry if you know nothing about cumulative distribution functions; we'll be able to talk about them when they come up.\n",
    "\n",
    "(**Bonus question**: can anyone tell what data type these images are encoded as, assuming their histograms are showing the entire range of possible pixel values? I.e., if we were to check the `.dtype` of these images, what would we see?)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Within `skimage`, there are multiple functions inside the `skimage.exposure` module for contrast adjustment and histogram observation. We will highlight a few of them here:\n",
    "\n",
    "* `skimage.exposure.histogram`\n",
    "* `skimage.exposure.equalize_hist`\n",
    "* `skimage.exposure.equalize_adapthist`\n",
    "* `skimage.exposure.rescale_intensity`\n",
    "* `skimage.exposure.adjust_gamma`\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **`skimage.exposure.histogram`**\n",
    "\n",
    "Let's use `skimage.exposure.histogram` to view the distribution of pixel intensities in the cat image.\n",
    "\n",
    "Note the use of indices on the cat image. It's often most informative to look at each channel in a color image separately, which is what we're doing here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.exposure import histogram\n",
    "\n",
    "# compute histogram in Red channel of cat image\n",
    "(counts, bins) = histogram(img2[:, :, 0])\n",
    "# view histogram\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.scatter(bins,counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Very cool. What do you all see here?\n",
    "\n",
    "Now, let's look at our cell image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute cell histogram\n",
    "#phase_image = imread(\"resources/img_000000000_Phase_000.png\")\n",
    "#phase_image_float = phase_image/phase_image.max()\n",
    "#(counts, bins) = histogram(phase_image_float)\n",
    "\n",
    "(counts, bins) = histogram(img)\n",
    "# view cell histogram\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.scatter(bins, counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that we have a lot of pixels with intermediate values and few (none?) the further we go out to the extremes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **`skimage.exposure.rescale_intensity`**\n",
    "\n",
    "A really basic way to adjust the contrast of an image is to stretch the min and max values found in the image to the smallest and \n",
    "largest possible intensity values, respectively, and then spread out all the intermediate values accordingly. That's what `skimage.exposure.rescale_intensity` does."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First check the cell image's min and max values\n",
    "print( \"Minimum pixel value: \", img.min())\n",
    "print( \"Maximum pixel value: \", img.max())\n",
    "\n",
    "# Also, check its dtype attribute to understand what its range of values is.\n",
    "print( \"Data type of image: \", img.dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now it's calculation time. Whoo!\n",
    "\n",
    "We can see, from looking at the data type of the image, that the minimum and maximum possible pixel values are _____ and _____, respectively.\n",
    "\n",
    "Since the observed pixel values run from 71 to 12342, we can see that our image is currently using a pixel intensity range of (12342 - 71 = 12271), which is _____% of the maximum possible intensity range of ______.\n",
    "\n",
    "This tells us that, if we just rescale the min and max values of the image to the max and min of their possible range, we will see a _____-fold increase in the pixel intensity range of the image.\n",
    "\n",
    "Let's see how much that helps now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import exposure\n",
    "from skimage.exposure import rescale_intensity\n",
    "\n",
    "# Rename initial cell image for conceptual clarity\n",
    "img_fitc = img\n",
    "\n",
    "#img_fitc = (((img_fitc-img_fitc.min())/(img_fitc.max()-img_fitc.min()))*255).astype(np.uint8)\n",
    "\n",
    "# adjust image contrast by rescaling pixel values\n",
    "img_rescaled = rescale_intensity(img_fitc)\n",
    "\n",
    "# compute histograms\n",
    "counts_fitc, bins_fitc = histogram(img_fitc)\n",
    "counts_rescaled, bins_rescaled = histogram(img_rescaled)\n",
    "\n",
    "# compute cumulative distribution functions\n",
    "img_cdf_fitc, bins_fitc = skimage.exposure.cumulative_distribution(img_fitc, 256)\n",
    "img_cdf_fitc = img_cdf_fitc\n",
    "img_cdf_rescaled, bins_rescaled = skimage.exposure.cumulative_distribution(img_rescaled, 256)\n",
    "img_cdf_rescaled = img_cdf_rescaled\n",
    "\n",
    "# view images and histograms\n",
    "fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(16, 16))\n",
    "\n",
    "## original image\n",
    "axes[0, 0].imshow(img_fitc, cmap='gray')\n",
    "\n",
    "## rescaled image\n",
    "axes[0, 1].imshow(img_rescaled, cmap='gray')\n",
    "\n",
    "## original histogram and cdf\n",
    "right_y_axis_0 = axes[1, 0].twinx()\n",
    "axes[1, 0].scatter(bins_fitc,counts_fitc)\n",
    "right_y_axis_0.set_ylim(0, 1)\n",
    "right_y_axis_0.plot(bins_fitc, img_cdf_fitc, 'r')\n",
    "\n",
    "## rescaled histogram and cdf\n",
    "right_y_axis_1 = axes[1, 1].twinx()\n",
    "axes[1, 1].scatter(bins_rescaled,counts_rescaled)\n",
    "right_y_axis_1.set_ylim(0, 1)\n",
    "right_y_axis_1.plot(bins_rescaled, img_cdf_rescaled, 'r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Any thoughts?\n",
    "\n",
    "What do we see if we look at the `min()` and `max()` of the rescaled image?\n",
    "\n",
    "How useful is this technique for this particular image?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **`skimage.exposure.equalize_hist`**\n",
    "\n",
    "Now, let's try a more sophisticated technique for contrast adjustment.\n",
    "\n",
    "Let's use `skimage.exposure.equalize_hist` to not only spread the range of pixel intensities out to its max, but to then move pixel values around so that the pixels more or less evenly occupy every possible value across this range. We won't be changing the relative intensities of any two pixels (if one pixel has a higher value than another before contrast adjustment, it will still have a higher value after), but we will be moving their absolute values around with that constraint in place."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.exposure import equalize_hist\n",
    "\n",
    "# adjust image contrast by equalizing the histogram\n",
    "img_equalized = equalize_hist(img_fitc)\n",
    "\n",
    "# compute histograms\n",
    "counts_fitc, bins_fitc = histogram(img_fitc)\n",
    "counts_equalized, bins_equalized = histogram(img_equalized)\n",
    "\n",
    "# compute cumulative distribution functions\n",
    "img_cdf_fitc, bins_fitc = exposure.cumulative_distribution(img_fitc, 256)\n",
    "img_cdf_fitc = img_cdf_fitc\n",
    "img_cdf_equalized, bins_equalized = exposure.cumulative_distribution(img_equalized, 256)\n",
    "img_cdf_equalized = img_cdf_equalized\n",
    "\n",
    "# view images and histograms\n",
    "fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(16,16))\n",
    "\n",
    "## original image\n",
    "axes[0, 0].imshow(img_fitc, cmap='gray')\n",
    "\n",
    "## rescaled image\n",
    "axes[0, 1].imshow(img_equalized, cmap='gray')\n",
    "\n",
    "## original histogram and cdf\n",
    "right_y_axis_0 = axes[1, 0].twinx()\n",
    "axes[1, 0].scatter(bins_fitc,counts_fitc)\n",
    "right_y_axis_0.set_ylim(0, 1)\n",
    "right_y_axis_0.plot(bins_fitc, img_cdf_fitc, 'r')\n",
    "\n",
    "## rescaled histogram and cdf\n",
    "right_y_axis_1 = axes[1, 1].twinx()\n",
    "axes[1, 1].scatter(bins_equalized,counts_equalized)\n",
    "right_y_axis_1.set_ylim(0, 1)\n",
    "right_y_axis_1.plot(bins_equalized, img_cdf_equalized, 'r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we have we moved from a uint16 pixel encoding to some sort of float encoding (probably float64) in the right image, due to the way `skimage.exposure.equalize_hist` operates. How can we tell that our data type has changed?\n",
    "\n",
    "Now, if we compare these two images' histograms, we see that the distribution of the pixels has been dramatically transformed!\n",
    "\n",
    "What have we done here? Let's discuss it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The net takeaway is that we've transformed the cumulative distribution function (CDF) of pixel values so that the pixels are now optimally spaced across the entire pixel intensity range. This maximizes the global contrast of the image, which generally makes things easier to see.\n",
    "\n",
    "To read more about the math of histogram equalization (which is pleasantly intuitive, compared to a lot of other stuff), visit  \n",
    "https://en.wikipedia.org/wiki/Histogram_equalization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **`skimage.exposure.equalize_adapthist`**\n",
    "\n",
    "As a more robust alternative to `skimage.exposure.equalize_hist`, we can use `skimage.exposure.equalize_adapthist`, which performs the same process of spreading pixel intensity values evenly over the available range and then remapping them to smooth out the histogram. However, whereas `skimage.exposure.equalize_hist` does this for the entire image at once, `skimage.exposure.equalize_adapthist` looks at each pixel within the context of that pixel's local neighborhood only. This is particularly useful for images with strong global contrasts between lights and darks across the entire image, but minimal local contrasts between lights and darks.\n",
    "\n",
    "(By default, the neighborhoods are 1/8x1/8 of the image. This neighborhood size can be changed by setting the `kernel_size` argument in `skimage.exposure.equalize_adapthist` to the desired pixel value of the neighborhood. Give it a try.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.exposure import equalize_adapthist\n",
    "\n",
    "# adjust image contrast by adaptively equalizing the histogram\n",
    "img_equalized_adapted = equalize_adapthist(img_fitc, kernel_size=None)\n",
    "\n",
    "# compute histograms\n",
    "counts_fitc, bins_fitc = histogram(img_fitc)\n",
    "counts_equalized_adapted, bins_equalized_adapted = histogram(img_equalized_adapted)\n",
    "\n",
    "# compute cumulative distribution functions\n",
    "img_cdf_fitc, bins_fitc = exposure.cumulative_distribution(img_fitc, 256)\n",
    "img_cdf_fitc = img_cdf_fitc\n",
    "img_cdf_equalized_adapted, bins_equalized_adapted = exposure.cumulative_distribution(img_equalized_adapted, 256)\n",
    "img_cdf_equalized_adapted = img_cdf_equalized_adapted\n",
    "\n",
    "# view images and histograms\n",
    "fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(16, 16))\n",
    "\n",
    "## original image\n",
    "axes[0, 0].imshow(img_fitc, cmap='gray')\n",
    "\n",
    "## rescaled image\n",
    "axes[0, 1].imshow(img_equalized_adapted, cmap='gray')\n",
    "\n",
    "## original histogram and cdf\n",
    "right_y_axis_0 = axes[1, 0].twinx()\n",
    "axes[1, 0].scatter(bins_fitc,counts_fitc)\n",
    "right_y_axis_0.set_ylim(0,1)\n",
    "right_y_axis_0.plot(bins_fitc, img_cdf_fitc, 'r')\n",
    "\n",
    "## rescaled histogram and cdf\n",
    "right_y_axis_1 = axes[1,1].twinx()\n",
    "axes[1,1].scatter(bins_equalized_adapted,counts_equalized_adapted)\n",
    "right_y_axis_1.set_ylim(0,1)\n",
    "right_y_axis_1.plot(bins_equalized_adapted, img_cdf_equalized_adapted, 'r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The net result of `skimage.exposure.equalize_adapthist` is to spread out really common pixel values somewhat, thereby increasing contrast among common pixel values, but not dramatically altering many values.\n",
    "\n",
    "Mathematical details can be found here:\n",
    "[https://en.wikipedia.org/wiki/CLAHE](https://en.wikipedia.org/wiki/CLAHE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **`skimage.exposure.adjust_gamma`**\n",
    "\n",
    "`skimage.exposure.adjust_gamma` spread pixel intensities out, like all the other contrast adjustment algorithms. After that, though, it normalizes all the pixel values so that they range from 0 to 1 (like it appears that the other `skimage` contrast adjustment functions are doing) and performs Gamme Correction, which just raises every pixel value to a certain exponent. Note that often this exponent is less than one. Since we normalized pixel values ahead of time, this transform doesn't cause us to sacrifice any of our pixel intensity range, which is part of what makes it appealing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.exposure import adjust_gamma\n",
    "\n",
    "# adjust image contrast according to the Gamma Correction algorithm\n",
    "img_gamma = adjust_gamma(img_fitc, gamma=1, gain=1)\n",
    "\n",
    "# compute histograms\n",
    "counts_fitc, bins_fitc = histogram(img_fitc)\n",
    "counts_gamma, bins_gamma = histogram(img_gamma)\n",
    "\n",
    "# compute cumulative distribution functions\n",
    "img_cdf_fitc, bins_fitc = exposure.cumulative_distribution(img_fitc, 256)\n",
    "img_cdf_fitc = img_cdf_fitc\n",
    "img_cdf_gamma, bins_gamma = exposure.cumulative_distribution(img_gamma, 256)\n",
    "img_cdf_gamma = img_cdf_gamma\n",
    "\n",
    "# view images and histograms\n",
    "fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(16, 16))\n",
    "\n",
    "## original image\n",
    "axes[0, 0].imshow(img_fitc, cmap='gray')\n",
    "\n",
    "## rescaled image\n",
    "axes[0, 1].imshow(img_gamma, cmap='gray')\n",
    "\n",
    "## original histogram and cdf\n",
    "right_y_axis_0 = axes[1,0].twinx()\n",
    "axes[1,0].scatter(bins_fitc,counts_fitc)\n",
    "right_y_axis_0.set_ylim(0,1)\n",
    "right_y_axis_0.plot(bins_fitc, img_cdf_fitc, 'r')\n",
    "\n",
    "## rescaled histogram and cdf\n",
    "right_y_axis_1 = axes[1, 1].twinx()\n",
    "axes[1,1].scatter(bins_gamma,counts_gamma)\n",
    "right_y_axis_1.set_ylim(0,1)\n",
    "right_y_axis_1.plot(bins_gamma, img_cdf_gamma, 'r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Summary**\n",
    "\n",
    "In this example only one of the contrast algorithms really worked well. In practice, maybe multiple of these will work well, maybe none will. There is a lot of trial and error in contrast adjustment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filtering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start simple with an application of Gaussian smoothing. To do so, we will use the Gaussian filter function `ndi.filters.gaussian_filter` from the image processing module `scipy.ndimage`. Make sure and review the SciPy documentation to see how to use this function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import\n",
    "from skimage.filters.thresholding import threshold_otsu\n",
    "\n",
    "# Import the image processing package scipy.ndimage as ndi\n",
    "import scipy.ndimage as ndi\n",
    "\n",
    "# The documentation tells us that the gaussian_filter function expects a smoothing factor sigma, \n",
    "# so we will arbitrarily define one (this can be changed later)\n",
    "sigma = 4\n",
    "\n",
    "# Apply the filter and allocate the output to a new variable.\n",
    "img_smooth = ndi.filters.gaussian_filter(img, sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize the result using 'plt.imshow'\n",
    "\n",
    "# Compare with the original image visualized above. Can you optimize sigma such that the image looks \n",
    "# smooth without blurring the membranes too much?\n",
    "plt.figure(figsize=(7,7))\n",
    "plt.imshow(img_smooth, interpolation='none', cmap='gray')\n",
    "plt.show()\n",
    "\n",
    "# To have a closer look at a specific region of the image, crop that region out and show it in a \n",
    "# separate plot. Remember that you can crop arrays by \"indexing\" or \"slicing\" them similar to lists.\n",
    "plt.figure(figsize=(6,6))\n",
    "plt.imshow(img_smooth[400:600, 200:400], interpolation='none', cmap='gray')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Look at the raw and smoothed images side by side using 'plt.subplots'\n",
    "fig, ax = plt.subplots(1, 2, figsize=(10,7))\n",
    "ax[0].imshow(img, interpolation='none', cmap='gray')\n",
    "ax[1].imshow(img_smooth, interpolation='none', cmap='gray')\n",
    "ax[0].set_title('Raw Image')\n",
    "ax[1].set_title('Smoothed Image')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Homework**  \n",
    "Try performing Canny Edge Detection tonight on your own. To do this, use the `skimage.feature.canny` function and feed it your smoothed image, `img_smooth`.\n",
    "\n",
    "If you can't get it to work, don't stress. We'll go over this tomorrow morning."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 3: Foreground Detection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Manual Thresholding & Threshold Detection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a variable (int) to hold a threshold value, which can be changed later to something more suitable.\n",
    "thresh = 700\n",
    "\n",
    "# Recall that relational (Boolean) expressions, such as 'smaller' (<), 'equal' (==) or 'greater or equal' (>=),\n",
    "# can be used with numpy arrays to directly assign the result to a new variable.\n",
    "mem = img_smooth > thresh\n",
    "\n",
    "# Check the dtype of your thresholded image - it should be boolean, meaning an array filled with 'True' and 'False',\n",
    "# where 'True' is the foreground (regions above the threshold) and 'False' is the background.\n",
    "print(mem.dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at the result. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(10, 7))\n",
    "ax[0].imshow(img_smooth, interpolation='none', cmap='gray')\n",
    "ax[1].imshow(mem, interpolation='none', cmap='gray')\n",
    "ax[0].set_title('Original Image')\n",
    "ax[1].set_title('Thresholded Membranes')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can iterate through various threshold values to select something more appropriate. To do so interactively, we will utilize a class of interactive functions called 'widgets.' These are incredibly useful in exploratory data analysis to create simplified 'User Interfaces' (UIs) on the fly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare widget\n",
    "from ipywidgets import interact\n",
    "\n",
    "@interact(thresh=(100,1500,50))\n",
    "def select_threshold(thresh=200):\n",
    "    \n",
    "    # Thresholding\n",
    "    mem = img_smooth > thresh\n",
    "    \n",
    "    # Visualization\n",
    "    plt.figure(figsize=(7,7))\n",
    "    plt.imshow(mem, interpolation='none', cmap='gray')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's say we get tired of searching for the perfect threshold by hand (regardless of how cool our widget is). The scikit-image module `skimage.filters.thresholding` provides several threshold detection algorithms. One of the most popular ones \n",
    "is [Otsu's method](https://en.wikipedia.org/wiki/Otsu%27s_method). We can import the module and use it to automatically \n",
    "determine a threshold for the smoothed image. Then we can apply the threshold and visualize the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import\n",
    "from skimage.filters.thresholding import threshold_otsu\n",
    "\n",
    "# Import the image processing package scipy.ndimage as ndi\n",
    "import scipy.ndimage as ndi\n",
    "\n",
    "# the gaussian_filter function expects a smoothing factor sigma, so we will arbitrarily define one (this can change later)\n",
    "sigma = 4\n",
    "img_smooth = ndi.filters.gaussian_filter(img, sigma)\n",
    "\n",
    "# Calculate and apply threshold\n",
    "thresh = threshold_otsu(img)\n",
    "mem = img > thresh\n",
    "    \n",
    "# Visualization\n",
    "plt.figure(figsize=(7, 7))\n",
    "plt.imshow(mem, interpolation='none', cmap='gray')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wait, what? Well, as it turns out, Otsu's method works best on smoothed images. Fortunately, scipy can help us there too. We can import the multi-dimensional image processing package from scipy and apply a Guassian filter to smooth the image before re-running Otsu's method. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import\n",
    "from skimage.filters.thresholding import threshold_otsu\n",
    "\n",
    "# Import the image processing package scipy.ndimage as ndi\n",
    "import scipy.ndimage as ndi\n",
    "\n",
    "# the gaussian_filter function expects a smoothing factor sigma, so we will arbitrarily define one (this can change later)\n",
    "sigma = 4\n",
    "img_smooth = ndi.filters.gaussian_filter(img, sigma)\n",
    "\n",
    "# Calculate and apply threshold\n",
    "thresh = threshold_otsu(img_smooth)\n",
    "mem = img_smooth > thresh\n",
    "    \n",
    "# Visualization\n",
    "plt.figure(figsize=(7,7))\n",
    "plt.imshow(mem, interpolation='none', cmap='gray')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example highlights a few important points: (1) Python (and the packages available for it) are very powerful in their ability to solve difficult problems quickly in very few lines of code; and (2) this can be dangerous - with great power comes great responsibility -> it is up to you to find these modules and understand what sort of data they are expecting.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Morphological Operations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Structuring Elements and Adaptive Thresholding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our earlier section on thresholding, we quickly discovered the limitations of naive thresholding by a fixed value accross an entire image. Adaptive thresholding methods improve on this approach by varying the threshold for certain parts of the image.  \n",
    "\n",
    "One way of accomplishing this is known as adpative background subtraction. In this approach, we use a strong \"mean filter\" with a circular [structuring element](https://en.wikipedia.org/wiki/Structuring_element) to create the background image. Then, we use this background image for thresholding. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a circular structuring element (SE) whose size depends on i\n",
    "i = 75\n",
    "SE = (np.mgrid[:i,:i][0] - np.floor(i/2))**2 + (np.mgrid[:i,:i][1] - np.floor(i/2))**2 <= np.floor(i/2)**2\n",
    "\n",
    "# Visualize the result\n",
    "plt.imshow(SE, interpolation='none')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the background by running a mean filter over the image using the disc SE and assign the output to a new variable\n",
    "# Use the function 'skimage.filters.rank.mean'\n",
    "from skimage.filters import rank \n",
    "\n",
    "# Test\n",
    "img8 = (((img - img.min()) / (img.max() - img.min())) * 255).astype('uint8')\n",
    "img_smooth8 = ndi.filters.gaussian_filter(img8, sigma)\n",
    "\n",
    "bg8 = rank.mean(img_smooth8, selem=SE)\n",
    "# End Test\n",
    "\n",
    "#bg = rank.mean(img_smooth, selem=SE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize the resulting background image. \n",
    "plt.figure(figsize=(7, 7))\n",
    "#plt.imshow(bg, interpolation='none', cmap='gray')\n",
    "plt.imshow(bg8, interpolation='none', cmap='gray')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Threshold the Gaussian-smoothed original image against the background image using a relational expression\n",
    "#mem = img_smooth > bg\n",
    "mem8 = img_smooth8 > bg8\n",
    "\n",
    "# Visualize the result \n",
    "fig, ax = plt.subplots(1, 2, figsize=(12, 12))\n",
    "ax[0].imshow(img_smooth, interpolation='none', cmap='gray')\n",
    "ax[1].imshow(mem8, interpolation='none', cmap='gray')\n",
    "ax[0].set_title('Original Image')\n",
    "ax[1].set_title('Adaptively Thresholded Cells')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Improving Masks with Binary Morphology"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While it is an improvment on our original results, our latest segmentation still leaves something to be desired. Let's see if we can improve the membrane segmentation with morphological operations. \n",
    "\n",
    "There are a number of operations available in the `ndimage` module that may be of use here (`ndi.binary_closing` or `ndi.binary_erosion` for example). Feel free to explore and see how the different functions affect the mask. Can you optimize the mask, for example by closing gaps?\n",
    "\n",
    "Also, note that the default SE for these functions is a square. Try to create another disc-shaped SE and see how that changes the outcome. Keep in mind that, for some funcitons, the documentation points out that the mask should be inverted. The function `np.logical_not` or the corresponding operator '~' could help there.\n",
    "Usage: \n",
    "    #mem_holefilled = np.logical_not(ndi.binary_fill_holes(np.logical_not(mem)))  # Long form\n",
    "    #mem_holefilled = ~ndi.binary_fill_holes(mem)                                 # Short form\n",
    "\n",
    "In our first attempt, we are going to attempt to improve the mask using binary erosion. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "mem_default = ndi.binary_opening(mem8)  # Short form\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# New circular SE of appropriate size \n",
    "i = 17\n",
    "SE = (np.mgrid[:i,:i][0] - np.floor(i/2))**2 + (np.mgrid[:i,:i][1] - np.floor(i/2))**2 <= np.floor(i/2)**2\n",
    "\n",
    "# Sometimes, SE can introduce boundary artifact issues that can be solved with 'padding.' 'Padding' refers to the \n",
    "# extension of the image at the boundaries, in this case using a 'reflection' of the pixel values next to the boundary. \n",
    "# If morphological operations are done on the padded image, the boundary artifacts will occur in the padded region\n",
    "# outside the original image, which can simply be cropped out again at the end.\n",
    "pad_size = i+1\n",
    "mem_padded = np.pad(mem8, pad_size, mode='reflect')\n",
    "\n",
    "# Binary closing works well to round off the membranes and close gaps\n",
    "mem_final = ndi.binary_opening(mem_padded, structure=SE)\n",
    "\n",
    "# This slicing operation crops the padded image back to the original size\n",
    "mem_final = mem_final[pad_size:-pad_size, pad_size:-pad_size]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize the final result\n",
    "fig, ax = plt.subplots(1, 2, figsize=(15, 15))\n",
    "ax[0].imshow(mem_default, interpolation='none', cmap='gray')\n",
    "ax[1].imshow(mem_final, interpolation='none', cmap='gray')\n",
    "ax[0].set_title('Default')\n",
    "ax[1].set_title('Circular SE with Padding')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 4: Object Detection and Segmentation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Connected Components Label"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the improved mask for the labeling of the connected components (cells).\n",
    "\n",
    "This is a straight-forward task thanks to the `ndi.label` function available in the `ndimage` module. Note that this function labels foreground pixels (1s, not 0s), so certain masks may need to be inverted. Also, note that 'ndi.label' returns another result in addition to the labeled image. Read up on this in the function's documention to make sure you don't mix up the two outputs.\n",
    "\n",
    "For functions with multiple outputs (in this case the labeled image and the number of detected objects), it is convention to unpack those outputs that will not be used in the remainder of code into the variable '_' (underscore). This makes it clear to those reading the code that the function returns multiple things but some of them are not important to this particular use case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell_labels, _ = ndi.label(mem_final)\n",
    "\n",
    "# Visualize the output\n",
    "\n",
    "# NB: It is no longer ideal to use a 'gray' colormap, since we want to visualize that each\n",
    "# cell has a unique ID. Explore various colormaps (check the docs to see what types of colormaps \n",
    "# are available) and choose one that works well.\n",
    "\n",
    "plt.figure(figsize=(10, 10))\n",
    "plt.imshow(cell_labels, interpolation='none', cmap='inferno')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cell Segmentation by Seeding & Expansion via the Watershed Transform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our results look much better, but we are still making some mistakes. The most egregious problem we need to address is when two cells are touching. We need to be able to distinguish where one cell ends and another begins. Let's see if we can use the watershed transform to help. \n",
    "\n",
    "We begin by running a distance transform to find seeds for our watershed transform. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the function 'ndi.distance_transform_edt'.\n",
    "dist_trans = ndi.distance_transform_edt(mem_final)\n",
    "\n",
    "# Smooth the distance transform\n",
    "dist_trans_smooth = ndi.filters.gaussian_filter(dist_trans, sigma=3)\n",
    "\n",
    "# Visualize the result\n",
    "plt.figure(figsize=(7, 7))\n",
    "plt.imshow(dist_trans_smooth, interpolation='none', cmap='viridis')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we need to retrieve the local maxima (or 'peaks') from the distance transform. To do so, we will utilize the `peak_local_max` function from the `skimage.feature` module.\n",
    "\n",
    "By default, this function returns the indices of the pixels where the local maxima are located. However, we need a boolean mast of the same shape as the original image, where the local maximum pixels are labeled as '1' and everything else is '0'. Fortunately, the function has a keyword argument for this (set indices = False)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.feature import peak_local_max\n",
    "seeds = peak_local_max(dist_trans_smooth, indices=False, min_distance=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we were to plot only the local maxima, it would just look like a bunch of distributed dots. To get an idea if the seeds are well-placed, we need to overlay these dots onto the original image.\n",
    "\n",
    "To do this, it is important to first understand a key point about how the 'pyplot' module works: the output of every plotting command is layered on top of the previous result, until everything is displayed at once when 'plt.show' is called. Therefore, we can first plot the raw (or smoothed) input image and then plot the seeds on top of it before showing both with 'plt.show'.\n",
    "\n",
    "Unfortunately, there is one other issue: the zero values in the seed array are painted in black and cover the image we hoped to use as the background. To solve this problem, you need to mask these zero values before plotting the seeds. You can do this by creating an appropriately masked array using the function 'np.ma.array' with the keyword argument 'mask'. Check the documentation or Stack Overflow to figure out how to do this.\n",
    "\n",
    "You may find it helpful to use 'ndi.filters.maximum_filter' to dilate the seeds a little bit, making them bigger and thus better visible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dilate seeds\n",
    "seeds_dil = ndi.filters.maximum_filter(seeds, size=10)\n",
    "\n",
    "# Visualize the output as an overlay on the raw (or smoothed) image\n",
    "plt.figure(figsize=(10, 10))\n",
    "plt.imshow(img_smooth, interpolation='none', cmap='gray')\n",
    "plt.imshow(np.ma.array(seeds_dil, mask=seeds_dil==0), interpolation='none', cmap='autumn')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Does this result make sense? Are all of the seeds actually representative of cells? If not, try adjusting the sigma for the smoothing filter or minimum distance on the peak local max function.\n",
    "\n",
    "If everything looks good, let's make each of the labels unique so that each cell has an ID"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use connected component labeling to give each cell seed a unique ID number.\n",
    "seeds_labeled = ndi.label(seeds)[0]\n",
    "\n",
    "# Visualize the final result (the labeled seeds) as an overlay on the raw (or smoothed) image\n",
    "seeds_labeled_dil = ndi.filters.maximum_filter(seeds_labeled, size=10)  # Expand a bit for visualization\n",
    "plt.figure(figsize=(10, 10))\n",
    "plt.imshow(img_smooth, interpolation='none', cmap='gray')\n",
    "plt.imshow(np.ma.array(seeds_labeled_dil, mask=seeds_labeled_dil==0), interpolation='none', cmap='prism')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Expansion via the Watershed Transform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are finally ready to perform the watershed transform using the function `watershed` from the module `skimage.morphology`. We will use the labeled cell seeds and the smoothed membrane image as our input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.morphology import watershed\n",
    "#ws = watershed(img_smooth8, seeds_labeled)\n",
    "#ws = watershed(~img_smooth8, seeds_labeled)\n",
    "ws = watershed(-img_smooth, seeds_labeled)\n",
    "\n",
    "# Show the result as transparent overlay over the smoothed input image\n",
    "\n",
    "# Like the masked overlay of the seeds, this can be achieved by making two calls to 'imshow',\n",
    "# one for the background image and one for the segmentation. Instead of masking away background,\n",
    "# this time we make the segmentation image semi-transparent by adjusting the keyword argument 'alpha' \n",
    "# of the 'imshow' function, which specifies opacity. Be sure to choose an appropriate colormap that \n",
    "# allows you to distinguish the segmented cells even if cells with a very similar ID are next to each other\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.imshow(img_smooth8, interpolation='none', cmap='gray')\n",
    "plt.imshow(ws, interpolation='none', cmap='prism', alpha=0.4)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is this what we were expecting? What about the background? \n",
    "\n",
    "Review the toy watershed example below (adapted from the `skimage` [documentation](http://scikit-image.org/docs/dev/api/skimage.morphology.html#skimage.morphology.watershed) and see if you can spot the problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an image with two overlapping circles\n",
    "x, y = np.indices((80, 80))\n",
    "x1, y1, x2, y2 = 28, 28, 44, 52\n",
    "r1, r2 = 16, 20\n",
    "mask_circle1 = (x - x1)**2 + (y - y1)**2 < r1**2\n",
    "mask_circle2 = (x - x2)**2 + (y - y2)**2 < r2**2\n",
    "image = np.logical_or(mask_circle1, mask_circle2)\n",
    "\n",
    "# Run the distance transformation and seeding\n",
    "from scipy import ndimage as ndi\n",
    "distance = ndi.distance_transform_edt(image)\n",
    "from skimage.feature import peak_local_max\n",
    "local_maxi = peak_local_max(distance, labels=image, footprint=np.ones((3, 3)), indices=False)\n",
    "markers = ndi.label(local_maxi)[0]\n",
    "\n",
    "# Run the watershed transform\n",
    "labels = watershed(-distance, markers, mask=image)\n",
    "#labels = watershed(-distance, markers)\n",
    "\n",
    "# Visualize\n",
    "fig, ax = plt.subplots(1, 2, figsize=(10, 10))\n",
    "ax[0].imshow(image, interpolation='none', cmap='gray')\n",
    "ax[1].imshow(labels, interpolation='none', cmap='gray')\n",
    "ax[0].set_title('Original Image')\n",
    "ax[1].set_title('Watershed Labels')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rerun the transform \n",
    "#ws = watershed(img_smooth8, seeds_labeled)\n",
    "#ws = watershed(~img_smooth8, seeds_labeled)\n",
    "ws = watershed(-img_smooth8, seeds_labeled, mask=mem_final)\n",
    "\n",
    "# Show the result as transparent overlay over the smoothed input image\n",
    "plt.figure(figsize=(10, 10))\n",
    "plt.imshow(img_smooth8, interpolation='none', cmap='gray')\n",
    "plt.imshow(ws, interpolation='none', cmap='prism', alpha=0.9)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 5: Postprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Affine Transformations\n",
    "\n",
    "Affine transformations are, intuitively, transformations of an image that preserve straight lines (so no wavy modifiers) and preserve relative distances along a line (so no stretching out of one half of an image).\n",
    "\n",
    "Here is an example of a function which can flip an image, which is an affine transformation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def flip_axis(x, axis):\n",
    "    x = np.asarray(x).swapaxes(axis, 0)\n",
    "    x = x[::-1, ...]\n",
    "    x = x.swapaxes(0, axis)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, in the case where we are rotating, we define the rotation matrix as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set rotation angle, in radians.\n",
    "theta = 3.14\n",
    "\n",
    "rotation_matrix = np.array([\n",
    "    [np.cos(theta), -np.sin(theta), 0],\n",
    "    [np.sin(theta), np.cos(theta), 0],\n",
    "    [0, 0, 1]\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Feel free to play around with the `scipy.ndimage.affine_transform` command to try out other affine transformations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 6: Writing Output to Files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of Python's strengths is how well-suited it is to every piece of the image processing pipeline. Arguably, one of the most important parts of this pipeline is providing an output. In this section, we will briefly discuss a few ways to write this information to disk."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write one or more of the images you produced to a tif file using the `imsave` function from the `skimage.io` module. \n",
    "\n",
    "Make sure that the array you are writing is of integer type. If necessary, you can use the method `astype` for conversions - e.g. `some_array.astype(np.uint8)` or `some_array.astype(np.uint16)`. Be careful when converting a segmentation to uint8; if there are more than 255 cells, the 8bit format won't have sufficient bit-depth to represent all of the cell IDs.\n",
    "\n",
    "You can also try adding the segmentation to the original image, creating an image with two channels, one of them being the segmentation. After writing the file, load it into Fiji and check that everything worked as intended."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.io import imsave\n",
    "imsave(\"watershedseg.tif\", ws.astype(np.uint16))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NPZ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another uesful file type is a NumPy file. Numpy files allow fast storage and reloading arrays.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the function 'np.save' to save the array and reload it using 'np.load'.\n",
    "np.save(\"example_seg\", ws)  # Save\n",
    "seg = np.load(\"example_seg.npy\")  # Load\n",
    "print(ws.shape, seg.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### JSON"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, explore the documation and review another popular file type, [JSON](https://realpython.com/python-json/). Attempt to store and reload the same information from the NPZ example again using JSON."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "seg2 = pd.DataFrame(seg).to_json('data.json', orient='split')\n",
    "\n",
    "import json\n",
    "with open('example_seg2.txt', 'w') as outfile:\n",
    "    json.dump(seg2, outfile)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}