Contents

4. Thresholding

4. Thresholding

When performing logical operations of the type image > value we have seen that this allowed us to perform thresholding, one the basic image processing methods to identify objects in an image. Here we come back to this topic and show how we can automated the threshold selection, combine masks etc.

import numpy as np
import matplotlib.pyplot as plt
import skimage
import skimage.io

Importing the image

Let’s import the image, plot it and check its dimensions:

#image_stack = skimage.io.imread('images/46658_784_B12_1.tif')
image_stack = skimage.io.imread('https://github.com/guiwitz/PyImageCourse_beginner/raw/master/images/46658_784_B12_1.tif')

image_stack.shape

(2048, 2048, 3)

plt.subplots(figsize=(10,10))
plt.imshow(image_stack);
_images/04-Thresholding_7_0.png

We have again a three channel image, i.e. our image is composed of three planes. For the moment we only want to work with the nuclei, so we isolate that channel, the third one:

image_nuclei = image_stack[:,:,2]

plt.subplots(figsize=(10,10))
plt.imshow(image_nuclei, cmap = 'gray');
_images/04-Thresholding_10_0.png

9.2 Manual investigation

We have seen in the previous chapter that if we want to isolate nuclei, we have to know more or less what intensity they have so that we can use a cutoff value on pixel intensities. We can look at the histogram of pixel values:

plt.hist(np.ravel(image_nuclei));
_images/04-Thresholding_13_0.png

Clearly we have background intensity between 0 and ~50, and the signal spread out between 50 and 255. We can create a mask of all the pixels with a value above 50:

mask = image_nuclei > 50

plt.imshow(mask, cmap = 'gray');
_images/04-Thresholding_16_0.png

For the sake of the experiment, let’s use a higher threshold:

mask = image_nuclei > 120
plt.imshow(mask, cmap = 'gray');
_images/04-Thresholding_18_0.png

Clearly, if the threshold is too high, we loose import parts of the image.

Automatic thresholding

Setting a threshold manually is not ideal for two reasons:

  1. your eyes might fool you

  2. if you have to do this operation on 100s of images you don’t want to do it by hand.

Several methods have been developed to detect an “ideal” threshold based on the pixel histogram. One of the most famous ones is the Otsu method: it sets a threshold by dividing the pixels in two classes and tries to minimize the variance of both classes.

Such classical algorithms are all implemented in scikit-image and we can just use them out of the box. All algorithm are grouped in logical groups (filtering, measuring etc.) Often the easiest is just to use Google to find them.

For example the Otsu algorithm:

import skimage.filters

help(skimage.filters.threshold_otsu)

Help on function threshold_otsu in module skimage.filters.thresholding:

threshold_otsu(image=None, nbins=256, *, hist=None)
    Return threshold value based on Otsu's method.
    
    Either image or hist must be provided. If hist is provided, the actual
    histogram of the image is ignored.
    
    Parameters
    ----------
    image : (N, M[, ..., P]) ndarray, optional
        Grayscale input image.
    nbins : int, optional
        Number of bins used to calculate histogram. This value is ignored for
        integer arrays.
    hist : array, or 2-tuple of arrays, optional
        Histogram from which to determine the threshold, and optionally a
        corresponding array of bin center intensities. If no hist provided,
        this function will compute it from the image.
    
    
    Returns
    -------
    threshold : float
        Upper threshold value. All pixels with an intensity higher than
        this value are assumed to be foreground.
    
    References
    ----------
    .. [1] Wikipedia, https://en.wikipedia.org/wiki/Otsu's_Method
    
    Examples
    --------
    >>> from skimage.data import camera
    >>> image = camera()
    >>> thresh = threshold_otsu(image)
    >>> binary = image <= thresh
    
    Notes
    -----
    The input image must be grayscale.

my_otsu_threshold = skimage.filters.threshold_otsu(image_nuclei)
print(my_otsu_threshold)

66

Using that value we can create a new mask:

mask_nuclei = image_nuclei > my_otsu_threshold

plt.imshow(mask_nuclei, cmap = 'gray');
_images/04-Thresholding_26_0.png

Combining masks

We can do the same thresholding operation on the other channel contained in the image. Let’s first “isolate” the first channel (red):

image_cells = image_stack[:,:,0]

plt.subplots(figsize=(10,10))
plt.imshow(image_cells, cmap = 'gray');
_images/04-Thresholding_30_0.png 
mask_cells = image_cells > skimage.filters.threshold_otsu(image_cells)

plt.subplots(figsize=(10,10))
plt.imshow(mask_cells, cmap = 'gray');
_images/04-Thresholding_32_0.png

We can in principle combine masks however we want. For example we can consider pixels that are ON (True) in both masks by multiplying them as only True x True = True

both_masks = mask_cells * mask_nuclei

plt.subplots(figsize=(10,10))
plt.imshow(both_masks, cmap = 'gray');
_images/04-Thresholding_35_0.png

or we can consider pixels that are ON (True) in one or the other image by adding them as True + False = True:

plt.subplots(figsize=(10,10))
plt.imshow(mask_cells + mask_nuclei, cmap = 'gray');
_images/04-Thresholding_37_0.png

To make things cleaner, instead of just using operators such as + or * we can be more explicit and use logical Numpy operations such as np.logical_or():

combine_or = np.logical_or(mask_cells, mask_nuclei)
plt.imshow(combine_or, cmap = 'gray');
_images/04-Thresholding_39_0.png

Now that we have segmented elements in our image, we want to be able to quantify those elements, which is what is done in chapter 11. In part 10 we can see how we can improve the masks by cleaning them up using binary/morphological operations.

previous

3. More on images

next

5. Filtering