%matplotlib inline import numpy as np import matplotlib.pyplot as plt
There are many tools and utilities in the package, far too many to cover in a tutorial. This notebook is designed as a road map, to guide you as you explore or search for additional tools for your applications. This is intended as a guide, not an exhaustive list.
Each submodule of scikit-image has its own section, which you can navigate to below in the table of contents.
import skimage.color as color
# Tab complete to see available functions in the color submodule color.rgb2 color.
from skimage import data from skimage import color original = data.astronaut() grayscale = color.rgb2gray(original) # Plot the results fig, axes = plt.subplots(1, 2, figsize=(8, 4)) ax = axes.ravel() ax.imshow(original) ax.set_title("Original") ax.axis('off') ax.imshow(grayscale, cmap='gray') ax.set_title("Grayscale") ax.axis('off') fig.tight_layout() plt.show();
Usually, objects in images have distinct colors (hues) and luminosities, so that these features can be used to separate different areas of the image. In the RGB representation the hue and the luminosity are expressed as a linear combination of the R,G,B channels, whereas they correspond to single channels of the HSV image (the Hue and the Value channels). A simple segmentation of the image can then be effectively performed by a mere thresholding of the HSV channels. See below link for additional details.
We first load the RGB image and extract the Hue and Value channels:
from skimage import data from skimage.color import rgb2hsv rgb_img = data.coffee() hsv_img = rgb2hsv(rgb_img) hue_img = hsv_img[:, :, 0] value_img = hsv_img[:, :, 2] fig, (ax0, ax1, ax2) = plt.subplots(ncols=3, figsize=(8, 2)) ax0.imshow(rgb_img) ax0.set_title("RGB image") ax0.axis('off') ax1.imshow(hue_img, cmap='hsv') ax1.set_title("Hue channel") ax1.axis('off') ax2.imshow(value_img) ax2.set_title("Value channel") ax2.axis('off') fig.tight_layout();
The cup and saucer have a Hue distinct from the remainder of the image, which can be isolated by thresholding
hue_threshold = 0.04 binary_img = hue_img > hue_threshold fig, (ax0, ax1) = plt.subplots(ncols=2, figsize=(8, 3)) ax0.hist(hue_img.ravel(), 512) ax0.set_title("Histogram of the Hue channel with threshold") ax0.axvline(x=hue_threshold, color='r', linestyle='dashed', linewidth=2) ax0.set_xbound(0, 0.12) ax1.imshow(binary_img) ax1.set_title("Hue-thresholded image") ax1.axis('off') fig.tight_layout();
An additional threshold in the value channel can remote most of the shadow
fig, ax0 = plt.subplots(figsize=(4, 3)) value_threshold = 0.10 binary_img = (hue_img > hue_threshold) | (value_img < value_threshold) ax0.imshow(binary_img) ax0.set_title("Hue and value thresholded image") ax0.axis('off') fig.tight_layout() plt.show();
data submodule includes standard test images useful for examples and testing the package. These images are shipped with the package.
There are scientific images, general test images, and a stereoscopic image.
from skimage import data
# Explore with tab completion example_image = data.camera()
fig, ax = plt.subplots(figsize=(6, 6)) ax.imshow(example_image) ax.axis('off');
# Room for experimentation
The majority of functions in this submodule return the coordinates of the specified shape/object in the image, rather than drawing it on the image directly. The coordinates can then be used as a mask to draw on the image, or you pass the image as well as those coordinates into the convenience function
Lines and circles can be drawn with antialiasing (these functions end in the suffix *_aa).
At the current time text is not supported; other libraries including matplotlib have robust support for overlaying text.
from skimage import draw
# Tab complete to see available options draw.
# Room for experimentation
fig, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(10, 6)) img = np.zeros((500, 500, 3), dtype=np.float64) # draw line rr, cc = draw.line(120, 123, 20, 400) img[rr, cc, 0] = 255 # fill polygon poly = np.array(( (300, 300), (480, 320), (380, 430), (220, 590), (300, 300), )) rr, cc = draw.polygon(poly[:, 0], poly[:, 1], img.shape) img[rr, cc, 1] = 1 # fill circle rr, cc = draw.circle(200, 200, 100, img.shape) img[rr, cc, :] = (1, 1, 0) # fill ellipse rr, cc = draw.ellipse(300, 300, 100, 200, img.shape) img[rr, cc, 2] = 1 # circle rr, cc = draw.circle_perimeter(120, 400, 15) img[rr, cc, :] = (1, 0, 0) # Bezier curve rr, cc = draw.bezier_curve(70, 100, 10, 10, 150, 100, 1) img[rr, cc, :] = (1, 0, 0) # ellipses rr, cc = draw.ellipse_perimeter(120, 400, 60, 20, orientation=np.pi / 4.) img[rr, cc, :] = (1, 0, 1) rr, cc = draw.ellipse_perimeter(120, 400, 60, 20, orientation=-np.pi / 4.) img[rr, cc, :] = (0, 0, 1) rr, cc = draw.ellipse_perimeter(120, 400, 60, 20, orientation=np.pi / 2.) img[rr, cc, :] = (1, 1, 1) ax1.imshow(img) ax1.set_title('No anti-aliasing') ax1.axis('off') img = np.zeros((100, 100), dtype=np.double) # anti-aliased line rr, cc, val = draw.line_aa(12, 12, 20, 50) img[rr, cc] = val # anti-aliased circle rr, cc, val = draw.circle_perimeter_aa(60, 40, 30) img[rr, cc] = val ax2.imshow(img, cmap=plt.cm.gray, interpolation='nearest') ax2.set_title('Anti-aliasing') ax2.axis('off');
One of the most common tools to evaluate exposure is the histogram, which plots the number of points which have a certain value against the values in order from lowest (dark) to highest (light). The function
exposure.histogram differs from
numpy.histogram in that there is no rebinnning; each value along the x-axis is preserved.
from skimage import data, img_as_float from skimage import exposure def plot_img_and_hist(image, axes, bins=256): """Plot an image along with its histogram and cumulative histogram. """ image = img_as_float(image) ax_img, ax_hist = axes ax_cdf = ax_hist.twinx() # Display image ax_img.imshow(image, cmap=plt.cm.gray) ax_img.set_axis_off() # Display histogram ax_hist.hist(image.ravel(), bins=bins, histtype='step', color='black') ax_hist.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0)) ax_hist.set_xlabel('Pixel intensity') ax_hist.set_xlim(0, 1) ax_hist.set_yticks() # Display cumulative distribution img_cdf, bins = exposure.cumulative_distribution(image, bins) ax_cdf.plot(bins, img_cdf, 'r') ax_cdf.set_yticks() return ax_img, ax_hist, ax_cdf # Load an example image img = data.moon() # Contrast stretching p2, p98 = np.percentile(img, (2, 98)) img_rescale = exposure.rescale_intensity(img, in_range=(p2, p98)) # Equalization img_eq = exposure.equalize_hist(img) # Adaptive Equalization img_adapteq = exposure.equalize_adapthist(img, clip_limit=0.03) # Display results fig = plt.figure(figsize=(8, 5)) axes = np.zeros((2, 4), dtype=np.object) axes[0, 0] = fig.add_subplot(2, 4, 1) for i in range(1, 4): axes[0, i] = fig.add_subplot(2, 4, 1+i, sharex=axes[0,0], sharey=axes[0,0]) for i in range(0, 4): axes[1, i] = fig.add_subplot(2, 4, 5+i) ax_img, ax_hist, ax_cdf = plot_img_and_hist(img, axes[:, 0]) ax_img.set_title('Low contrast image') y_min, y_max = ax_hist.get_ylim() ax_hist.set_ylabel('Number of pixels') ax_hist.set_yticks(np.linspace(0, y_max, 5)) ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_rescale, axes[:, 1]) ax_img.set_title('Contrast stretch') ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_eq, axes[:, 2]) ax_img.set_title('Histogram eq') ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_adapteq, axes[:, 3]) ax_img.set_title('Adaptive eq') ax_cdf.set_ylabel('Fraction of total intensity') ax_cdf.set_yticks(np.linspace(0, 1, 5)) # prevent overlap of y-axis labels fig.tight_layout();
# Explore with tab completion exposure.
# Room for experimentation
This submodule presents a diverse set of tools to identify or extract certain features from images, including tools for
from skimage import feature
# Explore with tab completion feature.
# Room for experimentation
This is a large submodule. For brevity here is a short example illustrating ORB feature matching, and additional examples can be explored in the online gallery.
from skimage import data from skimage import transform as tf from skimage import feature from skimage.color import rgb2gray # Import the astronaut then warp/rotate the image img1 = rgb2gray(data.astronaut()) img2 = tf.rotate(img1, 180) tform = tf.AffineTransform(scale=(1.3, 1.1), rotation=0.5, translation=(0, -200)) img3 = tf.warp(img1, tform) # Build ORB extractor and extract features descriptor_extractor = feature.ORB(n_keypoints=200) descriptor_extractor.detect_and_extract(img1) keypoints1 = descriptor_extractor.keypoints descriptors1 = descriptor_extractor.descriptors descriptor_extractor.detect_and_extract(img2) keypoints2 = descriptor_extractor.keypoints descriptors2 = descriptor_extractor.descriptors descriptor_extractor.detect_and_extract(img3) keypoints3 = descriptor_extractor.keypoints descriptors3 = descriptor_extractor.descriptors # Find matches between the extracted features matches12 = feature.match_descriptors(descriptors1, descriptors2, cross_check=True) matches13 = feature.match_descriptors(descriptors1, descriptors3, cross_check=True) # Plot the results fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(10, 10)) plt.gray() feature.plot_matches(ax, img1, img2, keypoints1, keypoints2, matches12) ax.axis('off') ax.set_title("Original Image vs. Transformed Image") feature.plot_matches(ax, img1, img3, keypoints1, keypoints3, matches13) ax.axis('off') ax.set_title("Original Image vs. Transformed Image");
# Room for experimentation
Filtering applies whole-image modifications such as sharpening or blurring. Thresholding methods also live in this submodule.
Notable functions include (links to relevant gallery examples)
from skimage import filters
# Explore with tab completion filters.
Graph theory. Currently this submodule primarily deals with a constructed "cost" image, and how to find the minimum cost path through it, with constraints if desired.
Reading your image and writing the results back out. There are multiple plugins available, which support multiple formats. The most commonly used functions include
Multiple algorithms to label images, or obtain information about discrete regions of an image.
label), find various properties of the labeled regions.
regionpropsis extremely useful
RANDom Sample Consensus fitting (RANSAC) - a powerful, robust approach to fitting a model to data. It exists here because its initial use was for fitting shapes, but it can also fit transforms.
from skimage import measure
# Explore with tab completion measure.
# Room to explore
Morphological image processing is a collection of non-linear operations related to the shape or morphology of features in an image, such as boundaries, skeletons, etc. In any given technique, we probe an image with a small shape or template called a structuring element, which defines the region of interest or neighborhood around a pixel.
from skimage import morphology as morph
# Explore with tab completion morph.
Flood fill is an algorithm to iteratively identify and/or change adjacent values in an image based on their similarity to an initial seed point. The conceptual analogy is the ‘paint bucket’ tool in many graphic editors.
flood function returns the binary mask of the flooded area.
flood_fill returns a modified image. Both of these can be set with a
tolerance keyword argument, within which the adjacent region will be filled.
Here we will experiment a bit on the cameraman, turning his coat from dark to light.
from skimage import data from skimage import morphology as morph cameraman = data.camera() # Change the cameraman's coat from dark to light (255). The seed point is # chosen as (200, 100), light_coat = morph.flood_fill(cameraman, (200, 100), 255, tolerance=10) fig, ax = plt.subplots(ncols=2, figsize=(10, 5)) ax.imshow(cameraman, cmap=plt.cm.gray) ax.set_title('Original') ax.axis('off') ax.imshow(light_coat, cmap=plt.cm.gray) ax.plot(100, 200, 'ro') # seed point ax.set_title('After flood fill') ax.axis('off');
Here we outline the following basic morphological operations:
To get started, let’s load an image using
io.imread. Note that morphology functions only work on gray-scale or binary images, so we set
import os from skimage.data import data_dir from skimage.util import img_as_ubyte from skimage import io orig_phantom = img_as_ubyte(io.imread(os.path.join(data_dir, "phantom.png"), as_gray=True)) fig, ax = plt.subplots(figsize=(5, 5)) ax.imshow(orig_phantom, cmap=plt.cm.gray) ax.axis('off');
def plot_comparison(original, filtered, filter_name): fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(8, 4), sharex=True, sharey=True) ax1.imshow(original, cmap=plt.cm.gray) ax1.set_title('original') ax1.axis('off') ax2.imshow(filtered, cmap=plt.cm.gray) ax2.set_title(filter_name) ax2.axis('off')
erosion sets a pixel at (i, j) to the minimum over all pixels in the neighborhood centered at (i, j). Erosion shrinks bright regions and enlarges dark regions.
The structuring element,
selem, passed to erosion is a boolean array that describes this neighborhood. Below, we use
disk to create a circular structuring element, which we use for most of the following examples.
from skimage import morphology as morph selem = morph.disk(6) eroded = morph.erosion(orig_phantom, selem) plot_comparison(orig_phantom, eroded, 'erosion')
dilation sets a pixel at (i, j) to the maximum over all pixels in the neighborhood centered at (i, j). Dilation enlarges bright regions and shrinks dark regions.
dilated = morph.dilation(orig_phantom, selem) plot_comparison(orig_phantom, dilated, 'dilation')
Notice how the white boundary of the image thickens, or gets dilated, as we increase the size of the disk. Also notice the decrease in size of the two black ellipses in the centre, and the thickening of the light grey circle in the center and the 3 patches in the lower part of the image.
opening on an image is defined as an erosion followed by a dilation. Opening can remove small bright spots (i.e. “salt”) and connect small dark cracks.
opened = morph.opening(orig_phantom, selem) plot_comparison(orig_phantom, opened, 'opening')
Since opening an image starts with an erosion operation, light regions that are smaller than the structuring element are removed. The dilation operation that follows ensures that light regions that are larger than the structuring element retain their original size. Notice how the light and dark shapes in the center their original thickness but the 3 lighter patches in the bottom get completely eroded. The size dependence is highlighted by the outer white ring: The parts of the ring thinner than the structuring element were completely erased, while the thicker region at the top retains its original thickness.
closing on an image is defined as a dilation followed by an erosion. Closing can remove small dark spots (i.e. “pepper”) and connect small bright cracks.
To illustrate this more clearly, let’s add a small crack to the white border:
phantom = orig_phantom.copy() phantom[10:30, 200:210] = 0 closed = morph.closing(phantom, selem) plot_comparison(phantom, closed, 'closing')
Since closing an image starts with an dilation operation, dark regions that are smaller than the structuring element are removed. The dilation operation that follows ensures that dark regions that are larger than the structuring element retain their original size. Notice how the white ellipses at the bottom get connected because of dilation, but other dark region retain their original sizes. Also notice how the crack we added is mostly removed.
white_tophat of an image is defined as the image minus its morphological opening. This operation returns the bright spots of the image that are smaller than the structuring element.
To make things interesting, we’ll add bright and dark spots to the image:
phantom = orig_phantom.copy() phantom[340:350, 200:210] = 255 phantom[100:110, 200:210] = 0 w_tophat = morph.white_tophat(phantom, selem) plot_comparison(phantom, w_tophat, 'white tophat')
As you can see, the 10-pixel wide white square is highlighted since it is smaller than the structuring element. Also, the thin, white edges around most of the ellipse are retained because they’re smaller than the structuring element, but the thicker region at the top disappears.
black_tophat of an image is defined as its morphological closing minus the original image. This operation returns the dark spots of the image that are smaller than the structuring element.
b_tophat = morph.black_tophat(phantom, selem) plot_comparison(phantom, b_tophat, 'black tophat')
As you can see, the 10-pixel wide black square is highlighted since it is smaller than the structuring element.
As you should have noticed, many of these operations are simply the reverse of another operation. This duality can be summarized as follows:
Thinning is used to reduce each connected component in a binary image to a single-pixel wide skeleton. It is important to note that this is performed on binary images only.
horse = io.imread(os.path.join(data_dir, "horse.png"), as_gray=True) sk = morph.skeletonize(horse == 0) plot_comparison(horse, sk, 'skeletonize')
As the name suggests, this technique is used to thin the image to 1-pixel wide skeleton by applying thinning successively.
The convex_hull_image is the set of pixels included in the smallest convex polygon that surround all white pixels in the input image. Again note that this is also performed on binary images.
hull1 = morph.convex_hull_image(horse == 0) plot_comparison(horse, hull1, 'convex hull')
This submodule includes routines to restore images. Currently these routines fall into four major categories. Links lead to topical gallery examples.
from skimage import restoration
# Explore with tab completion restoration.
# Space to experiment with restoration techniques
One of the key image analysis tasks is identifying regions of interest. These could be a person, an object, certain features of an animal, microscopic image, or stars. Segmenting an image is the process of determining where these things you want are in your images.
Segmentation has two overarching categories: Supervised and Unsupervised.
Supervised - must provide some guidance (seed points or initial conditions)
Unsupervised - no human input
from skimage import segmentation
# Explore with tab completion segmentation.
# Room for experimentation
This submodule has multiple features which fall under the umbrella of transformations.
radon) and inverse (
iradon) radon transforms, as well as some variants (
iradon_sart) and the finite versions of these transforms (
ifrt2). These are used for reconstructing medical computed tomography (CT) images.
Hough transforms for identifying lines, circles, and ellipses.
Changing image size, shape, or resolution with
warp_coordinates which take an image or set of coordinates and translate them through one of the defined
*Transforms in this submodule.
estimate_transform may be assist in estimating the parameters.
from skimage import transform
# Explore with tab completion transform.
# Room for experimentation
These are generally useful functions which have no definite other place in the package.
util.img_as_* are convenience functions for datatype conversion.
util.invert is a convenient way to invert any image, accounting for its datatype.
util.random_noise is a comprehensive function to apply any amount of many different types of noise to images. The seed may be set, resulting in pseudo-random noise for testing.
util.view_as_* allows for overlapping views into the same memory array, which is useful for elegant local computations with minimal memory impact.
util.apply_parallel uses Dask to apply a function across subsections of an image. This can result in dramatic performance or memory improvements, but depending on the algorithm edge effects or lack of knowledge of the remainder of the image may result in unexpected results.
util.crop pads or crops the edges of images.
util.pad is now a direct wrapper for
from skimage import util
# Explore with tab completion util.
# Room to experiment