Implementing Image Segmentation Techniques using Scikit-Image

20 September, 2020

Report –

Experiments –

Table of Contents

Implementing Image Segmentation Techniques using Scikit-ImageIntroductionSetupExperiments and Results1. Edge Detection using Roberts, Sobel and Prewitt Operators2. Canny Edge Detection3. ThresholdingBimodal HistogramOtsu’s ThresholdingLocal Thresholding4. Multi-Otsu Thresholding5. Active Contour ModelConclusionReferences

Introduction

Image segmentation is a critical process in computer vision. It involves dividing a visual input into segments to simplify image analysis. Segments represent objects or parts of objects, and comprise sets of pixels. Image segmentation sorts pixels into larger components, eliminating the need to consider individual pixels as units of observation.

scikit-image is a collection of algorithms for image processing. It is available free of charge and free of restriction. We pride ourselves on high-quality, peer-reviewed code, written by an active community of volunteers.

This report is an account of using Scikit-Image – an open source Python package designed for image preprocessing, to experiment with some of the basic Image Segmentation Techniques as listed below –

Edge Detection using Roberts, Sobel and Prewitt Operators
Canny Edge Detection
Thresholding - local and global
Multi-Otsu Thresholding
Active Contour Model

Setup

All the programs and code works were executed on Jupyter environment on my local machine.

The Jupyter Notebook is an open-source web application that allows you to create and share documents that contain live code, equations, visualizations and narrative text. Uses include: data cleaning and transformation, numerical simulation, statistical modeling, data visualization, machine learning, and much more.

Alternatively, one can opt for Google Colabaortory or Binder to reproduce the same experiments.

Experiments and Results

Importing the necessary libraries of Scikit-Image, NumPy and Matplotlib’s pyplot.


xxxxxxxxxx
3
1
import skimage
2
import numpy as np
3
import matplotlib.pyplot as plt

1. Edge Detection using Roberts, Sobel and Prewitt Operators

Edge operators are used in image processing within edge detection algorithms. They are discrete differentiation operators, computing an approximation of the gradient of the image intensity function. Different operators compute different finite-difference approximations of the gradient.


xxxxxxxxxx
6
1
from skimage import filters
2
from skimage.data import camera
3
from skimage.util import compare_images
4
5
image = camera()   #from skimage library
6
plt.imshow(image)

png


xxxxxxxxxx
21
1
edge_roberts = filters.roberts(image)
2
edge_sobel = filters.sobel(image)
3
edge_prewitt = filters.prewitt(image)
4
5
fig, axes = plt.subplots(ncols=3, sharex=True, sharey=True,
6
                         figsize=(8, 8))
7
8
axes[0].imshow(edge_roberts, cmap=plt.cm.gray)
9
axes[0].set_title('Roberts Edge Detection')
10
11
axes[1].imshow(edge_sobel, cmap=plt.cm.gray)
12
axes[1].set_title('Sobel Edge Detection')
13
14
axes[2].imshow(edge_prewitt, cmap=plt.cm.gray)
15
axes[2].set_title('Prewitt Edge Detection')
16
17
for ax in axes:
18
    ax.axis('off')
19
20
plt.tight_layout()
21
plt.show()

png

2. Canny Edge Detection

The Canny filter is a multi-stage edge detector. It uses a filter based on the derivative of a Gaussian in order to compute the intensity of the gradients.The Gaussian reduces the effect of noise present in the image. Then, potential edges are thinned down to 1-pixel curves by removing non-maximum pixels of the gradient magnitude. Finally, edge pixels are kept or removed using hysteresis thresholding on the gradient magnitude.

The Canny has three adjustable parameters: the width of the Gaussian (the noisier the image, the greater the width), and the low and high threshold for the hysteresis thresholding.


xxxxxxxxxx
35
1
from scipy import ndimage as ndi
2
from skimage import feature
3
4
5
# Generate noisy image of a square
6
im = np.zeros((128, 128))
7
im[32:-32, 32:-32] = 1
8
9
im = ndi.rotate(im, 15, mode='constant')
10
im = ndi.gaussian_filter(im, 4)
11
im += 0.2 * np.random.random(im.shape)
12
13
# Compute the Canny filter for two values of sigma
14
edges1 = feature.canny(im)
15
edges2 = feature.canny(im, sigma=3)
16
17
# display results
18
fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, figsize=(8, 3),
19
                                    sharex=True, sharey=True)
20
21
ax1.imshow(im, cmap=plt.cm.gray)
22
ax1.axis('off')
23
ax1.set_title('noisy image', fontsize=20)
24
25
ax2.imshow(edges1, cmap=plt.cm.gray)
26
ax2.axis('off')
27
ax2.set_title(r'Canny filter, $\sigma=1$', fontsize=20)
28
29
ax3.imshow(edges2, cmap=plt.cm.gray)
30
ax3.axis('off')
31
ax3.set_title(r'Canny filter, $\sigma=3$', fontsize=20)
32
33
fig.tight_layout()
34
35
plt.show()

png

3. Thresholding

Thresholding is used to create a binary image from a grayscale image. It is the simplest way to segment objects from a background.

Thresholding algorithms implemented in scikit-image can be separated in two categories:

Histogram-based. The histogram of the pixels’ intensity is used and certain assumptions are made on the properties of this histogram (e.g. bimodal).
Local. To process a pixel, only the neighboring pixels are used. These algorithms often require more computation time.

If you are not familiar with the details of the different algorithms and the underlying assumptions, it is often difficult to know which algorithm will give the best results. Therefore, Scikit-image includes a function to evaluate thresholding algorithms provided by the library. At a glance, you can select the best algorithm for you data without a deep understanding of their mechanisms.


xxxxxxxxxx
7
1
from skimage import data
2
from skimage.filters import try_all_threshold
3
4
img = data.page()
5
6
fig, ax = try_all_threshold(img, figsize=(10, 8), verbose=False)
7
plt.show()

png

Now, we illustrate how to apply one of these thresholding algorithms. This example uses the mean value of pixel intensities. It is a simple and naive threshold value, which is sometimes used as a guess value


xxxxxxxxxx
20
1
from skimage.filters import threshold_mean
2
3
4
image = data.camera()
5
thresh = threshold_mean(image)
6
binary = image > thresh
7
8
fig, axes = plt.subplots(ncols=2, figsize=(8, 3))
9
ax = axes.ravel()
10
11
ax[0].imshow(image, cmap=plt.cm.gray)
12
ax[0].set_title('Original image')
13
14
ax[1].imshow(binary, cmap=plt.cm.gray)
15
ax[1].set_title('Result')
16
17
for a in ax:
18
    a.axis('off')
19
20
plt.show()

png

Bimodal Histogram

For pictures with a bimodal histogram, more specific algorithms can be used. For instance, the minimum algorithm takes a histogram of the image and smooths it repeatedly until there are only two peaks in the histogram.


xxxxxxxxxx
25
1
from skimage.filters import threshold_minimum
2
3
4
image = data.camera()
5
6
thresh_min = threshold_minimum(image)
7
binary_min = image > thresh_min
8
9
fig, ax = plt.subplots(2, 2, figsize=(10, 10))
10
11
ax[0, 0].imshow(image, cmap=plt.cm.gray)
12
ax[0, 0].set_title('Original')
13
14
ax[0, 1].hist(image.ravel(), bins=256)
15
ax[0, 1].set_title('Histogram')
16
17
ax[1, 0].imshow(binary_min, cmap=plt.cm.gray)
18
ax[1, 0].set_title('Thresholded (min)')
19
20
ax[1, 1].hist(image.ravel(), bins=256)
21
ax[1, 1].axvline(thresh_min, color='r')
22
23
for a in ax[:, 0]:
24
    a.axis('off')
25
plt.show()

png

Otsu’s Thresholding

Otsu’s method calculates an “optimal” threshold (marked by a red line in the histogram below) by maximizing the variance between two classes of pixels, which are separated by the threshold. Equivalently, this threshold minimizes the intra-class variance”


xxxxxxxxxx
26
1
from skimage.filters import threshold_otsu
2
3
4
image = data.camera()
5
thresh = threshold_otsu(image)
6
binary = image > thresh
7
8
fig, axes = plt.subplots(ncols=3, figsize=(16, 5))
9
ax = axes.ravel()
10
ax[0] = plt.subplot(1, 3, 1)
11
ax[1] = plt.subplot(1, 3, 2)
12
ax[2] = plt.subplot(1, 3, 3, sharex=ax[0], sharey=ax[0])
13
14
ax[0].imshow(image, cmap=plt.cm.gray)
15
ax[0].set_title('Original')
16
ax[0].axis('off')
17
18
ax[1].hist(image.ravel(), bins=256)
19
ax[1].set_title('Histogram')
20
ax[1].axvline(thresh, color='r')
21
22
ax[2].imshow(binary, cmap=plt.cm.gray)
23
ax[2].set_title('Thresholded')
24
ax[2].axis('off')
25
26
plt.show()

png

Local Thresholding

If the image background is relatively uniform, then you can use a global threshold value as presented above. However, if there is large variation in the background intensity, adaptive thresholding (a.k.a. local or dynamic thresholding) may produce better results. Note that local is much slower than global thresholding.

Here, we binarize an image using the threshold_local function, which calculates thresholds in regions with a characteristic size block_size surrounding each pixel (i.e. local neighborhoods). Each threshold value is the weighted mean of the local neighborhood minus an offset value.


xxxxxxxxxx
29
1
from skimage.filters import threshold_otsu, threshold_local
2
3
4
image = data.page()
5
6
global_thresh = threshold_otsu(image)
7
binary_global = image > global_thresh
8
9
block_size = 35
10
local_thresh = threshold_local(image, block_size, offset=10)
11
binary_local = image > local_thresh
12
13
fig, axes = plt.subplots(ncols=3, figsize=(16, 16))
14
ax = axes.ravel()
15
plt.gray()
16
17
ax[0].imshow(image)
18
ax[0].set_title('Original')
19
20
ax[1].imshow(binary_global)
21
ax[1].set_title('Global thresholding')
22
23
ax[2].imshow(binary_local)
24
ax[2].set_title('Local thresholding')
25
26
for a in ax:
27
    a.axis('off')
28
29
plt.show()

png

Now, we show how Otsu’s threshold method can be applied locally. For each pixel, an “optimal” threshold is determined by maximizing the variance between two classes of pixels of the local neighborhood defined by a structuring element.

The example compares the local threshold with the global threshold.


xxxxxxxxxx
37
1
from skimage.morphology import disk
2
from skimage.filters import threshold_otsu, rank
3
from skimage.util import img_as_ubyte
4
5
6
img = img_as_ubyte(data.page())
7
8
radius = 15
9
selem = disk(radius)
10
11
local_otsu = rank.otsu(img, selem)
12
threshold_global_otsu = threshold_otsu(img)
13
global_otsu = img >= threshold_global_otsu
14
15
fig, axes = plt.subplots(2, 2, figsize=(10, 10), sharex=True, sharey=True)
16
ax = axes.ravel()
17
plt.tight_layout()
18
19
fig.colorbar(ax[0].imshow(img, cmap=plt.cm.gray),
20
             ax=ax[0], orientation='horizontal')
21
ax[0].set_title('Original')
22
ax[0].axis('off')
23
24
fig.colorbar(ax[1].imshow(local_otsu, cmap=plt.cm.gray),
25
             ax=ax[1], orientation='horizontal')
26
ax[1].set_title('Local Otsu (radius=%d)' % radius)
27
ax[1].axis('off')
28
29
ax[2].imshow(img >= local_otsu, cmap=plt.cm.gray)
30
ax[2].set_title('Original >= Local Otsu' % threshold_global_otsu)
31
ax[2].axis('off')
32
33
ax[3].imshow(global_otsu, cmap=plt.cm.gray)
34
ax[3].set_title('Global Otsu (threshold = %d)' % threshold_global_otsu)
35
ax[3].axis('off')
36
37
plt.show()

png

4. Multi-Otsu Thresholding

The multi-Otsu threshold is a thresholding algorithm that is used to separate the pixels of an input image into several different classes, each one obtained according to the intensity of the gray levels within the image.

Multi-Otsu calculates several thresholds, determined by the number of desired classes. The default number of classes is 3: for obtaining three classes, the algorithm returns two threshold values. They are represented by a red line in the histogram below.


xxxxxxxxxx
38
1
import matplotlib
2
from skimage.filters import threshold_multiotsu
3
4
# Setting the font size for all plots.
5
matplotlib.rcParams['font.size'] = 9
6
7
# The input image.
8
image = data.camera()
9
10
# Applying multi-Otsu threshold for the default value, generating
11
# three classes.
12
thresholds = threshold_multiotsu(image)
13
14
# Using the threshold values, we generate the three regions.
15
regions = np.digitize(image, bins=thresholds)
16
17
fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(16, 5))
18
19
# Plotting the original image.
20
ax[0].imshow(image, cmap='gray')
21
ax[0].set_title('Original')
22
ax[0].axis('off')
23
24
# Plotting the histogram and the two thresholds obtained from
25
# multi-Otsu.
26
ax[1].hist(image.ravel(), bins=255)
27
ax[1].set_title('Histogram')
28
for thresh in thresholds:
29
    ax[1].axvline(thresh, color='r')
30
31
# Plotting the Multi Otsu result.
32
ax[2].imshow(regions, cmap='Accent')
33
ax[2].set_title('Multi-Otsu result')
34
ax[2].axis('off')
35
36
plt.subplots_adjust()
37
38
plt.show()

png

5. Active Contour Model

The active contour model is a method to fit open or closed splines to lines or edges in an image. It works by minimising an energy that is in part defined by the image and part by the spline’s shape: length and smoothness. The minimization is done implicitly in the shape energy and explicitly in the image energy.

In the following two examples the active contour model is used (1) to segment the face of a person from the rest of an image by fitting a closed curve to the edges of the face and (2) to find the darkest curve between two fixed points while obeying smoothness considerations. Typically it is a good idea to smooth images a bit before analyzing, as done in the following examples.

We initialize a circle around the astronaut’s face and use the default boundary condition boundary_condition=’periodic’ to fit a closed curve. The default parameters w_line=0, w_edge=1 will make the curve search towards edges, such as the boundaries of the face.


xxxxxxxxxx
27
1
import numpy as np
2
import matplotlib.pyplot as plt
3
from skimage.color import rgb2gray
4
from skimage import data
5
from skimage.filters import gaussian
6
from skimage.segmentation import active_contour
7
8
9
img = data.astronaut()
10
img = rgb2gray(img)
11
12
s = np.linspace(0, 2*np.pi, 400)
13
r = 100 + 100*np.sin(s)
14
c = 220 + 100*np.cos(s)
15
init = np.array([r, c]).T
16
17
snake = active_contour(gaussian(img, 3),
18
                       init, alpha=0.015, beta=10, gamma=0.001)
19
20
fig, ax = plt.subplots(figsize=(7, 7))
21
ax.imshow(img, cmap=plt.cm.gray)
22
ax.plot(init[:, 1], init[:, 0], '--r', lw=3)
23
ax.plot(snake[:, 1], snake[:, 0], '-b', lw=3)
24
ax.set_xticks([]), ax.set_yticks([])
25
ax.axis([0, img.shape[1], img.shape[0], 0])
26
27
plt.show()

png

Here we initialize a straight line between two points, (5, 136) and (424, 50), and require that the spline has its end points there by giving the boundary condition boundary_condition=’fixed’. We furthermore make the algorithm search for dark lines by giving a negative w_line value.


xxxxxxxxxx
17
1
img = data.text()
2
3
r = np.linspace(136, 50, 100)
4
c = np.linspace(5, 424, 100)
5
init = np.array([r, c]).T
6
7
snake = active_contour(gaussian(img, 1), init, boundary_condition='fixed',
8
                       alpha=0.1, beta=1.0, w_line=-5, w_edge=0, gamma=0.1)
9
10
fig, ax = plt.subplots(figsize=(9, 5))
11
ax.imshow(img, cmap=plt.cm.gray)
12
ax.plot(init[:, 1], init[:, 0], '--r', lw=3)
13
ax.plot(snake[:, 1], snake[:, 0], '-b', lw=3)
14
ax.set_xticks([]), ax.set_yticks([])
15
ax.axis([0, img.shape[1], img.shape[0], 0])
16
17
plt.show()

png

Conclusion

The above Image Segmentation techniques were experimented using Scikit-Image library package to study and compare the underlying concepts.