Lesson 34: Introduction to image processing with scikit-image

This tutorial was generated from a Jupyter notebook. You can download the notebook here.

In [21]:
# Our workhorse
import numpy as np

# Our image processing tools
import skimage.filters
import skimage.io
import skimage.morphology
import skimage.exposure

# This is how we import the module of Matplotlib we'll be using
import matplotlib.pyplot as plt

# Seaborn makes plots pretty!
import seaborn as sns

# Magic function to make matplotlib inline; other style specs must come AFTER
%matplotlib inline

# This enables SVG graphics inline (only use with static plots (non-Bokeh))
%config InlineBackend.figure_format = 'svg'

# Set JB's favorite Seaborn settings
rc={'lines.linewidth': 2, 'axes.labelsize': 18, 'axes.titlesize': 18, 
    'axes.facecolor': 'DFDFE5'}
sns.set_context('notebook', rc=rc)

In this tutorial, we will learn some basic techniques for image processing using scikit-image with Python.

Image processing tools for Python

There are many image processing tools available for Python. Some of them, such as ITK and OpenCV are mature image processing packages that have bindings for Python, allowing easy use of their functionality. Others were developed specifically for Python. Some of the many packages are

The first two packages are standard with Anaconda. They provide a set of basic image processing tools, with more sophisticated packages such as ITK and Fiji supplying many more bells and whistles. If in the future you have more demanding image processing requirements, the other packages can prove very useful.

We will almost exclusively use scikit-image along with the standard tools from NumPy. The package scipy.ndimage is quite useful, but we will use scikit-image, since it has expanded functionality. We will call is skimage for short (which is how the package is imported anyhow). A potential annoyance with skimage is that the main package has minimal functionality, and you must import subpackages as needed. For example, to load and view images, you will need to import skimage.io. Importantly, skimage is well-documented, and you can access the documentation at http://scikit-image.org/.

We will explore skimage’s capabilities and some basic image processing techniques through example. In this lesson, we will take a brightfield and a fluorescent image of bacteria and perform segementation, that is the identification of each pixel in the image as being bacteria or background.

Loading and viewing 12-bit images

We will now load and view the test images we will use for segmentation. We load the image using the skimage.io.imread. The image is stored as a NumPy array. Each entry in the array is a pixel value. This is an important point: a digital image is data. It is a set of numbers with spatial positions.

You should have downloaded the images ecoli_phase.tif and ecoli_cfp.tif. We'll start with an image taken using phase contrast microscopy of E. coli. These images were acquired by Caltech graduate student Griffin Chure.

In [22]:
# Load image
im_phase = skimage.io.imread('../data/ecoli_phase.tif')

# Display the image, set Seaborn style 'dark' to avoid grid lines
with sns.axes_style('dark'):
    skimage.io.imshow(im_phase)
/Users/Justin/anaconda/lib/python3.4/site-packages/skimage/io/_plugins/matplotlib_plugin.py:74: UserWarning: Low image dynamic range; displaying image with stretched contrast.
  warnings.warn("Low image dynamic range; displaying image with "

We get a warning and then a strange looking picture! The warning is that we have a low dynamic range. We get this warning because skimage assumes that images with integer pixel values are either 8-bit, 16-bit, or 32-bit. This image, as is often the case with images acquired with scientific cameras, were acquired with a 12-bit camera. This means that the pixel values range from 0 to $2^{12}-1=4095$. This is much less than what we would expect for a maximal pixel value for a 16-bit image, $2^{16}-1 = 65535$.

So, to view the images in their full range, we should divide by the maximal pixel value.

In [23]:
# Load image
im_phase = skimage.io.imread('../data/ecoli_phase.tif')

# Display the image
with sns.axes_style('dark'):
    skimage.io.imshow(im_phase / im_phase.max())

Lookup tables

As an alternative to stretching the image and using skimage.io.imshow to view the image, we can use matploblib's image viewing function, which automatically adjusts the image. We will need to specify a gray colormap to look at it. In fact, while we're at it, we can specify whatever colormap we want.

In [24]:
with sns.axes_style('dark'):
    # Get subplots
    fig, ax = plt.subplots(2, 2, figsize=(8,6))

    # Display various LUTs
    ax[0,0].imshow(im_phase, cmap=plt.cm.gray)
    ax[0,1].imshow(im_phase, cmap=plt.cm.RdBu_r)
    ax[1,0].imshow(im_phase, cmap=plt.cm.YlGnBu_r)
    ax[1,1].imshow(im_phase, cmap=plt.cm.copper)

We did a few new things here. First, we use the plt.subplots function to generate subplots. Second, we used the imshow method of the axes object to display the images. Lastly, we specified the colormap to be used for showing the image. We used a grayscale colormap, plt.cm.gray. This specifies how the pixel values are interpreted as they are displayed. For a grayscale colormap, high pixel values are more white, and low pixel values are more black.

In image processing, a colormap is called a lookup table (LUT). A LUT is a mapping of pixel values to a color. This sometimes helps visualize images, especially when we use false coloring. Remember, a digital image is data, and false coloring an image if not manipulation of data. It is simply a different way of plotting it.

As we just saw, we specify a lookup table with a colormap. There are plenty available in matplotlib. There is lots of debate about that the best colormaps (LUTs) are. The data visualization community seems to universally reject using rainbow colormaps. See, e.g., D. Borland and R. M. Taylor, Rainbow Color Map (Still) Considered Harmful, IEEE Computer Graphics and Applications, 27,14-17, 2007. In the example, I use a hue-diverging colorscale, which goes from blue to red, as people accustomed to rainbow colormaps expect, but in a perceptually ordered fashion. In the soon-to-be released Matplotlib 2.0, the colormap viridis will be available and is beautiful!

Importantly, the false coloring helps use see that the intensity of the pixel values near the center of some of the bacteria are not that dissimilar from the background. This will become an issue, as we will see, when segmenting.

Introductory segmentation

As mentioned before, segmentation is the process by which we separate regions of an image according to their identity for easier analysis. E.g., if we have an image of bacteria and we want to determine what is "bacteria" and what is "not bacteria," we would do some segmentation. We will use bacterial test images for this purpose.

Histograms

As we begin segmentation, remember that viewing an image is just a way of plotting the digital image data. We can also plot a histogram. This helps use see some patterns in the pixel values and is often an important first step toward segmentation.

The histogram of an image is simply a list of counts of pixel values. When we plot the histogram, we can often readily see breaks in which pixel values are most frequently encountered. There are many ways of looking at histograms. I’ll show you my preferred way.

In [25]:
# Get the histogram data
hist_phase, bins_phase = skimage.exposure.histogram(im_phase)

# Use matplotlib to make a pretty plot of histogram data
plt.fill_between(bins_phase, hist_phase, alpha=0.5)

# Label axes
plt.xlabel('pixel value')
plt.ylabel('count')
Out[25]:
<matplotlib.text.Text at 0x11f52b588>