Image_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Image analysis notebook\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Image analysis is an important part of materials science. Microstructures, or morphologies are often imaged using microscopes and often we are interested in quantifying, that is, measuring different features. These could be pore areas, cell sizes, fibre lengths, grain sizes, crack widths etc. In some cases we might be interested in analysing the image of a spectrum or a diffraction pattern for example. Here is a gratuitously pretty image of the microstruture of a 3D printed aerospace titanium alloy:\n",
    "\n",
    "![3D printed Ti Alloy Microstructure](http://www.cl.eps.manchester.ac.uk/medialand/chem_eng/images/research/facilities-institutes/emc-1_340x140.jpg)\n",
    "\n",
    "Since it is now easy to digitise images, it is easy to exploit the power of computers to process and analyse them. Digital image processing is a big subject and there are many [multi-volume books](http://www.springer.com/us/book/9781848001909) dedicated to image processing techniques. For general image processing, it's often easier to use a dedicated image processing program like the excellent [imageJ](http://imagej.net/Welcome), for example. \n",
    "\n",
    "In this notebook, we will use python to do some simple image processing and analysis. The aim is to understand how digital images are structured and to explore how the power of python and its packages can be used to make scientific measurements on images.\n",
    "\n",
    "Special Thanks of this notebook goes for Dr. James Gibson."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As usual, we start by importing the usual useful python modules <code>numpy</code> and <code>matplotlib</code>:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#impot libraries \n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this notebook we also suppress warnings. This doesn't interfere with the code and avoids confusing warning (not error) messages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ignore warnings \n",
    "import warnings\n",
    "warnings.simplefilter('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Images as arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step in working with digital images is to get used to the idea that images are just arrays of numbers. For example, the next array might not mean much to you: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mysterious =np.array(\n",
    "[[125,118,115,117,122,118,122,131,130,121,122,130,126,121,123,120,129,128],\n",
    " [140,134,130,128,132,140,137,132,136,137,136,140,133,128,130,130,132,133],\n",
    " [143,147,141,136,146,117, 57, 24, 40, 69,117,146,136,130,136,138,139,144],\n",
    " [156,153,145,159,128, 59,105, 64, 13,  0,  0,129,153,139,144,151,152,145],\n",
    " [164,159,159,162, 91,179,198,163, 75, 17,  0, 57,157,158,159,158,159,156],\n",
    " [175,165,173,128,102,193,196,177,118, 61, 19,  0,115,141,152,156,162,165],\n",
    " [158,181,182, 72, 97,161,170,131,115, 68, 10,  0, 56, 92,117,131,160,167],\n",
    " [100,158,146, 51, 83,117,126,101,101, 89,  0,  9, 40, 67, 81, 81,155,168],\n",
    " [ 65, 70, 66,  0,133,178,143,166,168, 84,  0, 13, 32, 74, 83, 75,132,147],\n",
    " [ 54, 49, 37,  0, 96,162,108,145,123, 33, 14, 12, 19, 50, 68, 73, 88, 91],\n",
    " [ 53, 48, 38, 17, 45,137, 96,101, 97, 36, 12, 11, 18, 37, 40, 34, 45, 60],\n",
    " [ 44, 51, 54, 37,  0,105,146, 88, 62, 26, 15,  9, 18, 45, 53, 49, 54, 54],\n",
    " [ 60, 77, 78, 33, 16,  0, 70, 33, 12, 28, 19, 11, 13, 36, 50, 47, 42, 33],\n",
    " [ 77, 89, 79, 33, 27,  0, 53, 87, 56, 75, 41,  6, 18, 32, 58, 61, 60, 72],\n",
    " [ 50, 42, 38, 23,  4, 29,113,134,125,131, 77, 20, 31, 12, 34, 54, 78, 97],\n",
    " [ 38, 20, 25,  0, 47,148,182,174,175,161, 90, 41, 29,  5, 13, 87,111, 86],\n",
    " [ 50, 45, 26,  0,113,193,197,187,188,162, 88, 31, 27, 56, 33, 38, 82, 98],\n",
    " [ 77, 33,  0, 93,181,196,198,199,193,167, 88, 68, 94, 83, 46,  1, 16, 52],])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now use <code> matplotlib </code> to plot it. We are going to use the `imshow` function, which is similar to the `matshow` function we used in the first week to plot matrices. Do you recognise the image?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# show the array\n",
    "plt.imshow(mysterious,cmap='gray',interpolation='none');\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\"> \n",
    "Notice how, in the `imshow` command above, we *pass on* two arguments: `cmap` and `interpolation`. These are options in `imshow` that tell the function to use a grayscale to represent the values of the array and to not interpolate the image.\n",
    "\n",
    "Usually, you can find out what kind of arguments a function uses by calling up the documentation string of the function, with the question mark magic:\n",
    "<div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plt.imshow?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most functions will have a signature that lists all the possible arguments. As you can see `imshow` has many arguments, and they are all optional. You don't need to know all these options, not many people do! You will discover them as you use different functions more and more, and what you want to achieve changes.\n",
    "<div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Back to image analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this image each pixel corresponds to one value, ranging from 0 to 256. We can plot a histogram that shows the distribution of values making up the image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Histogram \n",
    "\n",
    "plt.hist(mysterious.flatten()); # note that we flatten the array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we use `.flatten` to flatten the array into a 1D array for plotting the histogram. Otherwise the `hist` function will try to plot a histogram for each row. Try it!\n",
    "\n",
    "The point is that an image is just an array of intensity values and can be plotted using any \"colour scale\" by using different colour maps. For example with <code>jet</code> it looks like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# colour maps\n",
    "\n",
    "plt.imshow(mysterious,cmap='jet',interpolation='none');\n",
    "plt.colorbar(); # colorbar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time we have added a `colorbar` to the plot. This shows how the colours relate to the intensity values in the arrays. \n",
    "As you can see the choice of color map can significantly change your image. Mona Lisa in <code>jet</code> doesn't really look like Mona Lisa. Let's try another color map, <code>bone</code>:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# colour bar\n",
    "plt.imshow(mysterious,cmap='bone',interpolation='none');\n",
    "plt.colorbar();\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can list all available colourmaps by running the `colormaps` command: \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plt.colormaps()\n",
    "plt.colormaps()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are rather a lot. You might prefer to browse them [here](http://matplotlib.org/examples/color/colormaps_reference.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: red\"> Task 1:</span> Display the `mysterious` array using a different colormap of your choice:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: blue\"> Solution:</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Image formats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although all images are all ultimately arrays of pixels, image files can store this information in a number of ways, or formats. In some cases all the values are stored without loss (e.g. .bmp and .tif) but these files can quickly become unfeasibly large. Instead, files are usually stored in \"lossy\" compressed formats (e.g. .jpg and .png). These formats contain the intensity values in encoded, compressed format and the image needs to be decoded using an algorithm. Opening such files can therefore be difficult, but because we are using Python ...wait for it... *there is a package for that.*\n",
    "\n",
    "\n",
    "We are going to use the scikit-image package, called [skimage](http://scikit-image.org/docs/dev/api/skimage.html). To open image files, we need the input/output <code>io</code> subpackage:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import libraries\n",
    "from skimage import io\n",
    "from skimage import color\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, opening an image in any standard format is only one command away: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read th image from the folder\n",
    "image=io.imread('images/Grey_Mona_lisa.jpg')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then plot it using <code>matplotlib</code>:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# show the image via matplotlib\n",
    "plt.imshow(image);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rgb2gray\n",
    "\n",
    "gs_image=color.rgb2gray(image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dimensions of the image\n",
    "\n",
    "np.shape(gs_image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot\n",
    "plt.imshow(gs_image,'plasma')\n",
    "plt.axis('off')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use the `imshow` function from `skimage` to plot the image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# plt vs. io\n",
    "io.imshow(image);\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function `imshow` from `skimage` does not rescale the intensity levels. What does this mean? Whereas `imshow` from `matplotlib` will give you the best view of the data by using the full colourscale and setting the minimum value to black and he maximum to white, `imshow` from `skimage` shows you the most faithful representation of the image by setting 0 to black and 255 (the highest intensity for an 8-bit image: $2^8=256$) to white. This way the original look of the image is preserved. In this case they are almost the same but dark and bright images would look very different."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: red\"> Task 2:</span> Crop the large Mona Lisa image above just below the neck and display it using <code>io.imshow</code>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Tip:** use the array slicing you learned in the \"Introduction to Python\" workbook to select the columns and rows you want."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: blue\"> Solution:</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Analysing images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Now we'll look at some images that might be more interesting to materials scientists.\n",
    "\n",
    "Many of the PhD students at the IMM have been investigating creep-resistant magnesium alloys. Indeed, Prof. Korte-Kerzel's funded SFB (SFB 1394) will also focus on this material system.\n",
    "\n",
    "Here is an SE image of the microsctructure of such a creep-resistant alloy (Mg-5Al-3Ca). Read the image file in with the `io.imread` function from `skimage` and display it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Load the image and show\n",
    "\n",
    "alloy=io.imread('images/MgAlCa.png')\n",
    "io.imshow(alloy);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The creep-resistance in this alloy comes from the intermetallic skeleton, which slows down diffusion and help prevents grain boundary sliding. However, this intermetallic is macroscopically brittle, so our research focusses on understanding its intrinsic mechanical properties, as well as how plastic deformation is transferred between the soft but ductile magnesium grains to the hard but brittle intermetallic Laves phases.\n",
    "\n",
    "These Laves phases are the light-coloured part of the image. For simplicity, we'll pretend there is only one type of Laves phase present, giving us a two-phase microsctructure.\n",
    "\n",
    "The strength and creep-resistance of this alloy depends in part on how much Laves phase is actually present, so we would like some way to quantify this.\n",
    "\n",
    "We can zoom into any region of the image to take a closer look. First it's useful to see how many rows and columns our image has. To do this we use the function `shape` from `numpy`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# dim.s of image\n",
    "np.shape(alloy)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The function tells us the image is 743 pixels high, 899 pixels wide and 3 high."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "#### <span style=\"color:red\">Question: What is the \"3\" referring to?</style>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: blue\"> Solution:</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This third dimension will cause some problems later on when we do the thresholding tasks. As there is no colour information that we care about at the moment, let's import an image saved as a grayscale image, with no colour information: \"MgAlCa_gray.png\".\n",
    "\n",
    "Read the image file in with the `io.imread` function from `skimage`, display it, and check out the shape:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# grey MgAlCa\n",
    "\n",
    "alloy=io.imread('images/MgAlCa_gray.png')\n",
    "io.imshow(alloy);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see the image looks identical to the one above (this makes sense as it's still gray). But, importantly our RGB values have been replaced by a simple intensity value, leaving us with a 2D array, 743 pixels high by 899 pixels wide.\n",
    "\n",
    "Note we have also overwritten our inital `alloy` variable with this new, gray one."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Let's now display the top-left 50$\\times$50 pixels, using the array slicing commands we learnt in the 'Introduction to Python' workbook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Zoom in\n",
    "\n",
    "io.imshow(alloy[:50,:50]);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### <span style=\"color: red\"> Task 3:</span>  Use array slicing to plot different regions of the image and identify a region where the two phases are clear and representative."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### <span style=\"color: blue\"> Solution:</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "As you can see the Laves phases are mostly white or light grey but have darker borders and are embedded in a matrix which has some other shades of grey.\n",
    "\n",
    "If we could find a way of counting just the light pixels, that would give us the area fraction of Laves phase, wouldn't it?\n",
    "\n",
    "Lets use the function <code>histogram</code> from the <code>exposure</code> subpackage from <code>skimage</code> to plot the distribution of grey scale values. First, we import the subpackage exposure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import exposure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can plot the histogram. We are using a new `matplolib` function called `fill_between`. Do you know what it does? How can you find out?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# histogram of exposure\n",
    "histogram=exposure.histogram(alloy)\n",
    "plt.fill_between(histogram[1],histogram[0],alpha=0.5, color='b')\n",
    "plt.xlim([0,256]);\n",
    "plt.ylabel('Number of pixels')\n",
    "plt.xlabel('Intensity (greyscale value)')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The histogram shows the distribution of grey scale values, with 0 corresponding to black and 255 to white. As you can see there are no black pixels and there is a spiked distribution around an intensity of 50. There is also a small but very sharp peak at 256. Can you explain where it comes from?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Separating the Laves phase from the matrix - thresholding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "In order to measure the amount of Laves phase, we need to first be able to separate it from the matrix. We know that it is generally light in colour, therefore we need to find a way of separating these pixels from the rest.\n",
    "\n",
    "We need to create a boolean (true/false) array, which is 1 when the pixels are light (intensity > (a certain value)) and 0 when the pixels are not. Using <code>numpy</code>, this is easy. We assign a new array, called `light_pixels` to the result of the comparison between our image and the value 70, which we'll use as a first guess."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# thresh\n",
    "light_pixels = alloy > 70\n",
    "light_pixels\n",
    "np.shape(light_pixels)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "This process is called _**thresholding**_. Lets display the resultant boolean array. This is not an image, but an array of true false. We can use `imshow` from `matplotlib` to display it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# show\n",
    "plt.imshow(light_pixels,cmap='gray')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "To see if we've been successful in separating the particles, we should compare the image of the boolean array with the original image, side by side using <code>subplots</code> from <code>matplotlib</code>:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# apply thresh\n",
    "f,images=plt.subplots(1,2, figsize=(16,9))# the figsize option allows us to change the size of the image\n",
    "images[0].imshow(alloy,cmap='gray',interpolation='none' )\n",
    "images[1].imshow(light_pixels, cmap='gray');\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<div class=\"alert alert-success\">\n",
    "In the cell above, we have used an advanced form of plotting that enables us to display more than one set of plots in one figure. First we create the subplots using `plt.subplots`, which specifies the number of rows and columns we need (1,2 in this case). We also use the `figsize` argument to change the size of the figure (16,9 in this case for widescreen!). When we call the `subplots` function we assign two new variables: `f`, which is a handle to the figure, and `images` which is a list of all the subplots in the figure. To plot in each of the subplots, we run the plotting functions on each subplot. So `image[0].imshow(array)` will display `array` in the first subplot, `image[1].imshow(array)` in the second etc.\n",
    "<div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### Back to image analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "It's easier to compare the images if we zoom in:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Zoom in\n",
    "f,images=plt.subplots(1,2, figsize=(16,9))\n",
    "images[0].imshow(alloy[-250:,-650:-250],cmap='gray',interpolation='none')\n",
    "images[1].imshow(light_pixels[-250:,-650:-250],cmap='gray',interpolation='none');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "It appears that although we have managed to separate some of the Laves phase, we have not picked up all of it. This means our threshold value of 70 is too high. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### <span style=\"color: red\"> Task 4:</span> Use different values of threshold to separate the Laves phase and choose the one you think best represents the microstructure in the original image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### <span style=\"color: blue\"> Solution:</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Automatic thresholding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "As you can see it's not easy to determine the \"ideal\" threshold and different people will choose different values. This makes thresholding problematic because it depends on the user, i.e. it is subjective.\n",
    "\n",
    "To avoid this, we can use an automatic thresholding algorithm, like [Otsu's method](https://en.wikipedia.org/wiki/Otsu%27s_method), which we can also find in <code>skimage</code>:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Automatic thresh lib\n",
    "from skimage.filters import threshold_otsu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Let's see what Otsu's method gives us:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# print thresh\n",
    "thresh=threshold_otsu(alloy)\n",
    "print(thresh)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Let's see what this threshold value gives now in terms of the images:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# subplots\n",
    "otsu_pixels = alloy > thresh\n",
    "f,images=plt.subplots(1,2, figsize=(16,9))\n",
    "images[0].imshow(alloy,cmap='gray',interpolation='none' )\n",
    "images[1].imshow(otsu_pixels, cmap='gray');\n",
    "f,images=plt.subplots(1,2, figsize=(16,9))\n",
    "images[0].imshow(alloy[-250:,-650:-250],cmap='gray',interpolation='none')\n",
    "images[1].imshow(otsu_pixels[-250:,-650:-250],cmap='gray',interpolation='none');\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How well do you think the Otsu thresholding does?\n",
    "\n",
    "#### <span style=\"color: red\"> Question:</span> Even if this is \"worse\" than the value you found, why might we want to use it anyway?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "#### <span style=\"color: blue\"> Answer:</span> \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Measuring area fraction "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have implemented thresholding, we can measure the area fraction of Laves phase in our alloy. We do this simply by dividing the number of white pixels in the thresholded, binary image by the total number of pixels. There are many different ways of doing this but since the white pixels are 1 and the black pixels are 0, we can just sum the binary image array and divide by the total number of pixels.\n",
    "\n",
    "So for a manual threshold value of 50 the fraction of particles is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_thresh = alloy > 50\n",
    "total_pixels=np.size(my_thresh) # total pixels in the image\n",
    "white_pixels=my_thresh.astype(int).sum() #sums all the values in my_tresh after converting it into an integer array.\n",
    "my_particle_frac=white_pixels/total_pixels #get pixel fraction\n",
    "print('My particle fraction is {0:.2f}%'.format(my_particle_frac*100)) #print it nicely"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Whereas for Otsu's method it is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_thresh = alloy > thresh # remember thresh is the value the automatic method (Otsu's method) found.\n",
    "total_pixels=np.size(my_thresh)\n",
    "white_pixels=my_thresh.astype(int).sum()\n",
    "my_particle_frac=white_pixels/total_pixels\n",
    "print('Otsu\\'s particle fraction is {0:.2f}%!'.format(my_particle_frac*100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "\n",
    "There are a few of new things happening in the cells above. First we are using the `np.size` function, which tells us the number of elements in an array. \n",
    "Then we have a chaining of methods on an array:\n",
    "\n",
    "...`= my_thresh.astype(int).sum()`\n",
    "\n",
    "what is happening here? Well it's easier then it seems. First we are using the `astype(integer)` method on the my_tresh array, which is boolean (True/False), to convert it into an array of integers containing 0s and 1s. Then we apply the method `sum`, which returns the sum of all the elements of the array. We do both these things because it is easier to understand what summing ones and zeros is, than what summing 'True' and 'False' is. This chaining of methods avoids defining an intermediate variable. For example, we could have done:\n",
    "\n",
    "`total_pixels=np.size(my_thresh)`\n",
    "\n",
    "`white_pixels_integers=my_thresh.astype(int)`\n",
    "\n",
    "`white_pixels_sum=white_pixels_integers.sum()`\n",
    "\n",
    "would achieve the same result.\n",
    "<div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The thresholding examples above show how, for an image with a very smooth distribution of greyscale values, getting an exact value of volume fraction is very difficult. With python we can quickly see how the calculated volume fraction depends on the selected threshold. To do this we first create an array of threshold values using `linspace`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thresh_values=np.linspace(30,256,226)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we iterate over each value and calculate the fraction of particles, using the same code we used earlier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "frac=[]\n",
    "for val in thresh_values:\n",
    "    my_thresh = alloy > val\n",
    "    total_pixels=np.size(my_thresh)\n",
    "    white_pixels=my_thresh.astype(int).sum()\n",
    "    my_particle_frac=white_pixels/total_pixels\n",
    "    frac.append(my_particle_frac*100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we plot the fraction of particles calculated using `matplotlib`. Notice the use of a new `matplotlib` function, `axvline`. What does it do?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(thresh_values,frac,'o')\n",
    "plt.axvline(60,color='red')\n",
    "plt.ylabel('Particle fraction (%)')\n",
    "plt.xlabel('Threshold');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Otsu's method picks a value where the change in particle fraction changes slope. Is that the \"right\" value? It's impossible to be sure. What this shows is that if we want to measure phases, particles or other interesting features of the microstructure, it is helpful to separate them during the image acquisition using the right imaging conditions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Particle size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To finish this session, we will look at another important parameter in materials science: the particle size and shape. To measure this we need to be able to separate individual particles. For this task we will use a different image.\n",
    "\n",
    "The image below is a high resolution image of Ni nanoparticles, imaged with a high resolution transmission electron microscope (TEM). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nanoparticles=io.imread('images/nanoparticles.tif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "io.imshow(nanoparticles)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step is to use thresholding to separate the images from the background. We'll use Otsu's method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nano_thresh = threshold_otsu(nanoparticles)\n",
    "nano_bw = nanoparticles > nano_thresh\n",
    "io.imshow(nano_bw)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: red\"> Task 5:</span> By plotting only smaller sections of the binary image, or by plotting the intensity profile across the image, estimate the size of the particles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: blue\"> Solution:</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Measuring particles by labeling\n",
    "from skimage import measure\n",
    "Now that we have segmented the image we need to find a way of measuring each particle. To do this we turn back to <code>skimage</code> and its <code>measure</code> subpackage:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import measure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function <code>label</code> finds all the separate regions in the image and labels them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tags=measure.label(nano_bw)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see how well the labelling worked using by plotting the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(tags, cmap='Paired')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It seems like it's done a reasonable job but not perfect. There are a number of very small regions that have been labeled as particles when they are just holes in the particles. There are ways to fix this but for now, we'll just have to remember that there are some small labels we need to ignore.\n",
    "\n",
    "The tags array has the same size as the original image but with labels where the particles are. We can see what these labels are by using the <code>unique</code> function of <code>numpy</code>:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels=np.unique(tags)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The label function has identified 44 particles. The first label is for the background:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "io.imshow(tags==0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now measure the particles individually. The particles can be selected by using the comparator <code>==</code>. For example the particle labeled 28:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "io.imshow((tags==28))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the size of particle 28? As you can see above, using <code>==</code> gives a boolean array which is only true where the selected tag exists. We need to convert it into a integer array made up of 0's and 1's, so that the size of the particle (in pixels) is just the sum of the array. We have done this before, when thresholding, we can just use the same method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "particle28=(tags==28).astype(int)\n",
    "particle_area=particle28.sum()\n",
    "print(particle_area)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can convert this area into an \"equivalent diameter\" by remembering that the area of a circle is $A=\\dfrac{\\pi d^2}{4}$. So the \"equivalent radius\" is $d=\\sqrt{\\dfrac{4A}{\\pi}}$ and the equivalent diameter of particle 28 is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "particle_diameter=np.sqrt(4*particle_area/(np.pi))\n",
    "print(particle_diameter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Actually, it is so common to want to know the size of particles that the <code>skimage</code> has a function to do just this. It's called <code>regionprops</code>:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "measurements=measure.regionprops(tags)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One important thing to note is that this function disregards the background. So to measure particle 28, we need to use the index 27:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "measurements[27].equivalent_diameter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the diameters obtained in the two different ways agree almost exactly!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: red\"> Task 6:</span> Calculate the mean and median particle diameter of all the nanoparticles\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: blue\"> Solution:</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: red\"> Bonus task:</span> Can you plot the distribution of diameters using a histogram?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: blue\"> Solution:</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: red\"> Bonus task:</span> Given the particle size distribution and how we calculated it, what is the mean or median particle size most accurate?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color: blue\"> Solution:</span>"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}