Segmenting cell images (fluorescent microscopy).ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Segmenting cell images\n",
      "\n",
      "We start with the same imports as last time"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import mahotas as mh\n",
      "import numpy as np\n",
      "from matplotlib import pyplot as plt\n",
      "from IPython.html.widgets import interact, fixed\n",
      "plt.rcParams['figure.figsize'] = (10.0, 8.0) # 10 x 8 inches\n",
      "plt.gray()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "mahotas ships with an image of nuclei for our analysis.\n",
      "\n",
      "Because of file format issues (many image loaders cannot correctly handle single channel PNG images), it is encoded as an RGB image, where all the channels have the same value. So, we convert it to a single channel image (using `max`, but `min` would give the same result):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dna = mh.demos.load('nuclear')\n",
      "print(dna.shape)\n",
      "dna = dna.max(axis=2)\n",
      "print(dna.shape)\n",
      "plt.imshow(dna)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Thresholding\n",
      "\n",
      "The first thing we try is Otsu thresholding, a venerable method in image processing:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "T_otsu = mh.otsu(dna)\n",
      "print(T_otsu)\n",
      "plt.imshow(dna > T_otsu)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Does not work so well on this data. What does work well is thresholding by the mean (this actually works well not just on this image, but for other similar images too, see [this paper of mine](http://ieeexplore.ieee.org/xpl/login.jsp?tp=&arnumber=5193098&url=http%3A%2F%2Fieeexplore.ieee.org%2Fxpls%2Fabs_all.jsp%3Farnumber%3D5193098):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "T_mean = dna.mean()\n",
      "print(T_mean)\n",
      "plt.imshow(dna > T_mean)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Only issue is that the results are very noisy around the edges, but we can add some Gaussian blurring to correct for this:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dnaf = mh.gaussian_filter(dna, 2.)\n",
      "T_mean = dnaf.mean()\n",
      "bin_image = dnaf > T_mean\n",
      "plt.imshow(bin_image)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Right now, we have a binary image, but it is useful to [label the image](http://en.wikipedia.org/wiki/Connected-component_labeling) in order to assign a different (integer) index to each component:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "labeled, nr_objects = mh.label(bin_image)\n",
      "print(nr_objects)\n",
      "\n",
      "plt.imshow(labeled)\n",
      "plt.jet()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Separating touching cells\n",
      "\n",
      "We now got the cells nicely separated from the background, but we have also merged several cells together.\n",
      "\n",
      "The general strategy we can use to separate the cells is the following:\n",
      "\n",
      "1. Smooth the image with a Gaussian filter (we need to specify the  $\\sigma$ parameter)\n",
      "2. Find regional maxima on this smoothed image to identify each cell\n",
      "3. Use [watershed](http://en.wikipedia.org/wiki/Watershed_(image_processing%29) on the [distance transformed image](http://en.wikipedia.org/wiki/Distance_transform) to separate cells.\n",
      "\n",
      "We are going to use an interactive mode function to find out the best value for $\\sigma$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "@interact(sigma=(1.,16.))\n",
      "def check_sigma(sigma):\n",
      "    dnaf = mh.gaussian_filter(dna.astype(float), sigma)\n",
      "    maxima = mh.regmax(mh.stretch(dnaf))\n",
      "    maxima = mh.dilate(maxima, np.ones((5,5)))\n",
      "    plt.imshow(mh.as_rgb(np.maximum(255*maxima, dnaf), dnaf, dna > T_mean))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "It seems $\\sigma = 12$ is a pretty good guess:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sigma = 12.0\n",
      "dnaf = mh.gaussian_filter(dna.astype(float), sigma)\n",
      "maxima = mh.regmax(mh.stretch(dnaf))\n",
      "maxima,_= mh.label(maxima)\n",
      "plt.imshow(maxima)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, we compute the distance transform. Again, this is a single call, `mh.distance`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dist = mh.distance(bin_image)\n",
      "plt.imshow(dist)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Because of the way that the watershed function is defined in mahotas, we need to invert the distance transform. Also, for technical reasons, we convert to `uint8`.\n",
      "\n",
      "Finally, we can call `mh.cwatershed` with the `dist` image and the `maxima` as seeds."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dist = 255 - mh.stretch(dist)\n",
      "watershed = mh.cwatershed(dist, maxima)\n",
      "plt.imshow(watershed)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We now take away the background to make things nicer:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "watershed *= bin_image\n",
      "plt.imshow(watershed)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Cleaning up regions\n",
      "\n",
      "Almost there.\n",
      "\n",
      "Let's first remove cells touching the border (the module `mh.labeled` has [several functions](http://mahotas.readthedocs.org/en/latest/labeled.html?highlight=labeled) to deal with labeled images):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "watershed = mh.labeled.remove_bordering(watershed)\n",
      "plt.imshow(watershed)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let us also remove things that are too small to be a cell. Again, we use an interactive display to choose the parameters:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sizes = mh.labeled.labeled_size(watershed)\n",
      "\n",
      "# The conversion below is not necessary in newer versions of mahotas:\n",
      "watershed = watershed.astype(np.intc)\n",
      "\n",
      "@interact(min_size=(100,4000,20))\n",
      "def do_plot(min_size):\n",
      "    filtered = mh.labeled.remove_regions_where(watershed, sizes < min_size)\n",
      "    print(\"filtering {}...\".format(min_size))\n",
      "    plt.imshow(filtered)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "2000 seems like it works.\n",
      "\n",
      "So, we remove these regions and relabel.\n",
      "\n",
      "Relabeling is necessary because the `remove_regions_where` call only sets the removed regions to zero, relabeling is necessary to make the result be in the range `0 .. N`. Relabel also returns the number of objects that are left."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "min_size = 2000\n",
      "filtered = mh.labeled.remove_regions_where(watershed, sizes < min_size)\n",
      "\n",
      "labeled,nr_objects = mh.labeled.relabel(filtered)\n",
      "print(\"Number of cells: {}\".format(nr_objects))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}