Build gallery with sphinx-gallery

Modified comments in some gallery examples for compatibility with
sphinx-gallery parsing. Also modified some links in the narrative doc
since image file names have changed.
This commit is contained in:
emmanuelle
2016-05-10 08:00:50 +02:00
parent ef9bcfb778
commit 844858f01b
23 changed files with 688 additions and 945 deletions
+37 -45
View File
@@ -35,11 +35,10 @@ def sobel_each(image):
def sobel_hsv(image):
return filters.sobel(image)
"""
We can use these functions as we would normally use them, but now they work
with both gray-scale and color images. Let's plot the results with a color
image:
"""
######################################################################
# We can use these functions as we would normally use them, but now they work
# with both gray-scale and color images. Let's plot the results with a color
# image:
from skimage import data
from skimage.exposure import rescale_intensity
@@ -63,30 +62,27 @@ ax_hsv.imshow(rescale_intensity(1 - sobel_hsv(image)))
ax_hsv.set_xticks([]), ax_hsv.set_yticks([])
ax_hsv.set_title("Sobel filter computed\n on (V)alue converted image (HSV)")
"""
.. image:: PLOT2RST.current_figure
######################################################################
# Notice that the result for the value-filtered image preserves the color of
# the original image, but channel filtered image combines in a more
# surprising way. In other common cases, smoothing for example, the channel
# filtered image will produce a better result than the value-filtered image.
#
# You can also create your own handler functions for ``adapt_rgb``. To do so,
# just create a function with the following signature::
#
# def handler(image_filter, image, *args, **kwargs):
# # Manipulate RGB image here...
# image = image_filter(image, *args, **kwargs)
# # Manipulate filtered image here...
# return image
#
# Note that ``adapt_rgb`` handlers are written for filters where the image is
# the first argument.
#
# As a very simple example, we can just convert any RGB image to grayscale
# and then return the filtered result:
Notice that the result for the value-filtered image preserves the color of the
original image, but channel filtered image combines in a more surprising way.
In other common cases, smoothing for example, the channel filtered image will
produce a better result than the value-filtered image.
You can also create your own handler functions for ``adapt_rgb``. To do so,
just create a function with the following signature::
def handler(image_filter, image, *args, **kwargs):
# Manipulate RGB image here...
image = image_filter(image, *args, **kwargs)
# Manipulate filtered image here...
return image
Note that ``adapt_rgb`` handlers are written for filters where the image is the
first argument.
As a very simple example, we can just convert any RGB image to grayscale and
then return the filtered result:
"""
from skimage.color import rgb2gray
@@ -94,13 +90,12 @@ def as_gray(image_filter, image, *args, **kwargs):
gray_image = rgb2gray(image)
return image_filter(gray_image, *args, **kwargs)
"""
It's important to create a signature that uses ``*args`` and ``**kwargs`` to
pass arguments along to the filter so that the decorated function is allowed to
have any number of positional and keyword arguments.
Finally, we can use this handler with ``adapt_rgb`` just as before:
"""
######################################################################
# It's important to create a signature that uses ``*args`` and ``**kwargs``
# to pass arguments along to the filter so that the decorated function is
# allowed to have any number of positional and keyword arguments.
#
# Finally, we can use this handler with ``adapt_rgb`` just as before:
@adapt_rgb(as_gray)
@@ -119,13 +114,10 @@ ax.set_title("Sobel filter computed\n on the converted grayscale image")
plt.show()
"""
.. image:: PLOT2RST.current_figure
.. note::
A very simple check of the array shape is used for detecting RGB images, so
``adapt_rgb`` is not recommended for functions that support 3D volumes or
color images in non-RGB spaces.
"""
######################################################################
#
# .. note::
#
# A very simple check of the array shape is used for detecting RGB
# images, so ``adapt_rgb`` is not recommended for functions that support
# 3D volumes or color images in non-RGB spaces.
@@ -48,11 +48,9 @@ for ax in axes.ravel():
fig.tight_layout()
"""
.. image:: PLOT2RST.current_figure
######################################################################
# Now we can easily manipulate the hematoxylin and DAB "channels":
Now we can easily manipulate the hematoxylin and DAB "channels":
"""
import numpy as np
from skimage.exposure import rescale_intensity
@@ -70,6 +68,3 @@ ax.set_title("Stain separated image (rescaled)")
ax.axis('off')
plt.show()
"""
.. image:: PLOT2RST.current_figure
"""
@@ -31,10 +31,9 @@ mask = image
dilated = reconstruction(seed, mask, method='dilation')
"""
Subtracting the dilated image leaves an image with just the coins and a flat,
black background, as shown below.
"""
######################################################################
# Subtracting the dilated image leaves an image with just the coins and a
# flat, black background, as shown below.
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(8, 2.5), sharex=True, sharey=True)
@@ -55,29 +54,25 @@ ax3.set_adjustable('box-forced')
fig.tight_layout()
"""
.. image:: PLOT2RST.current_figure
Although the features (i.e. the coins) are clearly isolated, the coins
surrounded by a bright background in the original image are dimmer in the
subtracted image. We can attempt to correct this using a different seed image.
Instead of creating a seed image with maxima along the image border, we can use
the features of the image itself to seed the reconstruction process. Here, the
seed image is the original image minus a fixed value, ``h``.
"""
######################################################################
# Although the features (i.e. the coins) are clearly isolated, the coins
# surrounded by a bright background in the original image are dimmer in the
# subtracted image. We can attempt to correct this using a different seed
# image.
#
# Instead of creating a seed image with maxima along the image border, we can
# use the features of the image itself to seed the reconstruction process.
# Here, the seed image is the original image minus a fixed value, ``h``.
h = 0.4
seed = image - h
dilated = reconstruction(seed, mask, method='dilation')
hdome = image - dilated
"""
To get a feel for the reconstruction process, we plot the intensity of the
mask, seed, and dilated images along a slice of the image (indicated by red
line).
"""
######################################################################
# To get a feel for the reconstruction process, we plot the intensity of the
# mask, seed, and dilated images along a slice of the image (indicated by red
# line).
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(8, 2.5))
@@ -104,14 +99,11 @@ ax3.axis('off')
fig.tight_layout()
plt.show()
"""
.. image:: PLOT2RST.current_figure
As you can see in the image slice, each coin is given a different baseline
intensity in the reconstructed image; this is because we used the local
intensity (shifted by ``h``) as a seed value. As a result, the coins in the
subtracted image have similar pixel intensities. The final result is known as
the h-dome of an image since this tends to isolate regional maxima of height
``h``. This operation is particularly useful when your images are unevenly
illuminated.
"""
######################################################################
# As you can see in the image slice, each coin is given a different baseline
# intensity in the reconstructed image; this is because we used the local
# intensity (shifted by ``h``) as a seed value. As a result, the coins in the
# subtracted image have similar pixel intensities. The final result is known
# as the h-dome of an image since this tends to isolate regional maxima of
# height ``h``. This operation is particularly useful when your images are
# unevenly illuminated.
@@ -34,26 +34,23 @@ ax2.imshow(yellow_multiplier * image)
ax1.set_adjustable('box-forced')
ax2.set_adjustable('box-forced')
"""
.. image:: PLOT2RST.current_figure
In many cases, dealing with RGB values may not be ideal. Because of that, there
are many other `color spaces`_ in which you can represent a color image. One
popular color space is called HSV, which represents hue (~the color),
saturation (~colorfulness), and value (~brightness). For example, a color
(hue) might be green, but its saturation is how intense that green is---where
olive is on the low end and neon on the high end.
In some implementations, the hue in HSV goes from 0 to 360, since hues wrap
around in a circle. In scikit-image, however, hues are float values from 0 to
1, so that hue, saturation, and value all share the same scale.
.. _color spaces:
http://en.wikipedia.org/wiki/List_of_color_spaces_and_their_uses
Below, we plot a linear gradient in the hue, with the saturation and value
turned all the way up:
"""
######################################################################
# In many cases, dealing with RGB values may not be ideal. Because of that,
# there are many other `color spaces`_ in which you can represent a color
# image. One popular color space is called HSV, which represents hue (~the
# color), saturation (~colorfulness), and value (~brightness). For example, a
# color (hue) might be green, but its saturation is how intense that green is
# ---where olive is on the low end and neon on the high end.
#
# In some implementations, the hue in HSV goes from 0 to 360, since hues wrap
# around in a circle. In scikit-image, however, hues are float values from 0
# to 1, so that hue, saturation, and value all share the same scale.
#
# .. _color spaces:
# http://en.wikipedia.org/wiki/List_of_color_spaces_and_their_uses
#
# Below, we plot a linear gradient in the hue, with the saturation and value
# turned all the way up:
import numpy as np
hue_gradient = np.linspace(0, 1)
@@ -67,22 +64,17 @@ fig, ax = plt.subplots(figsize=(5, 2))
ax.imshow(all_hues, extent=(0, 1, 0, 0.2))
ax.set_axis_off()
"""
.. image:: PLOT2RST.current_figure
Notice how the colors at the far left and far right are the same. That reflects
the fact that the hues wrap around like the color wheel (see HSV_ for more
info).
.. _HSV: http://en.wikipedia.org/wiki/HSL_and_HSV
Now, let's create a little utility function to take an RGB image and:
1. Transform the RGB image to HSV
2. Set the hue and saturation
3. Transform the HSV image back to RGB
"""
######################################################################
# Notice how the colors at the far left and far right are the same. That
# reflects the fact that the hues wrap around like the color wheel (see HSV_
# for more info).
#
# .. _HSV: http://en.wikipedia.org/wiki/HSL_and_HSV
#
# Now, let's create a little utility function to take an RGB image and:
#
# 1. Transform the RGB image to HSV 2. Set the hue and saturation 3.
# Transform the HSV image back to RGB
def colorize(image, hue, saturation=1):
@@ -96,13 +88,13 @@ def colorize(image, hue, saturation=1):
return color.hsv2rgb(hsv)
"""
Notice that we need to bump up the saturation; images with zero saturation are
grayscale, so we need to a non-zero value to actually see the color we've set.
Using the function above, we plot six images with a linear gradient in the hue
and a non-zero saturation:
"""
######################################################################
# Notice that we need to bump up the saturation; images with zero saturation
# are grayscale, so we need to a non-zero value to actually see the color
# we've set.
#
# Using the function above, we plot six images with a linear gradient in the
# hue and a non-zero saturation:
hue_rotations = np.linspace(0, 1, 6)
@@ -116,15 +108,12 @@ for ax, hue in zip(axes.flat, hue_rotations):
ax.set_adjustable('box-forced')
fig.tight_layout()
"""
.. image:: PLOT2RST.current_figure
You can combine this tinting effect with numpy slicing and fancy-indexing to
selectively tint your images. In the example below, we set the hue of some
rectangles using slicing and scale the RGB values of some pixels found by
thresholding. In practice, you might want to define a region for tinting based
on segmentation results or blob detection methods.
"""
######################################################################
# You can combine this tinting effect with numpy slicing and fancy-indexing
# to selectively tint your images. In the example below, we set the hue of
# some rectangles using slicing and scale the RGB values of some pixels found
# by thresholding. In practice, you might want to define a region for tinting
# based on segmentation results or blob detection methods.
from skimage.filters import rank
@@ -153,11 +142,7 @@ ax2.set_adjustable('box-forced')
plt.show()
"""
.. image:: PLOT2RST.current_figure
For coloring multiple regions, you may also be interested in
`skimage.color.label2rgb
<http://scikit-image.org/docs/0.9.x/api/skimage.color.html#label2rgb>`_.
"""
######################################################################
# For coloring multiple regions, you may also be interested in
# `skimage.color.label2rgb <http://scikit-
# image.org/docs/0.9.x/api/skimage.color.html#label2rgb>`_.
+5 -12
View File
@@ -64,14 +64,11 @@ if new_scipy:
ax.set_xticks([]), ax.set_yticks([])
ax.axis([0, img.shape[1], img.shape[0], 0])
"""
.. image:: PLOT2RST.current_figure
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 `bc='fixed'`. We furthermore make the algorithm search
for dark lines by giving a negative `w_line` value.
"""
######################################################################
# 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 `bc='fixed'`. We furthermore make the algorithm
# search for dark lines by giving a negative `w_line` value.
img = data.text()
@@ -93,7 +90,3 @@ if new_scipy:
ax.axis([0, img.shape[1], img.shape[0], 0])
plt.show()
"""
.. image:: PLOT2RST.current_figure
"""
+24 -32
View File
@@ -77,33 +77,30 @@ for idx in np.argsort(accums)[::-1][:5]:
ax.imshow(image, cmap=plt.cm.gray)
"""
.. image:: PLOT2RST.current_figure
Ellipse detection
=================
In this second example, the aim is to detect the edge of a coffee cup.
Basically, this is a projection of a circle, i.e. an ellipse.
The problem to solve is much more difficult because five parameters have to be
determined, instead of three for circles.
Algorithm overview
------------------
The algorithm takes two different points belonging to the ellipse. It assumes
that it is the main axis. A loop on all the other points determines how much
an ellipse passes to them. A good match corresponds to high accumulator values.
A full description of the algorithm can be found in reference [1]_.
References
----------
.. [1] Xie, Yonghong, and Qiang Ji. "A new efficient ellipse detection
method." Pattern Recognition, 2002. Proceedings. 16th International
Conference on. Vol. 2. IEEE, 2002
"""
######################################################################
# Ellipse detection
# =================
#
# In this second example, the aim is to detect the edge of a coffee cup.
# Basically, this is a projection of a circle, i.e. an ellipse. The problem
# to solve is much more difficult because five parameters have to be
# determined, instead of three for circles.
#
# Algorithm overview
# -------------------
#
# The algorithm takes two different points belonging to the ellipse. It
# assumes that it is the main axis. A loop on all the other points determines
# how much an ellipse passes to them. A good match corresponds to high
# accumulator values.
#
# A full description of the algorithm can be found in reference [1]_.
#
# References
# ----------
# .. [1] Xie, Yonghong, and Qiang Ji. "A new efficient
# ellipse detection method." Pattern Recognition, 2002. Proceedings.
# 16th International Conference on. Vol. 2. IEEE, 2002
import matplotlib.pyplot as plt
@@ -149,8 +146,3 @@ ax2.set_title('Edge (white) and result (red)')
ax2.imshow(edges)
plt.show()
"""
.. image:: PLOT2RST.current_figure
"""
+20 -26
View File
@@ -31,28 +31,26 @@ ax1.axis('off')
plt.tight_layout()
"""
.. image:: PLOT2RST.current_figure
Different operators compute different finite-difference approximations of the
gradient. For example, the Scharr filter results in a less rotational variance
than the Sobel filter that is in turn better than the Prewitt filter [1]_ [2]_
[3]_. The difference between the Prewitt and Sobel filters and the Scharr filter
is illustrated below with an image that is the discretization of a rotation-
invariant continuous function. The discrepancy between the Prewitt and Sobel
filters, and the Scharr filter is stronger for regions of the image where the
direction of the gradient is close to diagonal, and for regions with high
spatial frequencies. For the example image the differences between the filter
results are very small and the filter results are visually almost
indistinguishable.
.. [1] https://en.wikipedia.org/wiki/Sobel_operator#Alternative_operators
.. [2] B. Jaehne, H. Scharr, and S. Koerkel. Principles of filter design. In
Handbook of Computer Vision and Applications. Academic Press, 1999.
.. [3] https://en.wikipedia.org/wiki/Prewitt_operator
"""
######################################################################
# Different operators compute different finite-difference approximations of
# the gradient. For example, the Scharr filter results in a less rotational
# variance than the Sobel filter that is in turn better than the Prewitt
# filter [1]_ [2]_ [3]_. The difference between the Prewitt and Sobel filters
# and the Scharr filter is illustrated below with an image that is the
# discretization of a rotation- invariant continuous function. The
# discrepancy between the Prewitt and Sobel filters, and the Scharr filter is
# stronger for regions of the image where the direction of the gradient is
# close to diagonal, and for regions with high spatial frequencies. For the
# example image the differences between the filter results are very small and
# the filter results are visually almost indistinguishable.
#
# .. [1] https://en.wikipedia.org/wiki/Sobel_operator#Alternative_operators
#
# .. [2] B. Jaehne, H. Scharr, and S. Koerkel. Principles of filter design.
# In Handbook of Computer Vision and Applications. Academic Press,
# 1999.
#
# .. [3] https://en.wikipedia.org/wiki/Prewitt_operator
x, y = np.ogrid[:100, :100]
# Rotation-invariant image with different spatial frequencies
@@ -86,7 +84,3 @@ ax3.axis('off')
plt.tight_layout()
plt.show()
"""
.. image:: PLOT2RST.current_figure
"""
@@ -18,14 +18,14 @@ Usually, lines are parameterised as :math:`y = mx + c`, with a gradient
infinity for vertical lines. Instead, we therefore construct a segment
perpendicular to the line, leading to the origin. The line is represented by
the length of that segment, :math:`r`, and the angle it makes with the x-axis,
:math:`\theta`.
:math:`\\theta`.
The Hough transform constructs a histogram array representing the parameter
space (i.e., an :math:`M \times N` matrix, for :math:`M` different values of
the radius and :math:`N` different values of :math:`\theta`). For each
parameter combination, :math:`r` and :math:`\theta`, we then find the number of
space (i.e., an :math:`M \\times N` matrix, for :math:`M` different values of
the radius and :math:`N` different values of :math:`\\theta`). For each
parameter combination, :math:`r` and :math:`\\theta`, we then find the number of
non-zero pixels in the input image that would fall close to the corresponding
line, and increment the array at position :math:`(r, \theta)` appropriately.
line, and increment the array at position :math:`(r, \\theta)` appropriately.
We can think of each non-zero pixel "voting" for potential line candidates. The
local maxima in the resulting histogram indicates the parameters of the most
@@ -29,16 +29,13 @@ ax[0].imshow(image)
ax[0].set_title('Original image')
ax[0].axis('off')
"""
.. image:: PLOT2RST.current_figure
Now we need to create the seed image, where the minima represent the starting
points for erosion. To fill holes, we initialize the seed image to the maximum
value of the original image. Along the borders, however, we use the original
values of the image. These border pixels will be the starting points for the
erosion process. We then limit the erosion by setting the mask to the values
of the original image.
"""
######################################################################
# Now we need to create the seed image, where the minima represent the
# starting points for erosion. To fill holes, we initialize the seed image
# to the maximum value of the original image. Along the borders, however, we
# use the original values of the image. These border pixels will be the
# starting points for the erosion process. We then limit the erosion by
# setting the mask to the values of the original image.
import numpy as np
from skimage.morphology import reconstruction
@@ -52,28 +49,24 @@ filled = reconstruction(seed, mask, method='erosion')
ax[1].imshow(filled)
ax[1].set_title('after filling holes')
ax[1].axis('off')
"""
.. image:: PLOT2RST.current_figure
As shown above, eroding inward from the edges removes holes, since (by
definition) holes are surrounded by pixels of brighter value. Finally, we can
isolate the dark regions by subtracting the reconstructed image from the
original image.
"""
######################################################################
# As shown above, eroding inward from the edges removes holes, since (by
# definition) holes are surrounded by pixels of brighter value. Finally, we
# can isolate the dark regions by subtracting the reconstructed image from
# the original image.
ax[2].imshow(image-filled)
ax[2].set_title('holes')
ax[2].axis('off')
"""
.. image:: PLOT2RST.current_figure
Alternatively, we can find bright spots in an image using morphological
reconstruction by dilation. Dilation is the inverse of erosion and expands the
*maximal* values of the seed image until it encounters a mask image. Since this
is an inverse operation, we initialize the seed image to the minimum image
intensity instead of the maximum. The remainder of the process is the same.
"""
######################################################################
# Alternatively, we can find bright spots in an image using morphological
# reconstruction by dilation. Dilation is the inverse of erosion and expands
# the *maximal* values of the seed image until it encounters a mask image.
# Since this is an inverse operation, we initialize the seed image to the
# minimum image intensity instead of the maximum. The remainder of the
# process is the same.
seed = np.copy(image)
seed[1:-1, 1:-1] = image.min()
@@ -83,7 +76,3 @@ ax[3].imshow(image-rec)
ax[3].set_title('peaks')
ax[3].axis('off')
plt.show()
"""
.. image:: PLOT2RST.current_figure
"""
@@ -69,21 +69,18 @@ for ax, values, name in zip(axes, binary_patterns, titles):
plot_lbp_model(ax, values)
ax.set_title(name)
"""
.. image:: PLOT2RST.current_figure
The figure above shows example results with black (or white) representing
pixels that are less (or more) intense than the central pixel. When surrounding
pixels are all black or all white, then that image region is flat (i.e.
featureless). Groups of continuous black or white pixels are considered
"uniform" patterns that can be interpreted as corners or edges. If pixels
switch back-and-forth between black and white pixels, the pattern is considered
"non-uniform".
When using LBP to detect texture, you measure a collection of LBPs over an
image patch and look at the distribution of these LBPs. Lets apply LBP to
a brick texture.
"""
######################################################################
# The figure above shows example results with black (or white) representing
# pixels that are less (or more) intense than the central pixel. When
# surrounding pixels are all black or all white, then that image region is
# flat (i.e. featureless). Groups of continuous black or white pixels are
# considered "uniform" patterns that can be interpreted as corners or edges.
# If pixels switch back-and-forth between black and white pixels, the pattern
# is considered "non-uniform".
#
# When using LBP to detect texture, you measure a collection of LBPs over an
# image patch and look at the distribution of these LBPs. Lets apply LBP to a
# brick texture.
from skimage.transform import rotate
from skimage.feature import local_binary_pattern
@@ -145,16 +142,13 @@ for ax in ax_img:
ax.axis('off')
"""
.. image:: PLOT2RST.current_figure
The above plot highlights flat, edge-like, and corner-like regions of the
image.
The histogram of the LBP result is a good measure to classify textures. Here,
we test the histogram distributions against each other using the
Kullback-Leibler-Divergence.
"""
######################################################################
# The above plot highlights flat, edge-like, and corner-like regions of the
# image.
#
# The histogram of the LBP result is a good measure to classify textures.
# Here, we test the histogram distributions against each other using the
# Kullback-Leibler-Divergence.
# settings for LBP
radius = 2
@@ -222,8 +216,4 @@ ax3.imshow(wall)
ax3.axis('off')
hist(ax6, refs['wall'])
"""
.. image:: PLOT2RST.current_figure
"""
plt.show()
@@ -43,10 +43,10 @@ lbp_code = multiblock_lbp(int_img, 0, 0, 3, 3)
assert_equal(correct_answer, lbp_code)
"""
Now let's apply the operator to a real image and see how the
visualization works.
"""
######################################################################
# Now let's apply the operator to a real image and see how the visualization
# works.
from skimage import data
from matplotlib import pyplot as plt
from skimage.feature import draw_multiblock_lbp
@@ -65,11 +65,9 @@ plt.imshow(img, interpolation='nearest')
plt.show()
"""
.. image:: PLOT2RST.current_figure
On the above plot we see the result of computing a MB-LBP and visualization of
the computed feature. The rectangles that have less intensities' sum than the
central rectangle are marked in cyan. The ones that have higher intensity
values are marked in white. The central rectangle is left untouched.
"""
######################################################################
# On the above plot we see the result of computing a MB-LBP and visualization
# of the computed feature. The rectangles that have less intensities' sum
# than the central rectangle are marked in cyan. The ones that have higher
# intensity values are marked in white. The central rectangle is left
# untouched.
+58 -69
View File
@@ -76,23 +76,20 @@ ax2.imshow(sinogram, cmap=plt.cm.Greys_r,
fig.tight_layout()
plt.show()
"""
.. image:: PLOT2RST.current_figure
Reconstruction with the Filtered Back Projection (FBP)
======================================================
The mathematical foundation of the filtered back projection is the Fourier
slice theorem [2]_. It uses Fourier transform of the projection and
interpolation in Fourier space to obtain the 2D Fourier transform of the image,
which is then inverted to form the reconstructed image. The filtered back
projection is among the fastest methods of performing the inverse Radon
transform. The only tunable parameter for the FBP is the filter, which is
applied to the Fourier transformed projections. It may be used to suppress
high frequency noise in the reconstruction. ``skimage`` provides a few
different options for the filter.
"""
######################################################################
#
# Reconstruction with the Filtered Back Projection (FBP)
# ======================================================
#
# The mathematical foundation of the filtered back projection is the Fourier
# slice theorem [2]_. It uses Fourier transform of the projection and
# interpolation in Fourier space to obtain the 2D Fourier transform of the
# image, which is then inverted to form the reconstructed image. The filtered
# back projection is among the fastest methods of performing the inverse
# Radon transform. The only tunable parameter for the FBP is the filter,
# which is applied to the Fourier transformed projections. It may be used to
# suppress high frequency noise in the reconstruction. ``skimage`` provides a
# few different options for the filter.
from skimage.transform import iradon
@@ -110,42 +107,39 @@ ax2.set_title("Reconstruction error\nFiltered back projection")
ax2.imshow(reconstruction_fbp - image, cmap=plt.cm.Greys_r, **imkwargs)
plt.show()
"""
.. image:: PLOT2RST.current_figure
Reconstruction with the Simultaneous Algebraic Reconstruction Technique
=======================================================================
Algebraic reconstruction techniques for tomography are based on a
straightforward idea: for a pixelated image the value of a single ray in a
particular projection is simply a sum of all the pixels the ray passes through
on its way through the object. This is a way of expressing the forward Radon
transform. The inverse Radon transform can then be formulated as a (large) set
of linear equations. As each ray passes through a small fraction of the pixels
in the image, this set of equations is sparse, allowing iterative solvers for
sparse linear systems to tackle the system of equations. One iterative method
has been particularly popular, namely Kaczmarz' method [3]_, which has the
property that the solution will approach a least-squares solution of the
equation set.
The combination of the formulation of the reconstruction problem as a set
of linear equations and an iterative solver makes algebraic techniques
relatively flexible, hence some forms of prior knowledge can be incorporated
with relative ease.
``skimage`` provides one of the more popular variations of the algebraic
reconstruction techniques: the Simultaneous Algebraic Reconstruction Technique
(SART) [1]_ [4]_. It uses Kaczmarz' method [3]_ as the iterative solver. A good
reconstruction is normally obtained in a single iteration, making the method
computationally effective. Running one or more extra iterations will normally
improve the reconstruction of sharp, high frequency features and reduce the
mean squared error at the expense of increased high frequency noise (the user
will need to decide on what number of iterations is best suited to the problem
at hand. The implementation in ``skimage`` allows prior information of the
form of a lower and upper threshold on the reconstructed values to be supplied
to the reconstruction.
"""
######################################################################
#
# Reconstruction with the Simultaneous Algebraic Reconstruction Technique
# =======================================================================
#
# Algebraic reconstruction techniques for tomography are based on a
# straightforward idea: for a pixelated image the value of a single ray in a
# particular projection is simply a sum of all the pixels the ray passes
# through on its way through the object. This is a way of expressing the
# forward Radon transform. The inverse Radon transform can then be formulated
# as a (large) set of linear equations. As each ray passes through a small
# fraction of the pixels in the image, this set of equations is sparse,
# allowing iterative solvers for sparse linear systems to tackle the system
# of equations. One iterative method has been particularly popular, namely
# Kaczmarz' method [3]_, which has the property that the solution will
# approach a least-squares solution of the equation set.
#
# The combination of the formulation of the reconstruction problem as a set
# of linear equations and an iterative solver makes algebraic techniques
# relatively flexible, hence some forms of prior knowledge can be
# incorporated with relative ease.
#
# ``skimage`` provides one of the more popular variations of the algebraic
# reconstruction techniques: the Simultaneous Algebraic Reconstruction
# Technique (SART) [1]_ [4]_. It uses Kaczmarz' method [3]_ as the iterative
# solver. A good reconstruction is normally obtained in a single iteration,
# making the method computationally effective. Running one or more extra
# iterations will normally improve the reconstruction of sharp, high
# frequency features and reduce the mean squared error at the expense of
# increased high frequency noise (the user will need to decide on what number
# of iterations is best suited to the problem at hand. The implementation in
# ``skimage`` allows prior information of the form of a lower and upper
# threshold on the reconstructed values to be supplied to the reconstruction.
from skimage.transform import iradon_sart
@@ -176,19 +170,14 @@ ax4.set_title("Reconstruction error\nSART, 2 iterations")
ax4.imshow(reconstruction_sart2 - image, cmap=plt.cm.Greys_r, **imkwargs)
plt.show()
"""
.. image:: PLOT2RST.current_figure
.. [1] AC Kak, M Slaney, "Principles of Computerized Tomographic Imaging",
IEEE Press 1988. http://www.slaney.org/pct/pct-toc.html
.. [2] Wikipedia, Radon transform,
http://en.wikipedia.org/wiki/Radon_transform#Relationship_with_the_Fourier_transform
.. [3] S Kaczmarz, "Angenaeherte Aufloesung von Systemen linearer
Gleichungen", Bulletin International de l'Academie Polonaise des
Sciences et des Lettres 35 pp 355--357 (1937)
.. [4] AH Andersen, AC Kak, "Simultaneous algebraic reconstruction technique
(SART): a superior implementation of the ART algorithm", Ultrasonic
Imaging 6 pp 81--94 (1984)
"""
######################################################################
# References
#
# .. [1] AC Kak, M Slaney, "Principles of Computerized Tomographic Imaging", IEEE Press 1988. http://www.slaney.org/pct/pct-toc.html
#
# .. [2] Wikipedia, Radon transform, http://en.wikipedia.org/wiki/Radon_transform#Relationship_with_the_Fourier_transform
#
# .. [3] S Kaczmarz, "Angenaeherte Aufloesung von Systemen linearer Gleichungen", Bulletin International de l'Academie Polonaise des Sciences et des Lettres, 35 pp 355--357 (1937)
#
# .. [4] AH Andersen, AC Kak, "Simultaneous algebraic reconstruction
# technique (SART): a superior implementation of the ART algorithm", Ultrasonic Imaging 6 pp 81--94 (1984)
+15 -27
View File
@@ -30,39 +30,31 @@ eimg = filters.sobel(color.rgb2gray(img))
plt.title('Original Image')
plt.imshow(img)
"""
.. image:: PLOT2RST.current_figure
"""
######################################################################
resized = transform.resize(img, (img.shape[0], img.shape[1] - 200))
plt.figure()
plt.title('Resized Image')
plt.imshow(resized)
"""
.. image:: PLOT2RST.current_figure
"""
######################################################################
out = transform.seam_carve(img, eimg, 'vertical', 200)
plt.figure()
plt.title('Resized using Seam Carving')
plt.imshow(out)
"""
.. image:: PLOT2RST.current_figure
Resizing distorts the rocket and surrounding objects, whereas seam carving
removes empty spaces and preserves object proportions.
Object Removal
--------------
Seam carving can also be used to remove artifacts from images.
This requires weighting the artifact with low values. Recall lower weights are
preferentially removed in seam carving. The following code masks the rocket's
region with low weights, indicating it should be removed.
"""
######################################################################
# Resizing distorts the rocket and surrounding objects, whereas seam carving
# removes empty spaces and preserves object proportions.
#
# Object Removal
# --------------
#
# Seam carving can also be used to remove artifacts from images. This
# requires weighting the artifact with low values. Recall lower weights are
# preferentially removed in seam carving. The following code masks the
# rocket's region with low weights, indicating it should be removed.
masked_img = img.copy()
@@ -78,9 +70,7 @@ plt.title('Object Marked')
plt.imshow(masked_img)
"""
.. image:: PLOT2RST.current_figure
"""
######################################################################
eimg[rr, cc] -= 1000
@@ -91,6 +81,4 @@ resized = transform.resize(img, out.shape)
plt.imshow(out)
plt.show()
"""
.. image:: PLOT2RST.current_figure
"""
######################################################################
@@ -22,18 +22,15 @@ ax1.axis('off')
ax2.plot(hist[1][:-1], hist[0], lw=2)
ax2.set_title('histogram of grey values')
"""
.. image:: PLOT2RST.current_figure
Thresholding
============
A simple way to segment the coins is to choose a threshold based on the
histogram of grey values. Unfortunately, thresholding this image gives a binary
image that either misses significant parts of the coins or merges parts of the
background with the coins:
"""
######################################################################
#
# Thresholding
# ============
#
# A simple way to segment the coins is to choose a threshold based on the
# histogram of grey values. Unfortunately, thresholding this image gives a
# binary image that either misses significant parts of the coins or merges
# parts of the background with the coins:
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(6, 3), sharex=True, sharey=True)
ax1.imshow(coins > 100, cmap=plt.cm.gray, interpolation='nearest')
@@ -47,17 +44,13 @@ ax2.set_adjustable('box-forced')
margins = dict(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1)
fig.subplots_adjust(**margins)
"""
.. image:: PLOT2RST.current_figure
Edge-based segmentation
=======================
Next, we try to delineate the contours of the coins using edge-based
segmentation. To do this, we first get the edges of features using the Canny
edge-detector.
"""
######################################################################
# Edge-based segmentation
# =======================
#
# Next, we try to delineate the contours of the coins using edge-based
# segmentation. To do this, we first get the edges of features using the
# Canny edge-detector.
from skimage.feature import canny
edges = canny(coins/255.)
@@ -67,11 +60,8 @@ ax.imshow(edges, cmap=plt.cm.gray, interpolation='nearest')
ax.axis('off')
ax.set_title('Canny detector')
"""
.. image:: PLOT2RST.current_figure
These contours are then filled using mathematical morphology.
"""
######################################################################
# These contours are then filled using mathematical morphology.
from scipy import ndimage as ndi
@@ -82,12 +72,9 @@ ax.imshow(fill_coins, cmap=plt.cm.gray, interpolation='nearest')
ax.axis('off')
ax.set_title('Filling the holes')
"""
.. image:: PLOT2RST.current_figure
Small spurious objects are easily removed by setting a minimum size for valid
objects.
"""
######################################################################
# Small spurious objects are easily removed by setting a minimum size for
# valid objects.
from skimage import morphology
coins_cleaned = morphology.remove_small_objects(fill_coins, 21)
@@ -96,20 +83,16 @@ ax.imshow(coins_cleaned, cmap=plt.cm.gray, interpolation='nearest')
ax.axis('off')
ax.set_title('Removing small objects')
"""
.. image:: PLOT2RST.current_figure
However, this method is not very robust, since contours that are not perfectly
closed are not filled correctly, as is the case for one unfilled coin above.
Region-based segmentation
=========================
We therefore try a region-based method using the watershed transform. First, we
find an elevation map using the Sobel gradient of the image.
"""
######################################################################
# However, this method is not very robust, since contours that are not
# perfectly closed are not filled correctly, as is the case for one unfilled
# coin above.
#
#Region-based segmentation
#=========================
#
#We therefore try a region-based method using the watershed transform.
#First, we find an elevation map using the Sobel gradient of the image.
from skimage.filters import sobel
@@ -120,12 +103,9 @@ ax.imshow(elevation_map, cmap=plt.cm.gray, interpolation='nearest')
ax.axis('off')
ax.set_title('elevation_map')
"""
.. image:: PLOT2RST.current_figure
Next we find markers of the background and the coins based on the extreme parts
of the histogram of grey values.
"""
######################################################################
# Next we find markers of the background and the coins based on the extreme
# parts of the histogram of grey values.
markers = np.zeros_like(coins)
markers[coins < 30] = 1
@@ -136,13 +116,9 @@ ax.imshow(markers, cmap=plt.cm.spectral, interpolation='nearest')
ax.axis('off')
ax.set_title('markers')
"""
.. image:: PLOT2RST.current_figure
Finally, we use the watershed transform to fill regions of the elevation map
starting from the markers determined above:
"""
######################################################################
# Finally, we use the watershed transform to fill regions of the elevation
# map starting from the markers determined above:
segmentation = morphology.watershed(elevation_map, markers)
fig, ax = plt.subplots(figsize=(4, 3))
@@ -150,13 +126,9 @@ ax.imshow(segmentation, cmap=plt.cm.gray, interpolation='nearest')
ax.axis('off')
ax.set_title('segmentation')
"""
.. image:: PLOT2RST.current_figure
This last method works even better, and the coins can be segmented and labeled
individually.
"""
######################################################################
# This last method works even better, and the coins can be segmented and
# labeled individually.
from skimage.color import label2rgb
@@ -175,9 +147,3 @@ ax2.set_adjustable('box-forced')
fig.subplots_adjust(**margins)
"""
.. image:: PLOT2RST.current_figure
"""
plt.show()
+42 -52
View File
@@ -18,48 +18,45 @@ from skimage import transform as tf
margins = dict(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1)
"""
Basics
======
Several different geometric transformation types are supported: similarity,
affine, projective and polynomial.
Geometric transformations can either be created using the explicit parameters
(e.g. scale, shear, rotation and translation) or the transformation matrix:
First we create a transformation using explicit parameters:
"""
######################################################################
# Basics
# ======
#
# Several different geometric transformation types are supported: similarity,
# affine, projective and polynomial.
#
# Geometric transformations can either be created using the explicit
# parameters (e.g. scale, shear, rotation and translation) or the
# transformation matrix:
#
# First we create a transformation using explicit parameters:
tform = tf.SimilarityTransform(scale=1, rotation=math.pi / 2,
translation=(0, 1))
print(tform.params)
"""
Alternatively you can define a transformation by the transformation matrix
itself:
"""
######################################################################
# Alternatively you can define a transformation by the transformation matrix
# itself:
matrix = tform.params.copy()
matrix[1, 2] = 2
tform2 = tf.SimilarityTransform(matrix)
"""
These transformation objects can then be used to apply forward and inverse
coordinate transformations between the source and destination coordinate
systems:
"""
######################################################################
# These transformation objects can then be used to apply forward and inverse
# coordinate transformations between the source and destination coordinate
# systems:
coord = [1, 0]
print(tform2(coord))
print(tform2.inverse(tform(coord)))
"""
Image warping
=============
Geometric transformations can also be used to warp images:
"""
######################################################################
# Image warping
# =============
#
# Geometric transformations can also be used to warp images:
text = data.text()
@@ -79,25 +76,24 @@ ax2.axis('off')
ax3.imshow(back_rotated)
ax3.axis('off')
"""
.. image:: PLOT2RST.current_figure
Parameter estimation
====================
In addition to the basic functionality mentioned above you can also estimate the
parameters of a geometric transformation using the least-squares method.
This can amongst other things be used for image registration or rectification,
where you have a set of control points or homologous/corresponding points in two
images.
Let's assume we want to recognize letters on a photograph which was not taken
from the front but at a certain angle. In the simplest case of a plane paper
surface the letters are projectively distorted. Simple matching algorithms would
not be able to match such symbols. One solution to this problem would be to warp
the image so that the distortion is removed and then apply a matching algorithm:
"""
######################################################################
# Parameter estimation
# ====================
#
# In addition to the basic functionality mentioned above you can also
# estimate the parameters of a geometric transformation using the least-
# squares method.
#
# This can amongst other things be used for image registration or
# rectification, where you have a set of control points or
# homologous/corresponding points in two images.
#
# Let's assume we want to recognize letters on a photograph which was not
# taken from the front but at a certain angle. In the simplest case of a
# plane paper surface the letters are projectively distorted. Simple matching
# algorithms would not be able to match such symbols. One solution to this
# problem would be to warp the image so that the distortion is removed and
# then apply a matching algorithm:
text = data.text()
@@ -126,9 +122,3 @@ ax1.plot(dst[:, 0], dst[:, 1], '.r')
ax1.axis('off')
ax2.imshow(warped)
ax2.axis('off')
"""
.. image:: PLOT2RST.current_figure
"""
plt.show()
+143 -171
View File
@@ -12,12 +12,19 @@ neighborhood around a pixel.
In this document we outline the following basic morphological operations:
1. Erosion
2. Dilation
3. Opening
4. Closing
5. White Tophat
6. Black Tophat
7. Skeletonize
8. Convex Hull
@@ -36,11 +43,8 @@ orig_phantom = img_as_ubyte(io.imread(os.path.join(data_dir, "phantom.png"),
fig, ax = plt.subplots()
ax.imshow(orig_phantom, cmap=plt.cm.gray)
"""
.. image:: PLOT2RST.current_figure
Let's also define a convenience function for plotting comparisons:
"""
######################################################################
# Let's also define a convenience function for plotting comparisons:
def plot_comparison(original, filtered, filter_name):
@@ -55,16 +59,15 @@ def plot_comparison(original, filtered, filter_name):
ax2.axis('off')
ax2.set_adjustable('box-forced')
"""
Erosion
=======
Morphological ``erosion`` sets a pixel at (i, j) to the *minimum over all
pixels in the neighborhood centered at (i, j)*. 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.
"""
######################################################################
# Erosion
# =======
#
# Morphological ``erosion`` sets a pixel at (i, j) to the *minimum over all
# pixels in the neighborhood centered at (i, j)*. 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.morphology import erosion, dilation, opening, closing, white_tophat
from skimage.morphology import black_tophat, skeletonize, convex_hull_image
@@ -74,68 +77,58 @@ selem = disk(6)
eroded = erosion(orig_phantom, selem)
plot_comparison(orig_phantom, eroded, 'erosion')
"""
.. image:: PLOT2RST.current_figure
Notice how the white boundary of the image disappears or gets eroded as we
increase the size of the disk. Also notice the increase in size of the two
black ellipses in the center and the disappearance of the 3 light grey
patches in the lower part of the image.
Dilation
========
Morphological ``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.
"""
######################################################################
# Notice how the white boundary of the image disappears or gets eroded as we
# increase the size of the disk. Also notice the increase in size of the two
# black ellipses in the center and the disappearance of the 3 light grey
# patches in the lower part of the image.
#
#Dilation
#========
#
#Morphological ``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 = dilation(orig_phantom, selem)
plot_comparison(orig_phantom, dilated, 'dilation')
"""
.. image:: PLOT2RST.current_figure
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
=======
Morphological ``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.
"""
######################################################################
# 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
#=======
#
#Morphological ``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 = opening(orig_phantom, selem)
plot_comparison(orig_phantom, opened, 'opening')
"""
.. image:: PLOT2RST.current_figure
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
=======
Morphological ``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:
"""
######################################################################
#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
#=======
#
#Morphological ``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
@@ -143,26 +136,23 @@ phantom[10:30, 200:210] = 0
closed = closing(phantom, selem)
plot_comparison(phantom, closed, 'closing')
"""
.. image:: PLOT2RST.current_figure
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
============
The ``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:
"""
######################################################################
# 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
# ============
#
# The ``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
@@ -171,86 +161,70 @@ phantom[100:110, 200:210] = 0
w_tophat = white_tophat(phantom, selem)
plot_comparison(phantom, w_tophat, 'white tophat')
"""
.. image:: PLOT2RST.current_figure
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
============
The ``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*.
"""
######################################################################
# 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
# ============
#
# The ``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 = black_tophat(phantom, selem)
plot_comparison(phantom, b_tophat, 'black tophat')
"""
.. image:: PLOT2RST.current_figure
As you can see, the 10-pixel wide black square is highlighted since it is
smaller than the structuring element.
Duality
-------
As you should have noticed, many of these operations are simply the reverse
of another operation. This duality can be summarized as follows:
1. Erosion <-> Dilation
2. Opening <-> Closing
3. White tophat <-> Black tophat
Skeletonize
===========
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.
"""
######################################################################
#As you can see, the 10-pixel wide black square is highlighted since
#it is smaller than the structuring element.
#
#**Duality**
#
#As you should have noticed, many of these operations are simply the reverse
#of another operation. This duality can be summarized as follows:
#
# 1. Erosion <-> Dilation
#
# 2. Opening <-> Closing
#
# 3. White tophat <-> Black tophat
#
#Skeletonize
#===========
#
#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_grey=True)
sk = skeletonize(horse == 0)
plot_comparison(horse, sk, 'skeletonize')
"""
.. image:: PLOT2RST.current_figure
As the name suggests, this technique is used to thin the image to 1-pixel wide
skeleton by applying thinning successively.
Convex hull
===========
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.
"""
######################################################################
#
# As the name suggests, this technique is used to thin the image to 1-pixel
# wide skeleton by applying thinning successively.
#
# Convex hull
# ===========
#
# 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 = convex_hull_image(horse == 0)
plot_comparison(horse, hull1, 'convex hull')
"""
.. image:: PLOT2RST.current_figure
As the figure illustrates, ``convex_hull_image`` gives the smallest polygon
which covers the white or True completely in the image.
If we add a small grain to the image, we can see how the convex hull adapts to
enclose that grain:
"""
######################################################################
# As the figure illustrates, ``convex_hull_image`` gives the smallest polygon
# which covers the white or True completely in the image.
#
# If we add a small grain to the image, we can see how the convex hull adapts
# to enclose that grain:
import numpy as np
@@ -260,19 +234,17 @@ horse_mask[45:50, 75:80] = 1
hull2 = convex_hull_image(horse_mask)
plot_comparison(horse_mask, hull2, 'convex hull')
"""
.. image:: PLOT2RST.current_figure
Additional Resources
====================
1. `MathWorks tutorial on morphological processing
<http://www.mathworks.com/help/images/morphology-fundamentals-dilation-and-erosion.html>`_
2. `Auckland university's tutorial on Morphological Image Processing
<http://www.cs.auckland.ac.nz/courses/compsci773s1c/lectures/ImageProcessing-html/topic4.htm>`_
3. http://en.wikipedia.org/wiki/Mathematical_morphology
"""
plt.show()
######################################################################
#
# Additional Resources
# ====================
#
# 1. `MathWorks tutorial on morphological processing
# <http://www.mathworks.com/help/images/morphology-fundamentals-dilation-and-
# erosion.html>`_
#
# 2. `Auckland university's tutorial on Morphological Image
# Processing <http://www.cs.auckland.ac.nz/courses/compsci773s1c/lectures
# /ImageProcessing-html/topic4.htm>`_
#
# 3. http://en.wikipedia.org/wiki/Mathematical_morphology
+156 -241
View File
@@ -50,17 +50,14 @@ ax1.axis('off')
ax2.plot(hist[1][:-1], hist[0], lw=2)
ax2.set_title('Histogram of grey values')
"""
.. image:: PLOT2RST.current_figure
Noise removal
=============
Some noise is added to the image, 1% of pixels are randomly set to 255, 1% are
randomly set to 0. The **median** filter is applied to remove the noise.
"""
######################################################################
#
# Noise removal
# =============
#
# Some noise is added to the image, 1% of pixels are randomly set to 255, 1%
# are randomly set to 0. The **median** filter is applied to remove the
# noise.
from skimage.filters.rank import median
from skimage.morphology import disk
@@ -96,23 +93,20 @@ ax4.axis('off')
ax4.set_adjustable('box-forced')
"""
.. image:: PLOT2RST.current_figure
The added noise is efficiently removed, as the image defaults are small (1
pixel wide), a small filter radius is sufficient. As the radius is increasing,
objects with bigger sizes are filtered as well, such as the camera tripod. The
median filter is often used for noise removal because borders are preserved and
e.g. salt and pepper noise typically does not distort the gray-level.
Image smoothing
================
The example hereunder shows how a local **mean** filter smooths the camera man
image.
"""
######################################################################
#
# The added noise is efficiently removed, as the image defaults are small (1
# pixel wide), a small filter radius is sufficient. As the radius is
# increasing, objects with bigger sizes are filtered as well, such as the
# camera tripod. The median filter is often used for noise removal because
# borders are preserved and e.g. salt and pepper noise typically does not
# distort the gray-level.
#
# Image smoothing
# ===============
#
# The example hereunder shows how a local **mean** filter smooths the camera
# man image.
from skimage.filters.rank import mean
@@ -130,21 +124,17 @@ ax2.set_title('Local mean $r=10$')
ax2.axis('off')
ax2.set_adjustable('box-forced')
"""
.. image:: PLOT2RST.current_figure
One may be interested in smoothing an image while preserving important borders
(median filters already achieved this), here we use the **bilateral** filter
that restricts the local neighborhood to pixel having a gray-level similar to
the central one.
.. note::
A different implementation is available for color images in
`skimage.filters.denoise_bilateral`.
"""
######################################################################
#
# One may be interested in smoothing an image while preserving important
# borders (median filters already achieved this), here we use the
# **bilateral** filter that restricts the local neighborhood to pixel having
# a gray-level similar to the central one.
#
# .. note::
#
# A different implementation is available for color images in
# `skimage.filters.denoise_bilateral`.
from skimage.filters.rank import mean_bilateral
@@ -173,27 +163,21 @@ ax4.imshow(bilat[200:350, 350:450], cmap=plt.cm.gray)
ax4.axis('off')
ax4.set_adjustable('box-forced')
"""
.. image:: PLOT2RST.current_figure
One can see that the large continuous part of the image (e.g. sky) is smoothed
whereas other details are preserved.
Contrast enhancement
====================
We compare here how the global histogram equalization is applied locally.
The equalized image [2]_ has a roughly linear cumulative distribution function
for each pixel neighborhood. The local version [3]_ of the histogram
equalization emphasizes every local gray-level variations.
.. [2] http://en.wikipedia.org/wiki/Histogram_equalization
.. [3] http://en.wikipedia.org/wiki/Adaptive_histogram_equalization
"""
######################################################################
# One can see that the large continuous part of the image (e.g. sky) is
# smoothed whereas other details are preserved.
#
# Contrast enhancement
# ====================
#
# We compare here how the global histogram equalization is applied locally.
#
# The equalized image [2]_ has a roughly linear cumulative distribution
# function for each pixel neighborhood. The local version [3]_ of the
# histogram equalization emphasizes every local gray-level variations.
#
# .. [2] http://en.wikipedia.org/wiki/Histogram_equalization
# .. [3] http://en.wikipedia.org/wiki/Adaptive_histogram_equalization
from skimage import exposure
from skimage.filters import rank
@@ -230,18 +214,13 @@ ax5.axis('off')
ax6.plot(loc_hist[1][:-1], loc_hist[0], lw=2)
ax6.set_title('Histogram of gray values')
"""
.. image:: PLOT2RST.current_figure
Another way to maximize the number of gray-levels used for an image is to apply
a local auto-leveling, i.e. the gray-value of a pixel is proportionally
remapped between local minimum and local maximum.
The following example shows how local auto-level enhances the camara man
picture.
"""
######################################################################
# Another way to maximize the number of gray-levels used for an image is to
# apply a local auto-leveling, i.e. the gray-value of a pixel is
# proportionally remapped between local minimum and local maximum.
#
# The following example shows how local auto-level enhances the camara man
# picture.
from skimage.filters.rank import autolevel
@@ -261,19 +240,14 @@ ax2.set_title('Local autolevel')
ax2.axis('off')
ax2.set_adjustable('box-forced')
"""
.. image:: PLOT2RST.current_figure
This filter is very sensitive to local outliers, see the little white spot in
the left part of the sky. This is due to a local maximum which is very high
comparing to the rest of the neighborhood. One can moderate this using the
percentile version of the auto-level filter which uses given percentiles (one
inferior, one superior) in place of local minimum and maximum. The example
below illustrates how the percentile parameters influence the local auto-level
result.
"""
######################################################################
# This filter is very sensitive to local outliers, see the little white spot
# in the left part of the sky. This is due to a local maximum which is very
# high comparing to the rest of the neighborhood. One can moderate this using
# the percentile version of the auto-level filter which uses given
# percentiles (one inferior, one superior) in place of local minimum and
# maximum. The example below illustrates how the percentile parameters
# influence the local auto-level result.
from skimage.filters.rank import autolevel_percentile
@@ -310,15 +284,10 @@ for i in range(0,len(image_list)):
axes_list[i].axis('off')
axes_list[i].set_adjustable('box-forced')
"""
.. image:: PLOT2RST.current_figure
The morphological contrast enhancement filter replaces the central pixel by the
local maximum if the original pixel value is closest to local maximum,
otherwise by the minimum local.
"""
######################################################################
# The morphological contrast enhancement filter replaces the central pixel by
# the local maximum if the original pixel value is closest to local maximum,
# otherwise by the minimum local.
from skimage.filters.rank import enhance_contrast
@@ -347,14 +316,9 @@ ax4.imshow(enh[200:350, 350:450], cmap=plt.cm.gray)
ax4.axis('off')
ax4.set_adjustable('box-forced')
"""
.. image:: PLOT2RST.current_figure
The percentile version of the local morphological contrast enhancement uses
percentile *p0* and *p1* instead of the local minimum and maximum.
"""
######################################################################
# The percentile version of the local morphological contrast enhancement uses
# percentile *p0* and *p1* instead of the local minimum and maximum.
from skimage.filters.rank import enhance_contrast_percentile
@@ -379,29 +343,25 @@ for ax in ax.ravel():
ax.axis('off')
ax.set_adjustable('box-forced')
"""
.. image:: PLOT2RST.current_figure
Image threshold
===============
The Otsu threshold [1]_ method can be applied locally using the local gray-
level distribution. In the example below, 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
`skimage.filters.threshold_otsu`.
.. note::
Local is much slower than global thresholding. A function for global Otsu
thresholding can be found in : `skimage.filters.threshold_otsu`.
.. [4] http://en.wikipedia.org/wiki/Otsu's_method
"""
######################################################################
#
# Image threshold
# ===============
#
# The Otsu threshold [1]_ method can be applied locally using the local gray-
# level distribution. In the example below, 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
# `skimage.filters.threshold_otsu`.
#
# .. note::
#
# Local is much slower than global thresholding. A function for global
# Otsu thresholding can be found in : `skimage.filters.threshold_otsu`.
#
# .. [4] http://en.wikipedia.org/wiki/Otsu's_method
from skimage.filters.rank import otsu
from skimage.filters import threshold_otsu
@@ -438,14 +398,9 @@ for ax in ax.ravel():
ax.axis('off')
ax.set_adjustable('box-forced')
"""
.. image:: PLOT2RST.current_figure
The following example shows how local Otsu thresholding handles a global level
shift applied to a synthetic image.
"""
######################################################################
# The following example shows how local Otsu thresholding handles a global
# level shift applied to a synthetic image.
n = 100
theta = np.linspace(0, 10 * np.pi, n)
@@ -467,25 +422,20 @@ ax2.set_title('Local Otsu ($r=%d$)' % radius)
ax2.axis('off')
ax2.set_adjustable('box-forced')
"""
.. image:: PLOT2RST.current_figure
Image morphology
================
Local maximum and local minimum are the base operators for gray-level
morphology.
.. note::
`skimage.dilate` and `skimage.erode` are equivalent filters (see below for
comparison).
Here is an example of the classical morphological gray-level filters: opening,
closing and morphological gradient.
"""
######################################################################
# Image morphology
# ================
#
# Local maximum and local minimum are the base operators for gray-level
# morphology.
#
# .. note::
#
# `skimage.dilate` and `skimage.erode` are equivalent filters (see below
# for comparison).
#
# Here is an example of the classical morphological gray-level filters:
# opening, closing and morphological gradient.
from skimage.filters.rank import maximum, minimum, gradient
@@ -514,28 +464,24 @@ ax4.set_title('Morphological gradient')
for ax in ax.ravel():
ax.axis('off')
ax.set_adjustable('box-forced')
"""
.. image:: PLOT2RST.current_figure
Feature extraction
===================
Local histograms can be exploited to compute local entropy, which is related to
the local image complexity. Entropy is computed using base 2 logarithm i.e. the
filter returns the minimum number of bits needed to encode local gray-level
distribution.
`skimage.rank.entropy` returns the local entropy on a given structuring
element. The following example shows applies this filter on 8- and 16-bit
images.
.. note::
to better use the available image bit, the function returns 10x entropy for
8-bit images and 1000x entropy for 16-bit images.
"""
######################################################################
#
# Feature extraction
# ===================
#
# Local histograms can be exploited to compute local entropy, which is
# related to the local image complexity. Entropy is computed using base 2
# logarithm i.e. the filter returns the minimum number of bits needed to
# encode local gray-level distribution.
#
# `skimage.rank.entropy` returns the local entropy on a given structuring
# element. The following example shows applies this filter on 8- and 16-bit
# images.
#
# .. note::
#
# to better use the available image bit, the function returns 10x entropy
# for 8-bit images and 1000x entropy for 16-bit images.
from skimage import data
from skimage.filters.rank import entropy
@@ -557,22 +503,18 @@ ax2.set_title('Entropy')
ax2.axis('off')
ax2.set_adjustable('box-forced')
"""
.. image:: PLOT2RST.current_figure
Implementation
==============
The central part of the `skimage.rank` filters is build on a sliding window
that updates the local gray-level histogram. This approach limits the algorithm
complexity to O(n) where n is the number of image pixels. The complexity is
also limited with respect to the structuring element size.
In the following we compare the performance of different implementations
available in `skimage`.
"""
######################################################################
#
# Implementation
# ==============
#
# The central part of the `skimage.rank` filters is build on a sliding window
# that updates the local gray-level histogram. This approach limits the
# algorithm complexity to O(n) where n is the number of image pixels. The
# complexity is also limited with respect to the structuring element size.
#
# In the following we compare the performance of different implementations
# available in `skimage`.
from time import time
@@ -611,16 +553,13 @@ def cm_dil(image, selem):
def ndi_med(image, n):
return percentile_filter(image, 50, size=n * 2 - 1)
"""
Comparison between
* `filters.rank.maximum`
* `morphology.dilate`
on increasing structuring element size:
"""
######################################################################
# Comparison between
#
# * `filters.rank.maximum`
# * `morphology.dilate`
#
# on increasing structuring element size:
a = data.camera()
@@ -641,13 +580,8 @@ ax.set_xlabel('Element radius')
ax.plot(e_range, rec)
ax.legend(['filters.rank.maximum', 'morphology.dilate'])
"""
.. image:: PLOT2RST.current_figure
and increasing image size:
"""
######################################################################
# and increasing image size:
r = 9
elem = disk(r + 1)
@@ -670,18 +604,13 @@ ax.plot(s_range, rec)
ax.legend(['filters.rank.maximum', 'morphology.dilate'])
"""
.. image:: PLOT2RST.current_figure
Comparison between:
* `filters.rank.median`
* `scipy.ndimage.percentile`
on increasing structuring element size:
"""
######################################################################
# Comparison between:
#
# * `filters.rank.median`
# * `scipy.ndimage.percentile`
#
# on increasing structuring element size:
a = data.camera()
@@ -702,12 +631,8 @@ ax.legend(['filters.rank.median', 'scipy.ndimage.percentile'])
ax.set_ylabel('Time (ms)')
ax.set_xlabel('Element radius')
"""
.. image:: PLOT2RST.current_figure
Comparison of outcome of the three methods:
"""
######################################################################
# Comparison of outcome of the three methods:
fig, (ax0, ax1) = plt.subplots(ncols=2, sharex=True, sharey=True)
ax0.set_title('filters.rank.median')
@@ -719,12 +644,8 @@ ax1.imshow(rndi)
ax1.axis('off')
ax1.set_adjustable('box-forced')
"""
.. image:: PLOT2RST.current_figure
and increasing image size:
"""
######################################################################
# and increasing image size:
r = 9
elem = disk(r + 1)
@@ -746,9 +667,3 @@ ax.legend(['filters.rank.median', 'scipy.ndimage.percentile'])
ax.set_ylabel('Time (ms)')
ax.set_xlabel('Image size')
"""
.. image:: PLOT2RST.current_figure
"""
plt.show()
+4 -1
View File
@@ -14,7 +14,7 @@
import sys
import os
import skimage
import sphinx_gallery
# If extensions (or modules to document with autodoc) are in another directory,
# add these directories to sys.path here. If the directory is relative to the
# documentation root, use os.path.abspath to make it absolute, like shown here.
@@ -37,11 +37,14 @@ extensions = ['sphinx.ext.autodoc',
'sphinx_gallery.gen_gallery'
]
autosummary_generate = True
#------------------------------------------------------------------------
# Sphinx-gallery configuration
#------------------------------------------------------------------------
sphinx_gallery_conf = {
'doc_module' : 'skimage',
# path to your examples scripts
'examples_dirs' : '../examples',
# path where to save gallery generated examples
@@ -6,9 +6,9 @@
Release notes
=============
.. include:: ../release/_release_notes_for_docs.txt
.. include:: ../release/_release_notes_for_docs.rst
Installation
============
.. include:: install.txt
.. include:: install.rst
+1 -1
View File
@@ -81,7 +81,7 @@ disk: ::
... (nrows / 2)**2)
>>> camera[outer_disk_mask] = 0
.. image:: ../auto_examples/numpy_operations/images/plot_camera_numpy_1.png
.. image:: ../../_images/sphx_glr_plot_camera_numpy_001.png
:width: 45%
:target: ../auto_examples/numpy_operations/plot_camera_numpy.html
@@ -78,7 +78,7 @@ using an array of labels to encode the regions to be represented with the
same color.
.. image: ../auto_examples/segmentation/images/plot_join_segmentations_1.png
.. image:: ../../_images/sphx_glr_plot_join_segmentations_001.png
:target: ../auto_examples/segmentation/plot_join_segmentations.html
:align: center
:width: 80%
@@ -87,9 +87,9 @@ same color.
.. topic:: Examples:
* :ref:`example_color_exposure_plot_tinting_grayscale_images.py`
* :ref:`example_segmentation_plot_join_segmentations.py`
* :ref:`example_segmentation_plot_rag_mean_color.py`
* :ref:`sphx_glr_auto_examples_color_exposure_plot_tinting_grayscale_images.py`
* :ref:`sphx_glr_auto_examples_segmentation_plot_join_segmentations.py`
* :ref:`sphx_glr_auto_examples_segmentation_plot_rag_mean_color.py`
Contrast and exposure
@@ -157,9 +157,9 @@ details are enhanced in large regions with poor contrast. As a further
refinement, histogram equalization can be performed in subregions of the
image with :func:`equalize_adapthist`, in order to correct for exposure
gradients across the image. See the example
:ref:`example_color_exposure_plot_equalize.py`.
:ref:`sphx_glr_auto_examples_color_exposure_plot_equalize.py`.
.. image:: ../auto_examples/color_exposure/images/plot_equalize_1.png
.. image:: ../../_images/sphx_glr_plot_equalize_001.png
:target: ../auto_examples/color_exposure/plot_equalize.html
:align: center
:width: 90%
@@ -167,6 +167,6 @@ gradients across the image. See the example
.. topic:: Examples:
* :ref:`example_color_exposure_plot_equalize.py`
* :ref:`sphx_glr_auto_examples_color_exposure_plot_equalize.py`
@@ -11,7 +11,7 @@ the coins cannot be done directly from the histogram of grey values,
because the background shares enough grey levels with the coins that a
thresholding segmentation is not sufficient.
.. image:: ../auto_examples/xx_applications/images/plot_coins_segmentation_1.png
.. image:: ../../_images/sphx_glr_plot_coins_segmentation_001.png
:target: ../auto_examples/xx_applications/plot_coins_segmentation.html
:align: center
@@ -26,7 +26,7 @@ Simply thresholding the image leads either to missing significant parts
of the coins, or to merging parts of the background with the
coins. This is due to the inhomogeneous lighting of the image.
.. image:: ../auto_examples/xx_applications/images/plot_coins_segmentation_2.png
.. image:: ../../_images/sphx_glr_plot_coins_segmentation_002.png
:target: ../auto_examples/xx_applications/plot_coins_segmentation.html
:align: center
@@ -53,7 +53,7 @@ boundary of the coins, or inside the coins.
>>> from scipy import ndimage as ndi
>>> fill_coins = ndi.binary_fill_holes(edges)
.. image:: ../auto_examples/xx_applications/images/plot_coins_segmentation_3.png
.. image:: ../../_images/sphx_glr_plot_coins_segmentation_003.png
:target: ../auto_examples/xx_applications/plot_coins_segmentation.html
:align: center
@@ -62,7 +62,7 @@ we fill the inner part of the coins using the
``ndi.binary_fill_holes`` function, which uses mathematical morphology
to fill the holes.
.. image:: ../auto_examples/xx_applications/images/plot_coins_segmentation_4.png
.. image:: ../../_images/sphx_glr_plot_coins_segmentation_004.png
:target: ../auto_examples/xx_applications/plot_coins_segmentation.html
:align: center
@@ -83,7 +83,7 @@ has not been segmented correctly at all. The reason is that the contour
that we got from the Canny detector was not completely closed, therefore
the filling function did not fill the inner part of the coin.
.. image:: ../auto_examples/xx_applications/images/plot_coins_segmentation_5.png
.. image:: ../../_images/sphx_glr_plot_coins_segmentation_005.png
:target: ../auto_examples/xx_applications/plot_coins_segmentation.html
:align: center
@@ -128,7 +128,7 @@ separate the coins from the background.
and here is the corresponding 2-D plot:
.. image:: ../auto_examples/xx_applications/images/plot_coins_segmentation_6.png
.. image:: ../../_images/sphx_glr_plot_coins_segmentation_006.png
:target: ../auto_examples/xx_applications/plot_coins_segmentation.html
:align: center
@@ -139,7 +139,7 @@ extreme parts of the histogram of grey values::
>>> markers[coins < 30] = 1
>>> markers[coins > 150] = 2
.. image:: ../auto_examples/xx_applications/images/plot_coins_segmentation_7.png
.. image:: ../../_images/sphx_glr_plot_coins_segmentation_007.png
:target: ../auto_examples/xx_applications/plot_coins_segmentation.html
:align: center
@@ -148,7 +148,7 @@ Let us now compute the watershed transform::
>>> from skimage.morphology import watershed
>>> segmentation = watershed(elevation_map, markers)
.. image:: ../auto_examples/xx_applications/images/plot_coins_segmentation_8.png
.. image:: ../../_images/sphx_glr_plot_coins_segmentation_008.png
:target: ../auto_examples/xx_applications/plot_coins_segmentation.html
:align: center
@@ -165,7 +165,7 @@ We can now label all the coins one by one using ``ndi.label``::
>>> labeled_coins, _ = ndi.label(segmentation)
.. image:: ../auto_examples/xx_applications/images/plot_coins_segmentation_9.png
.. image:: ../../_images/sphx_glr_plot_coins_segmentation_009.png
:target: ../auto_examples/xx_applications/plot_coins_segmentation.html
:align: center
+1 -1
View File
@@ -36,7 +36,7 @@ class ApiDocWriter(object):
def __init__(self,
package_name,
rst_extension='.txt',
rst_extension='.rst',
package_skip_patterns=None,
module_skip_patterns=None,
):