Image Processing II¶

(c) 2024 Justin Bois. This work is licensed under a Creative Commons Attribution License CC-BY 4.0. All code contained herein is licensed under an MIT license.

This document was prepared at Caltech with support financial support from the Donna and Benjamin M. Rosen Bioengineering Center.

No description has been provided for this image

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

In [1]:
import numpy as np

import skimage.io
import skimage.exposure
import skimage.filters
import skimage.morphology

import bokeh.io
import bokeh.plotting

import colorcet

import bi1x

bokeh.io.output_notebook()
Loading BokehJS ...

Segmentation¶

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.

For our segmentation, we will use the phase and CFP images of bacteria from the first tutorial. To remind ourselves of the images, let's take a look at the phase image.

In [2]:
# Load image
im_phase = skimage.io.imread("ecoli_phase.tif")

# Update interpixel distance for these bacteria
ip_distance = 0.064

# Display the image
bokeh.io.show(
    bi1x.viz.imshow(
        im_phase,
        interpixel_distance=ip_distance,
        length_units="µm",
    )
)

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. This function is included in the bi1x module as bi1x.viz.im_hist(), but I show how to make the plot here so you can see how to use scikit-image to compute the histogram.

In [3]:
# Get the histogram data
hist, bins = skimage.exposure.histogram(im_phase)

p = bokeh.plotting.figure(
    frame_height=300, frame_width=400, x_axis_label="intensity", y_axis_label="count"
)
p.line(bins, hist, line_width=2, line_join="bevel")

bokeh.io.show(p)

We see that there are two peaks to the histogram of the phase image. The peak to the right is brighter, so likely represents the background. Therefore, if we can find where the valley between the two peaks is, we may take pixels with intensity below that value to be bacteria and those above to be background. Eyeballing it, I think this critical pixel value is about 182.

In [4]:
# Add a line where we eyeballed the threshold
p.ray(182, 0, angle=90, angle_units='deg', line_width=2, length=0, color='tomato')

bokeh.io.show(p)

Thresholding¶

The process of taking pixels above or below a certain value is called thresholding. It is one of the simplest ways to segment an image. We call every pixel with a value below 182 part of a bacterium and everything above not part of a bacterium.

In [5]:
# Gray colormap for thresholded image
gray_mapper = bokeh.models.LinearColorMapper(colorcet.gray)

# Threshold value, as obtained by eye
thresh_phase = 182

# Generate thresholded image
im_phase_bw = im_phase < thresh_phase

# Display phase and thresholded image
p_phase = bi1x.viz.imshow(
    im_phase,
    interpixel_distance=ip_distance,
    length_units="µm",
    frame_height=250,
)
p_phase_bw = bi1x.viz.imshow(
    im_phase_bw,
    color_mapper=gray_mapper,
    interpixel_distance=ip_distance,
    length_units="µm",
    frame_height=250,
)

# Link ranges for panning
p_phase_bw.x_range = p_phase.x_range
p_phase_bw.y_range = p_phase.y_range

# Displace
bokeh.io.show(bokeh.layouts.gridplot([p_phase, p_phase_bw], ncols=2))

We can overlay these images to get a good view. To do this, we will make an RGB image, and saturate the green channel where the thresholded image is white.

In [6]:
# Build RGB image by stacking grayscale images
im_phase_rgb = np.dstack(3 * [im_phase / im_phase.max()])

# Saturate green channel wherever there are white pixels in thresh image
im_phase_rgb[im_phase_bw, 1] = 1.0

# Show the result
bokeh.io.show(bi1x.viz.imshow(im_phase_rgb, color_mapper='rgb'))

We see that we did a decent job finding bacteria, but we do not effectively label the bacteria in the middle of colonies. This is because of the "halo" of high intensity signal near boundaries of the bacteria that we get from using phase contrast microscopy.

Using the CFP channel¶

One way around this issue is to use bacteria that constitutively express a fluorescent protein and to segment in using the fluorescent channel. Let's try the same procedure with the CFP channel. First, let's look at the image.

In [7]:
# Load image
im_cfp = skimage.io.imread('ecoli_cfp.tif')

# Display the image
bokeh.io.show(bi1x.viz.imshow(im_cfp))

We see that the bacteria are typically brighter than the background, so this might help us in segmentation.

Filtering noise: the median filter¶

It is strange that there do not appear to be any yellow (indicating high intensity) pixels in the display of the CFP image with the Viridis LUT. This is because there are some "hot" pixels in the image, resulting from noise or some other error in the detector. We can see this if we zoom in on one of the bad pixels.

In [8]:
bokeh.io.show(bi1x.viz.imshow(im_cfp[150:250,450:550]))

We see a single bright pixel. This will throw off our lookup table. We can remove this noise by using a median filter. The concept is simple. We take a shape of pixels, called a structuring element, and pass it over the image. The value of the center pixel in the max is replaced by the median value of all pixels within the structuring element. To do this, we first need to construct a mask. This is done using the skimage.morphology module. The filtering is then done using skimage.filters.median().

Let's try it with a 3$\times$3 square mask.

In [9]:
# Make the structuring element
selem = skimage.morphology.square(3)

# Perform the median filter
im_cfp_filt = skimage.filters.median(im_cfp, selem)

# Display image
bokeh.io.show(bi1x.viz.imshow(im_cfp_filt))

Now that we have dealt with the noisy pixels, we can now see more clearly that some cells are very bright (shown in yellow) compared with others. We also have an image that makes more sense; we have eliminated the noise.

It is important to note that several cameras in the Bi 1x lab, especially the Flir cameras, have many hot pixels, and median filtering is an important step in processing those images.

Thresholding in the CFP channel¶

We'll proceed by plotting the histogram and finding the threshold value. Eyeballing it, I get a threshold value of 196.

In [10]:
p = bi1x.viz.im_hist(im_cfp_filt)
p.ray(196, 0, angle=90, angle_units='deg', line_width=2, length=0, color='tomato')

bokeh.io.show(p)

Now let's try thresholding this image.

In [11]:
# Threshold value, as obtained by eye
thresh_cfp = 196

# Generate thresholded image
im_cfp_bw = im_cfp_filt > thresh_cfp

# Display phase and thresholded image
p_cfp = bi1x.viz.imshow(
    im_cfp, interpixel_distance=ip_distance, length_units="µm", frame_height=250
)
p_cfp_bw = bi1x.viz.imshow(
    im_cfp_bw,
    color_mapper=gray_mapper,
    interpixel_distance=ip_distance,
    length_units="µm",
    frame_height=250,
)

# Link ranges for panning
p_cfp_bw.x_range = p_cfp.x_range
p_cfp_bw.y_range = p_cfp.y_range

# Displace
bokeh.io.show(bokeh.layouts.gridplot([p_cfp, p_cfp_bw], ncols=2))

Looks like we're doing much better! Let's try overlapping the images now.

In [12]:
# Build RGB image by stacking grayscale images
im_rgb = np.dstack(3 * [im_phase / im_phase.max()])

# Saturate green channel wherever there are white pixels in thresh image
im_rgb[im_cfp_bw, 1] = 1.0

# Show the result
bokeh.io.show(
    bi1x.viz.imshow(
        im_rgb, color_mapper="rgb", interpixel_distance=ip_distance, length_units="µm"
    )
)

Much better, though we see that we overcount the spaces between bacteria.

Otsu's method for thresholding¶

It turns out that there are many automated ways to find the threshold value, as opposed to eyeballing it like we have been doing. Otsu's method provides this functionality.

In [13]:
# Compute Otsu thresholds for phase and cfp
thresh_phase_otsu = skimage.filters.threshold_otsu(im_phase)
thresh_cfp_otsu = skimage.filters.threshold_otsu(im_cfp_filt)

# Compare results to eyeballing it
print('Phase by eye: ', thresh_phase, '   CFP by eye: ', thresh_cfp)
print('Phase by Otsu:', thresh_phase_otsu,   
      '   CFP by Otsu:', thresh_cfp_otsu)

Phase by eye:  182    CFP by eye:  196
Phase by Otsu: 248    CFP by Otsu: 195

We see that for the CFP channel, the Otsu method did very well. However, for phase, we see a big difference. This is because the Otsu method assumes a bimodal distribution of pixels. If we look at the histograms on a log scale, we see more clearly that the phase image has a long tail, which will trip up the Otsu algorithm. The moral of the story is that you can use automated thresholding, but you should always do sanity checks to make sure it is working as expected.

In [14]:
# Plot the histograms together.
p = bi1x.viz.im_hist(im_phase, y_axis_type='log', legend_label='phase')
p = bi1x.viz.im_hist(im_cfp_filt, p=p, color='orange', legend_label='CFP')

bokeh.io.show(p)

Determining the bacterial area¶

Now that we have a thresholded image, we can determine the total area taken up by bacteria. It's as simple as summing up the pixel values of the thresholded image!

In [15]:
# Compute bacterial area
bacterial_area_pix = im_cfp_bw.sum()

# Print out the result
print('bacterial area =', bacterial_area_pix, 'pixels')

bacterial area = 264437 pixels

For growth curves, we really only need the pixel values, since we are trying to get a time constant, which is independent of the unit of measure for bacterial numbers. Nonetheless, if we want to get the total area that is bacterial in units of square µm, we could use the interpixel distances to get the area represented by each pixel. For this setup, the interpixel distance is 0.0636 µm. We can then compute the bacterial area as follows.

In [16]:
# Compute bacterial area
bacterial_area_micron = bacterial_area_pix * ip_distance**2

# Print total area
print('bacterial area =', bacterial_area_micron, 'square microns')

bacterial area = 1083.133952 square microns

Computing environment¶

In [17]:
%load_ext watermark
%watermark -v -p numpy,skimage,bi1x,bokeh,colorcet,jupyterlab

Python implementation: CPython
Python version       : 3.12.2
IPython version      : 8.20.0

numpy     : 1.26.4
skimage   : 0.22.0
bi1x      : 0.0.14
bokeh     : 3.4.0
colorcet  : 3.1.0
jupyterlab: 4.0.11