06-Image manipulation with Python.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Working with images in Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we seen briefly before, NumPy array can be used to represent image data, and therefore you can perform potentially arbitrarily complex image manipulation by manipulating NumPy arrays. \n",
    "\n",
    "Here we are going to take a deeper look at how image data is represented and handled in Python, and take a look at various basic image manipulation/processing operations that can be carried on. Let's get started! \n",
    "\n",
    "As always, the very first step should be importing of the essential packages: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Arbitrary 2D array as \"grayscale\" image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creating an \"image\" can be as simple as preparing any 2D NumPy array. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flat_image = np.ones((32, 32)) # creates 32x32 \"image\" filled with ones"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "can already be interpreted as an image. We plot the image content with `plt.imshow` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(flat_image)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When looking at an image, we would often want to remove the axis using `plt.axis('off')`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(flat_image)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This resulted in some flat image (we'll get to the question of the color very shortly)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Im order for an array to be interpreted as an image, it must be at least 2D in dimensions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gradation = np.linspace(0, 1, 5)\n",
    "\n",
    "# wouldn't work because it is just a 1D array\n",
    "plt.imshow(gradation)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can overcome this by giving the array an additional dimension with the `np.newaxis` trick:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = gradation[np.newaxis, :] # make 1 x 32 copy\n",
    "plt.imshow(x)\n",
    "x\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = gradation[:, np.newaxis] # make 32 x 1 copy\n",
    "plt.imshow(y)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A bit more interesting image can be used by doing a bit of broadcasting trick:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = x * y\n",
    "\n",
    "plt.imshow(image)\n",
    "plt.axis('off')\n",
    "\n",
    "image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Images to represent 2D function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just like a 1D array can be thought of representing a sequence of values read out from a function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(5)\n",
    "f = x**2\n",
    "\n",
    "plt.plot(f) # let matplotlib assume index as x-axis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can actually think of a 2D array to be representing readouts of a **2D function**. In the earlier example of the gradation image, the image could be thought of as representing a function:\n",
    "\n",
    "$z = x * y$\n",
    "\n",
    "such that the value of the $i^\\text{th}$ row and $j^\\text{th}$ column corresponds to the value of the function at $i^\\text{th}$ entry of the y-value and $j^\\text{th}$ entry of the x-value.\n",
    "\n",
    "To make things, concrete, let's take a look at another function:\n",
    "\n",
    "$z = (x + y) / 2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0, 1, 5)\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.linspace(0, 1, 5)[:, np.newaxis]\n",
    "y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = (x + y) / 2\n",
    "z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(z)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can compute arbitrary complex functions as you like. You can quickly get the appropriate x row vector and y column vector by using `np.ix_` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h, w = 11, 11\n",
    "cy, cx = h // 2, w // 2\n",
    "yv, xv = np.ix_(np.arange(h), np.arange(w))\n",
    "# z is the distance away from the center\n",
    "z = np.sqrt((yv - cy)**2 + (xv - cx)**2)\n",
    "\n",
    "\n",
    "plt.imshow(z)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### On the tangent: Plotting with interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When you think of a 2D array as representing your image, every grid location is essentially a **sample** from your 2D function. By default, the function will be plotted showing a **colored patch** centered at the sample, with the color corresponding to the function value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 5\n",
    "x = np.linspace(0, 1, N)\n",
    "y = np.linspace(0, 1, N)[:, np.newaxis]\n",
    "z = (x + y) / 2\n",
    "plt.imshow(z)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can achieve a **finer 2D function representation** by using more samples for both x and y axis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100\n",
    "x = np.linspace(0, 1, N)\n",
    "y = np.linspace(0, 1, N)[:, np.newaxis]\n",
    "z = (x + y) / 2\n",
    "plt.imshow(z)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, you can **estimate** finer mapping by **interpolating** between adjacent samples. To do so, pass in `interpolation` keyword argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 5\n",
    "x = np.linspace(0, 1, N)\n",
    "y = np.linspace(0, 1, N)[:, np.newaxis]\n",
    "z = (x + y) / 2\n",
    "plt.imshow(z, interpolation='bilinear')\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a more complex functions, the difference between actual finer grid of samples and an interpolated image becomes more apparent:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 10 # number of samples\n",
    "x = np.linspace(0, 1, N)\n",
    "y = np.linspace(0, 1, N)[:, np.newaxis]\n",
    "z = np.cos(np.sqrt((x-0.5)**2 + (y-0.5)**2) * 10 * np.pi)\n",
    "plt.imshow(z)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 10 # number of samples\n",
    "x = np.linspace(0, 1, N)\n",
    "y = np.linspace(0, 1, N)[:, np.newaxis]\n",
    "z = np.cos(np.sqrt((x-0.5)**2 + (y-0.5)**2) * 10 * np.pi)\n",
    "plt.imshow(z, interpolation='bilinear')\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100 # number of samples\n",
    "x = np.linspace(0, 1, N)\n",
    "y = np.linspace(0, 1, N)[:, np.newaxis]\n",
    "z = np.cos(np.sqrt((x-0.5)**2 + (y-0.5)**2) * 10 * np.pi)\n",
    "plt.imshow(z)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With some interesting values in the image, it's time to talk about what exactly the color means in our image plot."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Colormap for a Grayscale Image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When working with a 2D image, each pixel location (row, col location in the 2D array) contains a single value. Such image is typically referred to as a **grayscale**. When plotting a grayscale image, Matplotlib makes use of a **colormap** to decide exactly what color will be used to represent a particular value in your 2D array.\n",
    "\n",
    "You can get the the mapping displayed by adding a **color bar** to your figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = x + y\n",
    "\n",
    "plt.imshow(image)\n",
    "plt.axis('off')\n",
    "# add a colorbar to show the mapping\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default Matplotlib choses some colormap, here mapping values from 0 (smallest value found in the image) to 1 (largest value in the image) into a gradation of colors from blue to yellow, respectively. You can specificy what colormap should be used by passing in `cmap` keywrod argument to `imshow`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = x + y\n",
    "\n",
    "# use cmap argument to change the color map used\n",
    "plt.imshow(image, cmap='gray')\n",
    "plt.axis('off')\n",
    "\n",
    "# add a colorbar to show the mapping\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using \"gray\" colormap maps different shades of gray to the values of your image, thereby giving rise to the name \"grayscale\" image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aside from specifying some colormap with shorthand names like \"gray\" or \"cool\"..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = x + y\n",
    "\n",
    "# use \"cool\" cmap\n",
    "plt.imshow(image, cmap='cool')\n",
    "plt.axis('off')\n",
    "\n",
    "# add a colorbar to show the mapping\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can find listing of colormaps under `plt.cm`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = x + y\n",
    "\n",
    "# use colormap plt.cm.gist_earth\n",
    "plt.imshow(image, cmap=plt.cm.gist_earth)\n",
    "plt.axis('off')\n",
    "\n",
    "# add a colorbar to show the mapping\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And of course, you can find a good listing of all colormaps provided out-of-box in Matplotlib by [checking out their documentation](https://matplotlib.org/tutorials/colors/colormaps.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Automatic scaling of colormap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Noticed that the colormap was mapped to values between 0 and 1. What happens if our image contains values outside of this range? For instance, let's adjust how we compute the image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-1, 1, 5)\n",
    "y = np.linspace(-1, 1, 5)[:, np.newaxis]\n",
    "\n",
    "z = x + y\n",
    "\n",
    "z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that our image ranges in values from -2 to 2, let's see what happens if we plot the image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(z, cmap='gray')\n",
    "plt.axis('off')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the the whole range of color (back to white in grayscale) is mapped to values between -2 and 2!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, when you pass in an 2D image to Matplotlib to plot, it will map the entire range between the minimum and the maximum value of the image to the color map!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adjusting the mapping range"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although this **auto-scaling** of image value range can be useful, sometimes, you would want to control exactly what range of values gets mapped to the colormap, giving you some consistency across different images. For example, if I prepare scaled versions of the images:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-1, 1, 5)\n",
    "y = np.linspace(-1, 1, 5)[:, np.newaxis]\n",
    "\n",
    "image1 = x + y\n",
    "\n",
    "image2 = image1 * 0.5 # scaled to half of image1\n",
    "\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.imshow(image1)\n",
    "plt.axis('off')\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.imshow(image2)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because of autoscaling, you cannot tell that these two images actually have differnet values. Plotting the colorbars can help a bit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-1, 1, 5)\n",
    "y = np.linspace(-1, 1, 5)[:, np.newaxis]\n",
    "\n",
    "image1 = x + y\n",
    "\n",
    "image2 = image1 * 0.5 # scaled to half of image1\n",
    "\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.imshow(image1)\n",
    "plt.axis('off')\n",
    "plt.colorbar()\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.imshow(image2)\n",
    "plt.axis('off')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "...but the difference is not intuitive! You can actually control the range of color mapping with `vmin` and `vmax` keyword arguments to `imshow`!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-1, 1, 5)\n",
    "y = np.linspace(-1, 1, 5)[:, np.newaxis]\n",
    "\n",
    "image1 = x + y\n",
    "\n",
    "image2 = image1 * 0.5 # scaled to half of image1\n",
    "\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.imshow(image1, vmin=-2, vmax=2)\n",
    "plt.axis('off')\n",
    "plt.colorbar()\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.imshow(image2, vmin=-2, vmax=2)\n",
    "plt.axis('off')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the two images **are plotted with the same colorbar**, and the difference between the two images (datawise) is much more visually apparent!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have discussed quite a bit about 2D array as a gray scale image, let's move onto talking about **colored images**!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Working with colored images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We already saw how a 2D array can be interpreted as a \"grayscale\" image with the exact pixel values mapping into different color controlled by colormap. But what if I want to represent a **colored image**? This is where a 3D array as image comes into! To get a sense of this, let's go ahead and take a look at a real colored image:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading external images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although you could create an arbitrary NumPy arrays from scratch and treat them as images as we have been, you are far more likely to want to get images loaded from external sources - i.e. an image file!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ability to load in common image file formats is provided by a package called `imageio` (along with many other packages)!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import imageio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = imageio.imread('data/image_japan.jpg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a quick look at the image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(image)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now take a closer look at this **colored** image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3D arrays as images with color channels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you take a look at the shape of the loaded color image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "you should notice two things:\n",
    "\n",
    "* the array is now 3D - there is an extra trailing dimension\n",
    "* and the shape of the array in the third dimension is 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This has to do a lot with how colors are presented in computers: as combination of red, green and blue values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Colored images are thought of being composed of 3 distinct **channels**, each specifying the red, greeen and blue values for each pixel. You can actually see this by taking a colored image, and turning all but one channel 0:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "red_image = imageio.imread('data/image_japan.jpg')\n",
    "\n",
    "# turn Green and Blue channels 0\n",
    "red_image[:, :, 1] = 0\n",
    "red_image[:, :, 2] = 0\n",
    "\n",
    "plt.imshow(red_image)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "green_image = imageio.imread('data/image_japan.jpg')\n",
    "\n",
    "# turn Red and Blue channels 0\n",
    "green_image[:, :, 0] = 0\n",
    "green_image[:, :, 2] = 0\n",
    "\n",
    "plt.imshow(green_image)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "blue_image = imageio.imread('data/image_japan.jpg')\n",
    "\n",
    "# turn Red and Green channels 0\n",
    "blue_image[:, :, 0] = 0\n",
    "blue_image[:, :, 1] = 0\n",
    "\n",
    "plt.imshow(blue_image)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A full colored image results from combining these individual channels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(red_image + green_image + blue_image)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating your own colored image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you know the secret behind a colored image, you can imagine creating some fancy colored image by simply stacking up three separate images, corresponding to red, green and blue channel!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100 # number of samples\n",
    "x = np.linspace(0, 1, N)\n",
    "y = np.linspace(0, 1, N)[:, np.newaxis]\n",
    "\n",
    "red = x + 0 * y\n",
    "green = 0 * x + y\n",
    "blue = np.sqrt((x - 0.5)**2 + (y-0.5)**2)\n",
    "\n",
    "plt.subplot(1, 3, 1)\n",
    "plt.imshow(red, cmap='gray')\n",
    "plt.axis('off')\n",
    "plt.title('Red channel')\n",
    "\n",
    "plt.subplot(1, 3, 2)\n",
    "plt.imshow(green, cmap='gray')\n",
    "plt.axis('off')\n",
    "plt.title('Green channel')\n",
    "\n",
    "plt.subplot(1, 3, 3)\n",
    "plt.imshow(blue, cmap='gray')\n",
    "plt.axis('off')\n",
    "plt.title('Blue channel')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "combined = np.stack([red, green, blue], -1)\n",
    "plt.imshow(combined)\n",
    "plt.axis('off')\n",
    "plt.title('Combined channels')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### On tangent: more structured plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Did you notice that you were now repeating something very similar over and over? It's time to structure your plotting code to get it done more quickly!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100 # number of samples\n",
    "x = np.linspace(0, 1, N)\n",
    "y = np.linspace(0, 1, N)[:, np.newaxis]\n",
    "\n",
    "red = x + 0 * y\n",
    "green = 0 * x + y\n",
    "blue = np.sqrt((x - 0.5)**2 + (y-0.5)**2)\n",
    "\n",
    "channels = [red, green, blue]\n",
    "labels = ['Red', 'Green', 'Blue']\n",
    "\n",
    "fig, axs = plt.subplots(1, 3, dpi=300)\n",
    "\n",
    "for ax, channel, label in zip(axs, channels, labels):\n",
    "    ax.imshow(channel, cmap='gray')\n",
    "    ax.axis('off')\n",
    "    ax.set_title('{} channel'.format(label))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The new guy here is `zip` function: `zip` takes in one or more **iterables** and **zips together** same positioned elements from each and every iterable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "names = ['Maku', 'Chabo', 'Lulu']\n",
    "colors = ['Black', 'Brown', 'Gray']\n",
    "ages = [7, 6, 1]\n",
    "\n",
    "for name, color, age in zip(names, colors, ages):\n",
    "    print('{} is a {} dog and {} years old!'.format(name, color, age))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From time to time, using `zip` with Matplotlib's object based plotting (e.g. getting an array of axes) can save a lot of typing!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Image data types"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We saw that when working with a grayscale image (that is, a 2D image), then you can control how the value of a pixel maps to a color using a *color map*, and you can even control the range using `vmin` and `vmax`. How does it work for a colored image?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It turns out that Matplotlib understands two types of colored images:\n",
    "\n",
    "1. Integer images with values between 0 and 255 for each pixel and channel.\n",
    "2. Float images with values between 0.0 and 1.0 for each pixel and channel."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have already encountered both types. You tend to get an integer image when you load images from a file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = imageio.imread('data/image_japan.jpg')\n",
    "\n",
    "plt.imshow(image)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image.max(), image.min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image.dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the data type (`dtype`) of `uint8` means 8-bit unsigned integer. In other words, it can hold integer values between 0-255."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we constructed an image by stacking channels, we constructd a **float** value NumPy array that actually ranged between 0 and 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100 # number of samples\n",
    "x = np.linspace(0, 1, N)\n",
    "y = np.linspace(0, 1, N)[:, np.newaxis]\n",
    "\n",
    "red = x + 0 * y\n",
    "green = 0 * x + y\n",
    "blue = np.sqrt((x - 0.5)**2 + (y-0.5)**2)\n",
    "combined = np.stack([red, green, blue], -1)\n",
    "\n",
    "plt.imshow(combined)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "combined.min(), combined.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "combined.dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happens if we would have scaled the floating image outside of range [0, 1]?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scaled_image = combined * 2 # double the value\n",
    "scaled_image.min(), scaled_image.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(scaled_image)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You actually get a warning that the input value gets **clipped** between 0 and 1 for floating image or 0 and 255 for integer image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Image manipulation and processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have learned quite a bit on how to load an image as well as how an image is represented using NumPy array, it's time to perform some basic **manipulation/processing** on image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's go ahead an work with an image that you can find from `scipy.misc`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.misc import face\n",
    "image = face(gray=True) # get gray scale image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(image, cmap='gray')\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Manipulating image by NumPy array manipulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since an image is nothing but a NumPy array, you can perform any valid manipulation on the array, resulting in a manipulation of the image represented by the array."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cropping an image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can **crop** an image by simply performing slicing on the first two dimensions of the image, corresponding to the *height* and *width*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(image[200:-200, 200:-200], cmap='gray')\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Flipping an image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can use reversed slicing (`::-1`) to easily **flip** an image:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Flip horizontally (left-right flip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(image[:, ::-1], cmap='gray')\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Flip vertically (up-down flip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(image[::-1, :], cmap='gray')\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fill in values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Manipulation of array content, results in modification of the image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = face(gray=True) # get gray scale image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image_c = image.copy() # make a copy of image\n",
    "h, w = image_c.shape # get height and width\n",
    "image_c[:, w//2:] = 127 # replace the right half with 127\n",
    "plt.imshow(image_c, cmap='gray', vmin=0, vmax=255)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Making of an image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can achieve more interesting results by combining two or more images/arrays together."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First let's create an image of same size and fill it with 0s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = np.zeros(image.shape)\n",
    "\n",
    "plt.imshow(mask, cmap='gray')\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we multiply this mask with an image, then obviously everything becomes black (0)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "masked_image = image * mask\n",
    "plt.imshow(masked_image, cmap='gray', vmin=0, vmax=255)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead, let's change the mask, but filling some center rectangle with 1's."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vpad = 100\n",
    "hpad = 150\n",
    "\n",
    "mask[vpad:-vpad, hpad:-hpad] = 1\n",
    "\n",
    "plt.imshow(mask, cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use this as a **mask** on the image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "masked_image = image * mask\n",
    "plt.imshow(masked_image, cmap='gray')\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### More complex masks with boolean operations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using arrays of x and y values and NumPy boolean arrays, you can create fancier masks."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Circular mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h, w = image.shape\n",
    "yv, xv = np.ix_(np.arange(h), np.arange(w))\n",
    "yc, xc = h // 2, w // 2\n",
    "\n",
    "# compute distance from the center\n",
    "r = np.sqrt((xv - xc)**2 + (yv - yc)**2)\n",
    "\n",
    "plt.imshow(r)\n",
    "plt.axis('off')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above, each pixel value is the distance from the center. We can now only select pixels with distance less than say 300:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = r < 300"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(mask, cmap='gray')\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we can combine this with our image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "masked_image = mask * image\n",
    "\n",
    "plt.imshow(masked_image, cmap='gray')\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Thresholding an image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In fact use of boolean array can be a quite powerful technique in image manipulation. A very common technique is **thresholding** where you select only pixels that is either above or below a specific value, or threshold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thresholded_image = image > 170\n",
    "\n",
    "plt.imshow(thresholded_image)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that `thesholded_image` is a boolean array. We can now use this as a mask onto the original image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "masked_image = image * thresholded_image\n",
    "\n",
    "plt.imshow(masked_image, cmap='gray')\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can of course use the boolean mask in *fancy indexing*, changing values of only `True`'s in the mask:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image_c = image.copy()\n",
    "\n",
    "image_c[thresholded_image] = 0 # make bright parts black\n",
    "\n",
    "plt.imshow(image_c, cmap='gray', vmin=0, vmax=255)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filtering images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of a very common yet extremely powerful image processing technique is **filtering** whereby you apply a filter to every pixel of an image. Depending on the filter you use, you can achieve various interesting image manipulations such as blurring or edge enhancement."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although we can use functions like `np.convolve` to achieve the filtering ourselves, various interesting filtering operations can be found inside `scipy.ndimage` subpackage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import  ndimage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gaussian blur"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can **blur** or **smoothen out** your image by letting each pixel be replaced by an average of nearby pixels. A very common filter for this job is **Gaussian filter**, which performs weighted average of pixel values weighting pixels based on Gaussian distance profile - more weights are given to nearby pixels.\n",
    "\n",
    "You control the size of the Gaussian filter (how far out it looks out to average) by controling the standard deviation or **sigma** of the filter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = face(gray=True) # load image again\n",
    "\n",
    "blurred_image = ndimage.gaussian_filter(image, sigma=3)\n",
    "very_blurred = ndimage.gaussian_filter(image, sigma=10)\n",
    "\n",
    "images = [image, blurred_image, very_blurred]\n",
    "sigmas = [0, 3, 10]\n",
    "\n",
    "fig, axs = plt.subplots(1, len(images), dpi=300)\n",
    "\n",
    "for ax, image, sigma in zip(axs, images, sigmas):\n",
    "    ax.imshow(image, cmap='gray')\n",
    "    ax.axis('off')\n",
    "    ax.set_title('sigma={}'.format(sigma))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sharpening an image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can achieve a bit of an opposite of blurring - that is, image sharpening - true a clever use of a **sharpening filter**. Alternatively, it turns out that the same result can be subtracting off Gaussian blurred image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = face(gray=True).astype(float)\n",
    "blurred_image = ndimage.gaussian_filter(image, 3)\n",
    "filter_blurred_image = ndimage.gaussian_filter(blurred_image, 1)\n",
    "alpha = 30\n",
    "sharpened_image = blurred_image + alpha * (blurred_image - filter_blurred_image)\n",
    "\n",
    "images = [image, blurred_image, sharpened_image]\n",
    "titles = ['Original', 'Blurred', 'Re-sharpened']\n",
    "\n",
    "fig, axs = plt.subplots(1, len(images), dpi=300)\n",
    "\n",
    "for ax, image, title in zip(axs, images, titles):\n",
    "    ax.imshow(image, cmap='gray')\n",
    "    ax.axis('off')\n",
    "    ax.set_title(title)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Foreground-Background separation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When processing an image, you often encounter a need to separate the background from the foreground. This may be an essential step in performing **object detection** whereby you identify areas in the image that corresponds to different objects (e.g. cells). Consider the following image of some bacterial blobs or colonies:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = plt.imread('data/blobs.jpg')\n",
    "plt.imshow(image)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This image is loaded actually as RGB image. To make it easier to work with lets **flatten out** the color by simply average out the last dimension:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = image.mean(axis=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(image)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could try to detect the foreground by theresholding the image:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To come up with a suitable threshold, let's look at the distribution of the image intensity:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(image.ravel(), bins=60)\n",
    "plt.xlabel('Image intensity')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There appears to be nice separation between the two humps - likely, one is for the background and the other for the foreground. Let's take this to set the value of the threshold:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = image > 125\n",
    "\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.imshow(image)\n",
    "plt.axis('off')\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.imshow(mask)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally you can use `ndimage`'s `label` function to label each contiguous areas: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "label_im, nb_labels = ndimage.label(mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(label_im)"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}