mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-19 11:27:45 +08:00
Added sections to gallery of examples
Modified travis_script.sh to account for the new structure of the gallery Added README.txt files in directories of gallery examples Fixed references to gallery images in user guide pages Fixed broken links
This commit is contained in:
emmanuelle
@@ -0,0 +1,2 @@
|
||||
Geometrical transformations and registration
|
||||
--------------------------------------------
|
||||
@@ -0,0 +1,37 @@
|
||||
"""
|
||||
=========================
|
||||
Interpolation: Edge Modes
|
||||
=========================
|
||||
|
||||
This example illustrates the different edge modes available during
|
||||
interpolation in routines such as `skimage.transform.rescale` and
|
||||
`skimage.transform.resize`.
|
||||
"""
|
||||
from skimage._shared.interpolation import extend_image
|
||||
import skimage.data
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
|
||||
img = np.zeros((16, 16))
|
||||
img[:8, :8] += 1
|
||||
img[:4, :4] += 1
|
||||
img[:2, :2] += 1
|
||||
img[:1, :1] += 2
|
||||
img[8, 8] = 4
|
||||
|
||||
modes = ['constant', 'edge', 'wrap', 'reflect', 'symmetric']
|
||||
fig, axes = plt.subplots(1, 5, figsize=(15, 5))
|
||||
for n, mode in enumerate(modes):
|
||||
img_extended = extend_image(img, pad=img.shape[0], mode=mode)
|
||||
axes[n].imshow(img_extended, cmap=plt.cm.gray, interpolation='nearest')
|
||||
axes[n].plot([15.5, 15.5], [15.5, 31.5], 'y--', linewidth=0.5)
|
||||
axes[n].plot([31.5, 31.5], [15.5, 31.5], 'y--', linewidth=0.5)
|
||||
axes[n].plot([15.5, 31.5], [15.5, 15.5], 'y--', linewidth=0.5)
|
||||
axes[n].plot([15.5, 31.5], [31.5, 31.5], 'y--', linewidth=0.5)
|
||||
axes[n].set_axis_off()
|
||||
axes[n].set_aspect('equal')
|
||||
axes[n].set_title(mode)
|
||||
|
||||
plt.tight_layout()
|
||||
|
||||
plt.show()
|
||||
@@ -0,0 +1,139 @@
|
||||
"""
|
||||
============================
|
||||
Robust matching using RANSAC
|
||||
============================
|
||||
|
||||
In this simplified example we first generate two synthetic images as if they
|
||||
were taken from different view points.
|
||||
|
||||
In the next step we find interest points in both images and find
|
||||
correspondences based on a weighted sum of squared differences of a small
|
||||
neighborhood around them. Note, that this measure is only robust towards
|
||||
linear radiometric and not geometric distortions and is thus only usable with
|
||||
slight view point changes.
|
||||
|
||||
After finding the correspondences we end up having a set of source and
|
||||
destination coordinates which can be used to estimate the geometric
|
||||
transformation between both images. However, many of the correspondences are
|
||||
faulty and simply estimating the parameter set with all coordinates is not
|
||||
sufficient. Therefore, the RANSAC algorithm is used on top of the normal model
|
||||
to robustly estimate the parameter set by detecting outliers.
|
||||
|
||||
"""
|
||||
from __future__ import print_function
|
||||
|
||||
import numpy as np
|
||||
from matplotlib import pyplot as plt
|
||||
|
||||
from skimage import data
|
||||
from skimage.util import img_as_float
|
||||
from skimage.feature import (corner_harris, corner_subpix, corner_peaks,
|
||||
plot_matches)
|
||||
from skimage.transform import warp, AffineTransform
|
||||
from skimage.exposure import rescale_intensity
|
||||
from skimage.color import rgb2gray
|
||||
from skimage.measure import ransac
|
||||
|
||||
|
||||
# generate synthetic checkerboard image and add gradient for the later matching
|
||||
checkerboard = img_as_float(data.checkerboard())
|
||||
img_orig = np.zeros(list(checkerboard.shape) + [3])
|
||||
img_orig[..., 0] = checkerboard
|
||||
gradient_r, gradient_c = (np.mgrid[0:img_orig.shape[0],
|
||||
0:img_orig.shape[1]]
|
||||
/ float(img_orig.shape[0]))
|
||||
img_orig[..., 1] = gradient_r
|
||||
img_orig[..., 2] = gradient_c
|
||||
img_orig = rescale_intensity(img_orig)
|
||||
img_orig_gray = rgb2gray(img_orig)
|
||||
|
||||
# warp synthetic image
|
||||
tform = AffineTransform(scale=(0.9, 0.9), rotation=0.2, translation=(20, -10))
|
||||
img_warped = warp(img_orig, tform.inverse, output_shape=(200, 200))
|
||||
img_warped_gray = rgb2gray(img_warped)
|
||||
|
||||
# extract corners using Harris' corner measure
|
||||
coords_orig = corner_peaks(corner_harris(img_orig_gray), threshold_rel=0.001,
|
||||
min_distance=5)
|
||||
coords_warped = corner_peaks(corner_harris(img_warped_gray),
|
||||
threshold_rel=0.001, min_distance=5)
|
||||
|
||||
# determine sub-pixel corner position
|
||||
coords_orig_subpix = corner_subpix(img_orig_gray, coords_orig, window_size=9)
|
||||
coords_warped_subpix = corner_subpix(img_warped_gray, coords_warped,
|
||||
window_size=9)
|
||||
|
||||
|
||||
def gaussian_weights(window_ext, sigma=1):
|
||||
y, x = np.mgrid[-window_ext:window_ext+1, -window_ext:window_ext+1]
|
||||
g = np.zeros(y.shape, dtype=np.double)
|
||||
g[:] = np.exp(-0.5 * (x**2 / sigma**2 + y**2 / sigma**2))
|
||||
g /= 2 * np.pi * sigma * sigma
|
||||
return g
|
||||
|
||||
|
||||
def match_corner(coord, window_ext=5):
|
||||
r, c = np.round(coord).astype(np.intp)
|
||||
window_orig = img_orig[r-window_ext:r+window_ext+1,
|
||||
c-window_ext:c+window_ext+1, :]
|
||||
|
||||
# weight pixels depending on distance to center pixel
|
||||
weights = gaussian_weights(window_ext, 3)
|
||||
weights = np.dstack((weights, weights, weights))
|
||||
|
||||
# compute sum of squared differences to all corners in warped image
|
||||
SSDs = []
|
||||
for cr, cc in coords_warped:
|
||||
window_warped = img_warped[cr-window_ext:cr+window_ext+1,
|
||||
cc-window_ext:cc+window_ext+1, :]
|
||||
SSD = np.sum(weights * (window_orig - window_warped)**2)
|
||||
SSDs.append(SSD)
|
||||
|
||||
# use corner with minimum SSD as correspondence
|
||||
min_idx = np.argmin(SSDs)
|
||||
return coords_warped_subpix[min_idx]
|
||||
|
||||
|
||||
# find correspondences using simple weighted sum of squared differences
|
||||
src = []
|
||||
dst = []
|
||||
for coord in coords_orig_subpix:
|
||||
src.append(coord)
|
||||
dst.append(match_corner(coord))
|
||||
src = np.array(src)
|
||||
dst = np.array(dst)
|
||||
|
||||
|
||||
# estimate affine transform model using all coordinates
|
||||
model = AffineTransform()
|
||||
model.estimate(src, dst)
|
||||
|
||||
# robustly estimate affine transform model with RANSAC
|
||||
model_robust, inliers = ransac((src, dst), AffineTransform, min_samples=3,
|
||||
residual_threshold=2, max_trials=100)
|
||||
outliers = inliers == False
|
||||
|
||||
|
||||
# compare "true" and estimated transform parameters
|
||||
print(tform.scale, tform.translation, tform.rotation)
|
||||
print(model.scale, model.translation, model.rotation)
|
||||
print(model_robust.scale, model_robust.translation, model_robust.rotation)
|
||||
|
||||
# visualize correspondence
|
||||
fig, ax = plt.subplots(nrows=2, ncols=1)
|
||||
|
||||
plt.gray()
|
||||
|
||||
inlier_idxs = np.nonzero(inliers)[0]
|
||||
plot_matches(ax[0], img_orig_gray, img_warped_gray, src, dst,
|
||||
np.column_stack((inlier_idxs, inlier_idxs)), matches_color='b')
|
||||
ax[0].axis('off')
|
||||
ax[0].set_title('Correct correspondences')
|
||||
|
||||
outlier_idxs = np.nonzero(outliers)[0]
|
||||
plot_matches(ax[1], img_orig_gray, img_warped_gray, src, dst,
|
||||
np.column_stack((outlier_idxs, outlier_idxs)), matches_color='r')
|
||||
ax[1].axis('off')
|
||||
ax[1].set_title('Faulty correspondences')
|
||||
|
||||
plt.show()
|
||||
@@ -0,0 +1,42 @@
|
||||
"""
|
||||
===============================
|
||||
Piecewise Affine Transformation
|
||||
===============================
|
||||
|
||||
This example shows how to use the Piecewise Affine Transformation.
|
||||
|
||||
"""
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from skimage.transform import PiecewiseAffineTransform, warp
|
||||
from skimage import data
|
||||
|
||||
|
||||
image = data.astronaut()
|
||||
rows, cols = image.shape[0], image.shape[1]
|
||||
|
||||
src_cols = np.linspace(0, cols, 20)
|
||||
src_rows = np.linspace(0, rows, 10)
|
||||
src_rows, src_cols = np.meshgrid(src_rows, src_cols)
|
||||
src = np.dstack([src_cols.flat, src_rows.flat])[0]
|
||||
|
||||
# add sinusoidal oscillation to row coordinates
|
||||
dst_rows = src[:, 1] - np.sin(np.linspace(0, 3 * np.pi, src.shape[0])) * 50
|
||||
dst_cols = src[:, 0]
|
||||
dst_rows *= 1.5
|
||||
dst_rows -= 1.5 * 50
|
||||
dst = np.vstack([dst_cols, dst_rows]).T
|
||||
|
||||
|
||||
tform = PiecewiseAffineTransform()
|
||||
tform.estimate(src, dst)
|
||||
|
||||
out_rows = image.shape[0] - 1.5 * 50
|
||||
out_cols = cols
|
||||
out = warp(image, tform, output_shape=(out_rows, out_cols))
|
||||
|
||||
fig, ax = plt.subplots()
|
||||
ax.imshow(out)
|
||||
ax.plot(tform.inverse(src)[:, 0], tform.inverse(src)[:, 1], '.b')
|
||||
ax.axis((0, out_cols, out_rows, 0))
|
||||
plt.show()
|
||||
@@ -0,0 +1,35 @@
|
||||
"""
|
||||
====================
|
||||
Build image pyramids
|
||||
====================
|
||||
|
||||
The `pyramid_gaussian` function takes an image and yields successive images
|
||||
shrunk by a constant scale factor. Image pyramids are often used, e.g., to
|
||||
implement algorithms for denoising, texture discrimination, and scale-
|
||||
invariant detection.
|
||||
|
||||
"""
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
from skimage import data
|
||||
from skimage.transform import pyramid_gaussian
|
||||
|
||||
|
||||
image = data.astronaut()
|
||||
rows, cols, dim = image.shape
|
||||
pyramid = tuple(pyramid_gaussian(image, downscale=2))
|
||||
|
||||
composite_image = np.zeros((rows, cols + cols / 2, 3), dtype=np.double)
|
||||
|
||||
composite_image[:rows, :cols, :] = pyramid[0]
|
||||
|
||||
i_row = 0
|
||||
for p in pyramid[1:]:
|
||||
n_rows, n_cols = p.shape[:2]
|
||||
composite_image[i_row:i_row + n_rows, cols:cols + n_cols] = p
|
||||
i_row += n_rows
|
||||
|
||||
fig, ax = plt.subplots()
|
||||
ax.imshow(composite_image)
|
||||
plt.show()
|
||||
@@ -0,0 +1,191 @@
|
||||
"""
|
||||
===============
|
||||
Radon transform
|
||||
===============
|
||||
|
||||
In computed tomography, the tomography reconstruction problem is to obtain
|
||||
a tomographic slice image from a set of projections [1]_. A projection is
|
||||
formed by drawing a set of parallel rays through the 2D object of interest,
|
||||
assigning the integral of the object's contrast along each ray to a single
|
||||
pixel in the projection. A single projection of a 2D object is one dimensional.
|
||||
To enable computed tomography reconstruction of the object, several projections
|
||||
must be acquired, each of them corresponding to a different angle between the
|
||||
rays with respect to the object. A collection of projections at several angles
|
||||
is called a sinogram, which is a linear transform of the original image.
|
||||
|
||||
The inverse Radon transform is used in computed tomography to reconstruct
|
||||
a 2D image from the measured projections (the sinogram). A practical, exact
|
||||
implementation of the inverse Radon transform does not exist, but there are
|
||||
several good approximate algorithms available.
|
||||
|
||||
As the inverse Radon transform reconstructs the object from a set of
|
||||
projections, the (forward) Radon transform can be used to simulate a
|
||||
tomography experiment.
|
||||
|
||||
This script performs the Radon transform to simulate a tomography experiment
|
||||
and reconstructs the input image based on the resulting sinogram formed by
|
||||
the simulation. Two methods for performing the inverse Radon transform
|
||||
and reconstructing the original image are compared: The Filtered Back
|
||||
Projection (FBP) and the Simultaneous Algebraic Reconstruction
|
||||
Technique (SART).
|
||||
|
||||
For further information on tomographic reconstruction, see
|
||||
|
||||
- AC Kak, M Slaney, "Principles of Computerized Tomographic Imaging",
|
||||
http://www.slaney.org/pct/pct-toc.html
|
||||
- http://en.wikipedia.org/wiki/Radon_transform
|
||||
|
||||
The forward transform
|
||||
=====================
|
||||
|
||||
As our original image, we will use the Shepp-Logan phantom. When calculating
|
||||
the Radon transform, we need to decide how many projection angles we wish
|
||||
to use. As a rule of thumb, the number of projections should be about the
|
||||
same as the number of pixels there are across the object (to see why this
|
||||
is so, consider how many unknown pixel values must be determined in the
|
||||
reconstruction process and compare this to the number of measurements
|
||||
provided by the projections), and we follow that rule here. Below is the
|
||||
original image and its Radon transform, often known as its _sinogram_:
|
||||
"""
|
||||
|
||||
from __future__ import print_function, division
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
from skimage.io import imread
|
||||
from skimage import data_dir
|
||||
from skimage.transform import radon, rescale
|
||||
|
||||
image = imread(data_dir + "/phantom.png", as_grey=True)
|
||||
image = rescale(image, scale=0.4)
|
||||
|
||||
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5))
|
||||
|
||||
ax1.set_title("Original")
|
||||
ax1.imshow(image, cmap=plt.cm.Greys_r)
|
||||
|
||||
theta = np.linspace(0., 180., max(image.shape), endpoint=False)
|
||||
sinogram = radon(image, theta=theta, circle=True)
|
||||
ax2.set_title("Radon transform\n(Sinogram)")
|
||||
ax2.set_xlabel("Projection angle (deg)")
|
||||
ax2.set_ylabel("Projection position (pixels)")
|
||||
ax2.imshow(sinogram, cmap=plt.cm.Greys_r,
|
||||
extent=(0, 180, 0, sinogram.shape[0]), aspect='auto')
|
||||
|
||||
fig.subplots_adjust(hspace=0.4, wspace=0.5)
|
||||
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.
|
||||
|
||||
"""
|
||||
|
||||
from skimage.transform import iradon
|
||||
|
||||
reconstruction_fbp = iradon(sinogram, theta=theta, circle=True)
|
||||
error = reconstruction_fbp - image
|
||||
print('FBP rms reconstruction error: %.3g' % np.sqrt(np.mean(error**2)))
|
||||
|
||||
imkwargs = dict(vmin=-0.2, vmax=0.2)
|
||||
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'})
|
||||
ax1.set_title("Reconstruction\nFiltered back projection")
|
||||
ax1.imshow(reconstruction_fbp, cmap=plt.cm.Greys_r)
|
||||
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.
|
||||
|
||||
"""
|
||||
|
||||
from skimage.transform import iradon_sart
|
||||
|
||||
reconstruction_sart = iradon_sart(sinogram, theta=theta)
|
||||
error = reconstruction_sart - image
|
||||
print('SART (1 iteration) rms reconstruction error: %.3g'
|
||||
% np.sqrt(np.mean(error**2)))
|
||||
|
||||
fig, ax = plt.subplots(2, 2, figsize=(8, 8.5), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'})
|
||||
ax1, ax2, ax3, ax4 = ax.ravel()
|
||||
ax1.set_title("Reconstruction\nSART")
|
||||
ax1.imshow(reconstruction_sart, cmap=plt.cm.Greys_r)
|
||||
ax2.set_title("Reconstruction error\nSART")
|
||||
ax2.imshow(reconstruction_sart - image, cmap=plt.cm.Greys_r, **imkwargs)
|
||||
|
||||
# Run a second iteration of SART by supplying the reconstruction
|
||||
# from the first iteration as an initial estimate
|
||||
reconstruction_sart2 = iradon_sart(sinogram, theta=theta,
|
||||
image=reconstruction_sart)
|
||||
error = reconstruction_sart2 - image
|
||||
print('SART (2 iterations) rms reconstruction error: %.3g'
|
||||
% np.sqrt(np.mean(error**2)))
|
||||
|
||||
ax3.set_title("Reconstruction\nSART, 2 iterations")
|
||||
ax3.imshow(reconstruction_sart2, cmap=plt.cm.Greys_r)
|
||||
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)
|
||||
|
||||
"""
|
||||
@@ -0,0 +1,56 @@
|
||||
"""
|
||||
=========================================
|
||||
Robust line model estimation using RANSAC
|
||||
=========================================
|
||||
|
||||
In this example we see how to robustly fit a line model to faulty data using
|
||||
the RANSAC algorithm.
|
||||
|
||||
"""
|
||||
import numpy as np
|
||||
from matplotlib import pyplot as plt
|
||||
|
||||
from skimage.measure import LineModel, ransac
|
||||
|
||||
|
||||
np.random.seed(seed=1)
|
||||
|
||||
# generate coordinates of line
|
||||
x = np.arange(-200, 200)
|
||||
y = 0.2 * x + 20
|
||||
data = np.column_stack([x, y])
|
||||
|
||||
# add faulty data
|
||||
faulty = np.array(30 * [(180., -100)])
|
||||
faulty += 5 * np.random.normal(size=faulty.shape)
|
||||
data[:faulty.shape[0]] = faulty
|
||||
|
||||
# add gaussian noise to coordinates
|
||||
noise = np.random.normal(size=data.shape)
|
||||
data += 0.5 * noise
|
||||
data[::2] += 5 * noise[::2]
|
||||
data[::4] += 20 * noise[::4]
|
||||
|
||||
# fit line using all data
|
||||
model = LineModel()
|
||||
model.estimate(data)
|
||||
|
||||
# robustly fit line only using inlier data with RANSAC algorithm
|
||||
model_robust, inliers = ransac(data, LineModel, min_samples=2,
|
||||
residual_threshold=1, max_trials=1000)
|
||||
outliers = inliers == False
|
||||
|
||||
# generate coordinates of estimated models
|
||||
line_x = np.arange(-250, 250)
|
||||
line_y = model.predict_y(line_x)
|
||||
line_y_robust = model_robust.predict_y(line_x)
|
||||
|
||||
fig, ax = plt.subplots()
|
||||
ax.plot(data[inliers, 0], data[inliers, 1], '.b', alpha=0.6,
|
||||
label='Inlier data')
|
||||
ax.plot(data[outliers, 0], data[outliers, 1], '.r', alpha=0.6,
|
||||
label='Outlier data')
|
||||
ax.plot(line_x, line_y, '-k', label='Line model from all data')
|
||||
ax.plot(line_x, line_y_robust, '-b', label='Robust line model')
|
||||
ax.legend(loc='lower left')
|
||||
plt.show()
|
||||
@@ -0,0 +1,40 @@
|
||||
"""
|
||||
============================================
|
||||
Robust 3D line model estimation using RANSAC
|
||||
============================================
|
||||
|
||||
In this example we see how to robustly fit a 3D line model to faulty data using
|
||||
the RANSAC algorithm.
|
||||
|
||||
"""
|
||||
import numpy as np
|
||||
from matplotlib import pyplot as plt
|
||||
from mpl_toolkits.mplot3d import Axes3D
|
||||
from skimage.measure import LineModelND, ransac
|
||||
|
||||
np.random.seed(seed=1)
|
||||
|
||||
# generate coordinates of line
|
||||
point = np.array([0, 0, 0], dtype='float')
|
||||
direction = np.array([1, 1, 1], dtype='float') / np.sqrt(3)
|
||||
xyz = point + 10 * np.arange(-100, 100)[..., np.newaxis] * direction
|
||||
|
||||
# add gaussian noise to coordinates
|
||||
noise = np.random.normal(size=xyz.shape)
|
||||
xyz += 0.5 * noise
|
||||
xyz[::2] += 20 * noise[::2]
|
||||
xyz[::4] += 100 * noise[::4]
|
||||
|
||||
# robustly fit line only using inlier data with RANSAC algorithm
|
||||
model_robust, inliers = ransac(xyz, LineModelND, min_samples=2,
|
||||
residual_threshold=1, max_trials=1000)
|
||||
outliers = inliers == False
|
||||
|
||||
fig = plt.figure()
|
||||
ax = fig.add_subplot(111, projection='3d')
|
||||
ax.scatter(xyz[inliers][:, 0], xyz[inliers][:, 1], xyz[inliers][:, 2], c='b',
|
||||
marker='o', label='Inlier data')
|
||||
ax.scatter(xyz[outliers][:, 0], xyz[outliers][:, 1], xyz[outliers][:, 2], c='r',
|
||||
marker='o', label='Outlier data')
|
||||
ax.legend(loc='lower left')
|
||||
plt.show()
|
||||
@@ -0,0 +1,91 @@
|
||||
"""
|
||||
=====================================
|
||||
Cross-Correlation (Phase Correlation)
|
||||
=====================================
|
||||
|
||||
In this example, we use phase correlation to identify the relative shift
|
||||
between two similar-sized images.
|
||||
|
||||
The ``register_translation`` function uses cross-correlation in Fourier space,
|
||||
optionally employing an upsampled matrix-multiplication DFT to achieve
|
||||
arbitrary subpixel precision. [1]_
|
||||
|
||||
.. [1] Manuel Guizar-Sicairos, Samuel T. Thurman, and James R. Fienup,
|
||||
"Efficient subpixel image registration algorithms," Optics Letters 33,
|
||||
156-158 (2008).
|
||||
|
||||
"""
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
from skimage import data
|
||||
from skimage.feature import register_translation
|
||||
from skimage.feature.register_translation import _upsampled_dft
|
||||
from scipy.ndimage import fourier_shift
|
||||
|
||||
image = data.camera()
|
||||
shift = (-2.4, 1.32)
|
||||
# (-2.4, 1.32) pixel offset relative to reference coin
|
||||
offset_image = fourier_shift(np.fft.fftn(image), shift)
|
||||
offset_image = np.fft.ifftn(offset_image)
|
||||
print("Known offset (y, x):")
|
||||
print(shift)
|
||||
|
||||
# pixel precision first
|
||||
shift, error, diffphase = register_translation(image, offset_image)
|
||||
|
||||
fig = plt.figure(figsize=(8, 3))
|
||||
ax1 = plt.subplot(1, 3, 1, adjustable='box-forced')
|
||||
ax2 = plt.subplot(1, 3, 2, sharex=ax1, sharey=ax1, adjustable='box-forced')
|
||||
ax3 = plt.subplot(1, 3, 3)
|
||||
|
||||
ax1.imshow(image)
|
||||
ax1.set_axis_off()
|
||||
ax1.set_title('Reference image')
|
||||
|
||||
ax2.imshow(offset_image.real)
|
||||
ax2.set_axis_off()
|
||||
ax2.set_title('Offset image')
|
||||
|
||||
# View the output of a cross-correlation to show what the algorithm is
|
||||
# doing behind the scenes
|
||||
image_product = np.fft.fft2(image) * np.fft.fft2(offset_image).conj()
|
||||
cc_image = np.fft.fftshift(np.fft.ifft2(image_product))
|
||||
ax3.imshow(cc_image.real)
|
||||
ax3.set_axis_off()
|
||||
ax3.set_title("Cross-correlation")
|
||||
|
||||
plt.show()
|
||||
|
||||
print("Detected pixel offset (y, x):")
|
||||
print(shift)
|
||||
|
||||
# subpixel precision
|
||||
shift, error, diffphase = register_translation(image, offset_image, 100)
|
||||
|
||||
fig = plt.figure(figsize=(8, 3))
|
||||
ax1 = plt.subplot(1, 3, 1, adjustable='box-forced')
|
||||
ax2 = plt.subplot(1, 3, 2, sharex=ax1, sharey=ax1, adjustable='box-forced')
|
||||
ax3 = plt.subplot(1, 3, 3)
|
||||
|
||||
ax1.imshow(image)
|
||||
ax1.set_axis_off()
|
||||
ax1.set_title('Reference image')
|
||||
|
||||
ax2.imshow(offset_image.real)
|
||||
ax2.set_axis_off()
|
||||
ax2.set_title('Offset image')
|
||||
|
||||
# Calculate the upsampled DFT, again to show what the algorithm is doing
|
||||
# behind the scenes. Constants correspond to calculated values in routine.
|
||||
# See source code for details.
|
||||
cc_image = _upsampled_dft(image_product, 150, 100, (shift*100)+75).conj()
|
||||
ax3.imshow(cc_image.real)
|
||||
ax3.set_axis_off()
|
||||
ax3.set_title("Supersampled XC sub-area")
|
||||
|
||||
|
||||
plt.show()
|
||||
|
||||
print("Detected subpixel offset (y, x):")
|
||||
print(shift)
|
||||
@@ -0,0 +1,96 @@
|
||||
"""
|
||||
============
|
||||
Seam Carving
|
||||
============
|
||||
|
||||
This example demonstrates how images can be resized using seam carving [1]_.
|
||||
Resizing to a new aspect ratio distorts image contents. Seam carving attempts
|
||||
to resize *without* distortion, by removing regions of an image which are less
|
||||
important. In this example we are using the Sobel filter to signify the
|
||||
importance of each pixel.
|
||||
|
||||
.. [1] Shai Avidan and Ariel Shamir
|
||||
"Seam Carving for Content-Aware Image Resizing"
|
||||
http://www.cs.jhu.edu/~misha/ReadingSeminar/Papers/Avidan07.pdf
|
||||
|
||||
"""
|
||||
from skimage import data, draw
|
||||
from skimage import transform, util
|
||||
import numpy as np
|
||||
from skimage import filters, color
|
||||
from matplotlib import pyplot as plt
|
||||
|
||||
|
||||
hl_color = np.array([0, 1, 0])
|
||||
|
||||
img = data.rocket()
|
||||
img = util.img_as_float(img)
|
||||
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.
|
||||
|
||||
"""
|
||||
|
||||
masked_img = img.copy()
|
||||
|
||||
poly = [(404, 281), (404, 360), (359, 364), (338, 337), (145, 337), (120, 322),
|
||||
(145, 304), (340, 306), (362, 284)]
|
||||
pr = np.array([p[0] for p in poly])
|
||||
pc = np.array([p[1] for p in poly])
|
||||
rr, cc = draw.polygon(pr, pc)
|
||||
|
||||
masked_img[rr, cc, :] = masked_img[rr, cc, :]*0.5 + hl_color*.5
|
||||
plt.figure()
|
||||
plt.title('Object Marked')
|
||||
|
||||
plt.imshow(masked_img)
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
"""
|
||||
|
||||
eimg[rr, cc] -= 1000
|
||||
|
||||
plt.figure()
|
||||
plt.title('Object Removed')
|
||||
out = transform.seam_carve(img, eimg, 'vertical', 90)
|
||||
resized = transform.resize(img, out.shape)
|
||||
plt.imshow(out)
|
||||
plt.show()
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
"""
|
||||
@@ -0,0 +1,75 @@
|
||||
"""
|
||||
===========================
|
||||
Structural similarity index
|
||||
===========================
|
||||
|
||||
When comparing images, the mean squared error (MSE)--while simple to
|
||||
implement--is not highly indicative of perceived similarity. Structural
|
||||
similarity aims to address this shortcoming by taking texture into account
|
||||
[1]_, [2]_.
|
||||
|
||||
The example shows two modifications of the input image, each with the same MSE,
|
||||
but with very different mean structural similarity indices.
|
||||
|
||||
.. [1] Zhou Wang; Bovik, A.C.; ,"Mean squared error: Love it or leave it? A new
|
||||
look at Signal Fidelity Measures," Signal Processing Magazine, IEEE,
|
||||
vol. 26, no. 1, pp. 98-117, Jan. 2009.
|
||||
|
||||
.. [2] Z. Wang, A. C. Bovik, H. R. Sheikh and E. P. Simoncelli, "Image quality
|
||||
assessment: From error visibility to structural similarity," IEEE
|
||||
Transactions on Image Processing, vol. 13, no. 4, pp. 600-612,
|
||||
Apr. 2004.
|
||||
|
||||
"""
|
||||
import numpy as np
|
||||
import matplotlib
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
from skimage import data, img_as_float
|
||||
from skimage.measure import structural_similarity as ssim
|
||||
|
||||
|
||||
matplotlib.rcParams['font.size'] = 9
|
||||
|
||||
|
||||
img = img_as_float(data.camera())
|
||||
rows, cols = img.shape
|
||||
|
||||
noise = np.ones_like(img) * 0.2 * (img.max() - img.min())
|
||||
noise[np.random.random(size=noise.shape) > 0.5] *= -1
|
||||
|
||||
|
||||
def mse(x, y):
|
||||
return np.linalg.norm(x - y)
|
||||
|
||||
img_noise = img + noise
|
||||
img_const = img + abs(noise)
|
||||
|
||||
fig, (ax0, ax1, ax2) = plt.subplots(nrows=1, ncols=3, figsize=(8, 4), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'})
|
||||
|
||||
mse_none = mse(img, img)
|
||||
ssim_none = ssim(img, img, dynamic_range=img.max() - img.min())
|
||||
|
||||
mse_noise = mse(img, img_noise)
|
||||
ssim_noise = ssim(img, img_noise,
|
||||
dynamic_range=img_const.max() - img_const.min())
|
||||
|
||||
mse_const = mse(img, img_const)
|
||||
ssim_const = ssim(img, img_const,
|
||||
dynamic_range=img_noise.max() - img_noise.min())
|
||||
|
||||
label = 'MSE: %2.f, SSIM: %.2f'
|
||||
|
||||
ax0.imshow(img, cmap=plt.cm.gray, vmin=0, vmax=1)
|
||||
ax0.set_xlabel(label % (mse_none, ssim_none))
|
||||
ax0.set_title('Original image')
|
||||
|
||||
ax1.imshow(img_noise, cmap=plt.cm.gray, vmin=0, vmax=1)
|
||||
ax1.set_xlabel(label % (mse_noise, ssim_noise))
|
||||
ax1.set_title('Image with noise')
|
||||
|
||||
ax2.imshow(img_const, cmap=plt.cm.gray, vmin=0, vmax=1)
|
||||
ax2.set_xlabel(label % (mse_const, ssim_const))
|
||||
ax2.set_title('Image plus constant')
|
||||
|
||||
plt.show()
|
||||
@@ -0,0 +1,84 @@
|
||||
"""
|
||||
=====
|
||||
Swirl
|
||||
=====
|
||||
|
||||
Image swirling is a non-linear image deformation that creates a whirlpool
|
||||
effect. This example describes the implementation of this transform in
|
||||
``skimage``, as well as the underlying warp mechanism.
|
||||
|
||||
Image warping
|
||||
-------------
|
||||
When applying a geometric transformation on an image, we typically make use of
|
||||
a reverse mapping, i.e., for each pixel in the output image, we compute its
|
||||
corresponding position in the input. The reason is that, if we were to do it
|
||||
the other way around (map each input pixel to its new output position), some
|
||||
pixels in the output may be left empty. On the other hand, each output
|
||||
coordinate has exactly one corresponding location in (or outside) the input
|
||||
image, and even if that position is non-integer, we may use interpolation to
|
||||
compute the corresponding image value.
|
||||
|
||||
Performing a reverse mapping
|
||||
----------------------------
|
||||
To perform a geometric warp in ``skimage``, you simply need to provide the
|
||||
reverse mapping to the ``skimage.transform.warp`` function. E.g., consider the
|
||||
case where we would like to shift an image 50 pixels to the left. The reverse
|
||||
mapping for such a shift would be::
|
||||
|
||||
def shift_left(xy):
|
||||
xy[:, 0] += 50
|
||||
return xy
|
||||
|
||||
The corresponding call to warp is::
|
||||
|
||||
from skimage.transform import warp
|
||||
warp(image, shift_left)
|
||||
|
||||
The swirl transformation
|
||||
------------------------
|
||||
Consider the coordinate :math:`(x, y)` in the output image. The reverse
|
||||
mapping for the swirl transformation first computes, relative to a center
|
||||
:math:`(x_0, y_0)`, its polar coordinates,
|
||||
|
||||
.. math::
|
||||
|
||||
\\theta = \\arctan(y/x)
|
||||
|
||||
\\rho = \sqrt{(x - x_0)^2 + (y - y_0)^2},
|
||||
|
||||
and then transforms them according to
|
||||
|
||||
.. math::
|
||||
|
||||
r = \ln(2) \, \mathtt{radius} / 5
|
||||
|
||||
\phi = \mathtt{rotation}
|
||||
|
||||
s = \mathtt{strength}
|
||||
|
||||
\\theta' = \phi + s \, e^{-\\rho / r + \\theta}
|
||||
|
||||
where ``strength`` is a parameter for the amount of swirl, ``radius`` indicates
|
||||
the swirl extent in pixels, and ``rotation`` adds a rotation angle. The
|
||||
transformation of ``radius`` into :math:`r` is to ensure that the
|
||||
transformation decays to :math:`\\approx 1/1000^{\mathsf{th}}` within the
|
||||
specified radius.
|
||||
|
||||
"""
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
from skimage import data
|
||||
from skimage.transform import swirl
|
||||
|
||||
|
||||
image = data.checkerboard()
|
||||
swirled = swirl(image, rotation=0, strength=10, radius=120, order=2)
|
||||
|
||||
fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(8, 3), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'})
|
||||
|
||||
ax0.imshow(image, cmap=plt.cm.gray, interpolation='none')
|
||||
ax0.axis('off')
|
||||
ax1.imshow(swirled, cmap=plt.cm.gray, interpolation='none')
|
||||
ax1.axis('off')
|
||||
|
||||
plt.show()
|
||||
Reference in New Issue
Block a user