Biomedical Image Analysis in Python (Part - 2).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1>Biomedical Image Analysis in Python (Part - 2)</h1>\n",
    "<h3>by Sakib Reza<h3>\n",
    "<h2>Intensity</h2>\n",
    "we will work with a hand <b>radiograph</b> from a 2017 <b>Radiological Society of North America competition</b>. X-ray absorption is highest in dense tissue such as bone, so the resulting intensities should be high. Consequently, images like this can be used to predict \"bone age\" in children.<br>\n",
    "To start, let's load the image and check its intensity range.<br>\n",
    "The image datatype determines the range of possible intensities: e.g., 8-bit unsigned integers (uint8) can take values in the range of 0 to 255. A colorbar can be helpful for connecting these values to the visualized image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import imageio\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Load the hand radiograph\n",
    "im = imageio.imread(\"hand-xray.JPG\")\n",
    "print('Data type:', im.dtype)\n",
    "print('Min. value:', im.min())\n",
    "print('Max value:', im.max())\n",
    "\n",
    "# Plot the grayscale image\n",
    "plt.imshow(im, cmap='gray', vmin=0, vmax=255)\n",
    "plt.colorbar()\n",
    "format_and_render_plot()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Output:</b><br>\n",
    "Data type: uint8<br>\n",
    "Min. value: 1<br>\n",
    "Max value: 230<br>\n",
    "<img src=\"img5.JPG\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Histograms</h2>\n",
    "Histograms display the distribution of values in your image by binning each element by its intensity then measuring the size of each bin.\n",
    "\n",
    "The area under a histogram is called the cumulative distribution function. It measures the frequency with which a given range of pixel intensities occurs.\n",
    "\n",
    "Here we will describe the intensity distribution in im by calculating the histogram and cumulative distribution function and displaying them together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import SciPy's \"ndimage\" module\n",
    "import scipy.ndimage as ndi \n",
    "\n",
    "# Create a histogram, binned at each possible value\n",
    "hist = ndi.histogram(im, min=0, max=255, bins=256)\n",
    "\n",
    "# Create a cumulative distribution function\n",
    "cdf = hist.cumsum() / hist.sum()\n",
    "\n",
    "# Plot the histogram and CDF\n",
    "fig, axes = plt.subplots(2, 1, sharex=True)\n",
    "axes[0].plot(hist, label='Histogram')\n",
    "axes[1].plot(cdf, label='CDF')\n",
    "format_and_render_plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Output:</b>\n",
    "<img src=\"img6.JPG\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Create a mask</h2>\n",
    "Masks are the primary method for removing or selecting specific parts of an image. They are binary arrays that indicate whether a value should be included in an analysis. Typically, masks are created by applying one or more logical operations to an image.\n",
    "\n",
    "We will try to use a simple intensity threshold to differentiate between skin and bone in the hand radiograph. (im has been equalized to utilize the whole intensity range.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create skin and bone masks\n",
    "mask_bone = im>=145\n",
    "mask_skin = (im>=45) & (im<145)\n",
    "\n",
    "# Plot the skin (0) and bone (1) masks\n",
    "fig, axes = plt.subplots(1,2)\n",
    "axes[0].imshow(mask_skin, cmap='gray')\n",
    "axes[1].imshow(mask_bone, cmap='gray')\n",
    "format_and_render_plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Output:</b>\n",
    "<img src=\"img7.JPG\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Apply a mask</h2>\n",
    "Although masks are binary, they can be applied to images to filter out pixels where the mask is False.\n",
    "\n",
    "NumPy's where() function is a flexible way of applying masks. It takes three arguments:\n",
    "\n",
    "np.where(condition, x, y)<br><br>\n",
    "condition, x and y can be either arrays or single values. This allows you to pass through original image values while setting masked values to 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import SciPy's \"ndimage\" module\n",
    "import scipy.ndimage as ndi \n",
    "\n",
    "# Screen out non-bone pixels from \"im\"\n",
    "mask_bone = im>=145\n",
    "im_bone = np.where(mask_bone, im, 0)\n",
    "\n",
    "# Get the histogram of bone intensities\n",
    "hist = ndi.histogram(im_bone, min=1, max=255, bins=255)\n",
    "\n",
    "# Plot masked image and histogram\n",
    "fig, axes = plt.subplots(2,1)\n",
    "axes[0].imshow(im_bone)\n",
    "axes[1].plot(hist)\n",
    "format_and_render_plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Output:</b>\n",
    "<img src=\"img8.JPG\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Tune a mask</h2>\n",
    "Imperfect masks can be tuned through the addition and subtraction of pixels. SciPy includes several useful methods for accomplishing these ends. These include:\n",
    "\n",
    "binary_dilation: Add pixels along edges<br>\n",
    "binary_erosion: Remove pixels along edges<br>\n",
    "binary_opening: Erode then dilate, \"opening\" areas near edges<br>\n",
    "binary_closing: Dilate then erode, \"filling in\" holes<br><br>\n",
    "Here we will create a bone mask then tune it to include additional pixels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create and tune bone mask\n",
    "mask_bone = im>=145\n",
    "mask_dilate = ndi.binary_dilation (mask_bone, iterations=5)\n",
    "mask_closed = ndi.binary_closing (mask_bone, iterations=5)\n",
    "\n",
    "# Plot masked images\n",
    "fig, axes = plt.subplots(1,3)\n",
    "axes[0].imshow(mask_bone)\n",
    "axes[1].imshow(mask_dilate)\n",
    "axes[2].imshow(mask_closed)\n",
    "format_and_render_plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Output:</b>\n",
    "<img src=\"img9.JPG\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Filter Convolutions</h2>\n",
    "Filters are an essential tool in image processing. They allow you to transform images based on intensity values surrounding a pixel, rather than globally.\n",
    "\n",
    "Here we will smooth the foot radiograph. First, specify the weights to be used. (These are called \"footprints\" and \"kernels\" as well.) Then, convolve the filter with im and plot the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set filter weights\n",
    "weights = [[0.11, 0.11, 0.11],\n",
    "           [0.11, 0.11, 0.11], \n",
    "           [0.11, 0.11, 0.11]]\n",
    "\n",
    "# Convolve the image with the filter\n",
    "im_filt = ndi.convolve(im, weights)\n",
    "\n",
    "# Plot the images\n",
    "fig, axes = plt.subplots(1,2)\n",
    "axes[0].imshow(im)\n",
    "axes[1].imshow(im_filt)\n",
    "format_and_render_plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Output:</b>\n",
    "<img src=\"img10.JPG\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Smoothing</h2>\n",
    "Smoothing can improve the signal-to-noise ratio of your image by blurring out small variations in intensity. The Gaussian filter is excellent for this: it is a circular (or spherical) smoothing kernel that weights nearby pixels higher than distant ones.\n",
    "\n",
    "The width of the distribution is controlled by the sigma argument, with higher values leading to larger smoothing effects.\n",
    "\n",
    "Here we will test the effects of applying Gaussian filters to the foot x-ray before creating a bone mask."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Smooth \"im\" with Gaussian filters\n",
    "im_s1 = ndi.gaussian_filter(im, sigma=1)\n",
    "im_s3 = ndi.gaussian_filter(im, sigma=3)\n",
    "\n",
    "# Draw bone masks of each image\n",
    "fig, axes = plt.subplots(1,3)\n",
    "axes[0].imshow(im >= 145)\n",
    "axes[1].imshow(im_s1 >= 145)\n",
    "axes[2].imshow(im_s3 >= 145)\n",
    "format_and_render_plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Output:</b>\n",
    "<img src=\"img11.JPG\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Detect Edges (1)</h2>\n",
    "Filters can also be used as \"detectors.\" If a part of the image fits the weighting pattern, the returned value will be very high (or very low).\n",
    "\n",
    "In the case of edge detection, that pattern is a change in intensity along a plane. Here we will create a vertical edge detector and see how well it performs on the hand x-ray (im)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set weights to detect vertical edges\n",
    "weights = [[+1, 0, -1], [+1, 0, -1], [+1, 0, -1]]\n",
    "\n",
    "# Convolve \"im\" with filter weights\n",
    "edges = ndi.convolve(im, weights)\n",
    "\n",
    "# Draw the image in color\n",
    "plt.imshow(edges, cmap='seismic', vmin=-150, vmax=150)\n",
    "plt.colorbar()\n",
    "format_and_render_plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Output:</b>\n",
    "<img src=\"img12.JPG\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Detect Edges (2)</h2>\n",
    "Edge detection can be performed along multiple axes, then combined into a single edge value. For 2D images, the horizontal and vertical \"edge maps\" can be combined using the Pythagorean theorem:\n",
    "\n",
    "z=√(x^2+y^2)\n",
    "\n",
    "One popular edge detector is the Sobel filter. The Sobel filter provides extra weight to the center pixels of the detector:\n",
    "\n",
    "weights = [[ 1,  2,  1], \n",
    "           [ 0,  0,  0],\n",
    "           [-1, -2, -1]]\n",
    "Here we will improve upon our previous detection effort by merging the results of two Sobel-filtered images into a composite edge map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply Sobel filter along both axes\n",
    "sobel_ax0 = ndi.sobel(im, axis=0)\n",
    "sobel_ax1 = ndi.sobel(im, axis=1)\n",
    "\n",
    "# Calculate edge magnitude \n",
    "edges = np.sqrt(np.square(sobel_ax0) + np.square(sobel_ax1))\n",
    "\n",
    "# Plot edge magnitude\n",
    "plt.imshow(edges, cmap='gray', vmax=75)\n",
    "format_and_render_plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Output:</b>\n",
    "<img src=\"img13.JPG\">"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}