wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-08-07 11:28:14 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1691 from juliusbierk/active-contour-model

Browse Source 
New feature: Active contour model
This commit is contained in:
Steven Silvester
2015-12-13 14:57:10 -06:00
parent 827e4b11c9 fb9e9e6f7a
commit 5b8fd77b80
4 changed files with 451 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/examples/plot_active_contours.py
+99
View File
@@ -0,0 +1,99 @@
"""
====================
Active Contour Model
====================
The active contour model is a method to fit open or closed splines to lines or
edges in an image. It works by minimising an energy that is in part defined by
the image and part by the spline's shape: length and smoothness. The
minimization is done implicitly in the shape energy and explicitly in the
image energy.
In the following two examples the active contour model is used (1) to segment
the face of a person from the rest of an image by fitting a closed curve
to the edges of the face and (2) to find the darkest curve between two fixed
points while obeying smoothness considerations. Typically it is a good idea to
smooth images a bit before analyzing, as done in the following examples.
.. [1] *Snakes: Active contour models*. Kass, M.; Witkin, A.; Terzopoulos, D.
International Journal of Computer Vision 1 (4): 321 (1988).
We initialize a circle around the astronaut's face and use the default boundary
condition ``bc='periodic'`` to fit a closed curve. The default parameters
``w_line=0, w_edge=1`` will make the curve search towards edges, such as the
boundaries of the face.
"""
import numpy as np
import matplotlib.pyplot as plt
from skimage.color import rgb2gray
from skimage import data
from skimage.filters import gaussian_filter
from skimage.segmentation import active_contour
# Test scipy version, since active contour is only possible
# with recent scipy version
import scipy
scipy_version = list(map(int, scipy.__version__.split('.')))
new_scipy = scipy_version[0] > 0 or \
(scipy_version[0] == 0 and scipy_version[1] >= 14)
img = data.astronaut()
img = rgb2gray(img)
s = np.linspace(0, 2*np.pi, 400)
x = 220 + 100*np.cos(s)
y = 100 + 100*np.sin(s)
init = np.array([x, y]).T
if not new_scipy:
print('You are using an old version of scipy. '
'Active contours is implemented for scipy versions '
'0.14.0 and above.')
if new_scipy:
snake = active_contour(gaussian_filter(img, 3),
init, alpha=0.015, beta=10, gamma=0.001)
fig = plt.figure(figsize=(7, 7))
ax = fig.add_subplot(111)
plt.gray()
ax.imshow(img)
ax.plot(init[:, 0], init[:, 1], '--r')
ax.plot(snake[:, 0], snake[:, 1], '-b')
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.
"""
img = data.text()
x = np.linspace(5, 424, 100)
y = np.linspace(136, 50, 100)
init = np.array([x, y]).T
if new_scipy:
snake = active_contour(gaussian_filter(img, 1), init, bc='fixed',
alpha=0.1, beta=1.0, w_line=-5, w_edge=0, gamma=0.1)
fig = plt.figure(figsize=(9, 5))
ax = fig.add_subplot(111)
plt.gray()
ax.imshow(img)
ax.plot(init[:, 0], init[:, 1], '--r')
ax.plot(snake[:, 0], snake[:, 1], '-b')
ax.set_xticks([]), ax.set_yticks([])
ax.axis([0, img.shape[1], img.shape[0], 0])
plt.show()
"""
.. image:: PLOT2RST.current_figure
"""
skimage/segmentation/__init__.py
+2
View File
@@ -1,4 +1,5 @@
from .random_walker_segmentation import random_walker
from .active_contour_model import active_contour
from ._felzenszwalb import felzenszwalb
from .slic_superpixels import slic
from ._quickshift import quickshift
@@ -8,6 +9,7 @@ from ._join import join_segmentations, relabel_from_one, relabel_sequential
__all__ = ['random_walker',
'active_contour',
'felzenszwalb',
'slic',
'quickshift',
skimage/segmentation/active_contour_model.py
+234
View File
@@ -0,0 +1,234 @@
import numpy as np
from skimage import img_as_float
import scipy
import scipy.linalg
from scipy.interpolate import RectBivariateSpline, interp2d
from skimage.filters import sobel
def active_contour(image, snake, alpha=0.01, beta=0.1,
w_line=0, w_edge=1, gamma=0.01,
bc='periodic', max_px_move=1.0,
max_iterations=2500, convergence=0.1):
"""Active contour model.
Active contours by fitting snakes to features of images. Supports single
and multichannel 2D images. Snakes can be periodic (for segmentation) or
have fixed and/or free ends.
Parameters
----------
image : (N, M) or (N, M, 3) ndarray
Input image.
snake : (N, 2) ndarray
Initialisation coordinates of snake. For periodic snakes, it should
not include duplicate endpoints.
alpha : float, optional
Snake length shape parameter. Higher values makes snake contract
faster.
beta : float, optional
Snake smoothness shape parameter. Higher values makes snake smoother.
w_line : float, optional
Controls attraction to brightness. Use negative values to attract to
dark regions.
w_edge : float, optional
Controls attraction to edges. Use negative values to repel snake from
edges.
gamma : float, optional
Explicit time stepping parameter.
bc : {'periodic', 'free', 'fixed'}, optional
Boundary conditions for worm. 'periodic' attaches the two ends of the
snake, 'fixed' holds the end-points in place, and'free' allows free
movement of the ends. 'fixed' and 'free' can be combined by parsing
'fixed-free', 'free-fixed'. Parsing 'fixed-fixed' or 'free-free'
yields same behaviour as 'fixed' and 'free', respectively.
max_px_move : float, optional
Maximum pixel distance to move per iteration.
max_iterations : int, optional
Maximum iterations to optimize snake shape.
convergence: float, optional
Convergence criteria.
Returns
-------
snake : (N, 2) ndarray
Optimised snake, same shape as input parameter.
References
----------
.. [1] Kass, M.; Witkin, A.; Terzopoulos, D. "Snakes: Active contour
models". International Journal of Computer Vision 1 (4): 321
(1988).
Examples
--------
>>> from skimage.draw import circle_perimeter
>>> from skimage.filters import gaussian_filter
Create and smooth image:
>>> img = np.zeros((100, 100))
>>> rr, cc = circle_perimeter(35, 45, 25)
>>> img[rr, cc] = 1
>>> img = gaussian_filter(img, 2)
Initiliaze spline:
>>> s = np.linspace(0, 2*np.pi,100)
>>> init = 50*np.array([np.cos(s), np.sin(s)]).T+50
Fit spline to image:
>>> snake = active_contour(img, init, w_edge=0, w_line=1) #doctest: +SKIP
>>> dist = np.sqrt((45-snake[:, 0])**2 +(35-snake[:, 1])**2) #doctest: +SKIP
>>> int(np.mean(dist)) #doctest: +SKIP
25
"""
scipy_version = list(map(int, scipy.__version__.split('.')))
new_scipy = scipy_version[0] > 0 or \
(scipy_version[0] == 0 and scipy_version[1] >= 14)
if not new_scipy:
raise NotImplementedError('You are using an old version of scipy. '
'Active contours is implemented for scipy versions '
'0.14.0 and above.')
max_iterations = int(max_iterations)
if max_iterations <= 0:
raise ValueError("max_iterations should be >0.")
convergence_order = 10
valid_bcs = ['periodic', 'free', 'fixed', 'free-fixed',
'fixed-free', 'fixed-fixed', 'free-free']
if bc not in valid_bcs:
raise ValueError("Invalid boundary condition.\n" +
"Should be one of: "+", ".join(valid_bcs)+'.')
img = img_as_float(image)
RGB = img.ndim == 3
# Find edges using sobel:
if w_edge != 0:
if RGB:
edge = [sobel(img[:, :, 0]), sobel(img[:, :, 1]),
sobel(img[:, :, 2])]
else:
edge = [sobel(img)]
for i in range(3 if RGB else 1):
edge[i][0, :] = edge[i][1, :]
edge[i][-1, :] = edge[i][-2, :]
edge[i][:, 0] = edge[i][:, 1]
edge[i][:, -1] = edge[i][:, -2]
else:
edge = [0]
# Superimpose intensity and edge images:
if RGB:
img = w_line*np.sum(img, axis=2) \
+ w_edge*sum(edge)
else:
img = w_line*img + w_edge*edge[0]
# Interpolate for smoothness:
if new_scipy:
intp = RectBivariateSpline(np.arange(img.shape[1]),
np.arange(img.shape[0]),
img.T, kx=2, ky=2, s=0)
else:
intp = np.vectorize(interp2d(np.arange(img.shape[1]),
np.arange(img.shape[0]), img, kind='cubic',
copy=False, bounds_error=False, fill_value=0))
x, y = snake[:, 0].copy(), snake[:, 1].copy()
xsave = np.empty((convergence_order, len(x)))
ysave = np.empty((convergence_order, len(x)))
# Build snake shape matrix for Euler equation
n = len(x)
a = np.roll(np.eye(n), -1, axis=0) + \
np.roll(np.eye(n), -1, axis=1) - \
2*np.eye(n) # second order derivative, central difference
b = np.roll(np.eye(n), -2, axis=0) + \
np.roll(np.eye(n), -2, axis=1) - \
4*np.roll(np.eye(n), -1, axis=0) - \
4*np.roll(np.eye(n), -1, axis=1) + \
6*np.eye(n) # fourth order derivative, central difference
A = -alpha*a + beta*b
# Impose boundary conditions different from periodic:
sfixed = False
if bc.startswith('fixed'):
A[0, :] = 0
A[1, :] = 0
A[1, :3] = [1, -2, 1]
sfixed = True
efixed = False
if bc.endswith('fixed'):
A[-1, :] = 0
A[-2, :] = 0
A[-2, -3:] = [1, -2, 1]
efixed = True
sfree = False
if bc.startswith('free'):
A[0, :] = 0
A[0, :3] = [1, -2, 1]
A[1, :] = 0
A[1, :4] = [-1, 3, -3, 1]
sfree = True
efree = False
if bc.endswith('free'):
A[-1, :] = 0
A[-1, -3:] = [1, -2, 1]
A[-2, :] = 0
A[-2, -4:] = [-1, 3, -3, 1]
efree = True
# Only one inversion is needed for implicit spline energy minimization:
inv = scipy.linalg.inv(A+gamma*np.eye(n))
# Explicit time stepping for image energy minimization:
for i in range(max_iterations):
if new_scipy:
fx = intp(x, y, dx=1, grid=False)
fy = intp(x, y, dy=1, grid=False)
else:
fx = intp(x, y, dx=1)
fy = intp(x, y, dy=1)
if sfixed:
fx[0] = 0
fy[0] = 0
if efixed:
fx[-1] = 0
fy[-1] = 0
if sfree:
fx[0] *= 2
fy[0] *= 2
if efree:
fx[-1] *= 2
fy[-1] *= 2
xn = np.dot(inv, gamma*x + fx)
yn = np.dot(inv, gamma*y + fy)
# Movements are capped to max_px_move per iteration:
dx = max_px_move*np.tanh(xn-x)
dy = max_px_move*np.tanh(yn-y)
if sfixed:
dx[0] = 0
dy[0] = 0
if efixed:
dx[-1] = 0
dy[-1] = 0
x[:] += dx
y[:] += dy
# Convergence criteria needs to compare to a number of previous
# configurations since oscillations can occur.
j = i % (convergence_order+1)
if j < convergence_order:
xsave[j, :] = x
ysave[j, :] = y
else:
dist = np.min(np.max(np.abs(xsave-x[None, :]) +
np.abs(ysave-y[None, :]), 1))
if dist < convergence:
break
return np.array([x, y]).T
skimage/segmentation/tests/test_active_contour_model.py
+116
View File
@@ -0,0 +1,116 @@
import numpy as np
from skimage import data
from skimage.color import rgb2gray
from skimage.filters import gaussian_filter
from skimage.segmentation import active_contour
from numpy.testing import assert_equal, assert_allclose, assert_raises
from numpy.testing.decorators import skipif
import scipy
scipy_version = list(map(int, scipy.__version__.split('.')))
new_scipy = scipy_version[0] > 0 or \
(scipy_version[0] == 0 and scipy_version[1] >= 14)
@skipif(not new_scipy)
def test_periodic_reference():
img = data.astronaut()
img = rgb2gray(img)
s = np.linspace(0, 2*np.pi, 400)
x = 220 + 100*np.cos(s)
y = 100 + 100*np.sin(s)
init = np.array([x, y]).T
snake = active_contour(gaussian_filter(img, 3), init,
alpha=0.015, beta=10, w_line=0, w_edge=1, gamma=0.001)
refx = [299, 298, 298, 298, 298, 297, 297, 296, 296, 295]
refy = [98, 99, 100, 101, 102, 103, 104, 105, 106, 108]
assert_equal(np.array(snake[:10, 0], dtype=np.int32), refx)
assert_equal(np.array(snake[:10, 1], dtype=np.int32), refy)
@skipif(not new_scipy)
def test_fixed_reference():
img = data.text()
x = np.linspace(5, 424, 100)
y = np.linspace(136, 50, 100)
init = np.array([x, y]).T
snake = active_contour(gaussian_filter(img, 1), init, bc='fixed',
alpha=0.1, beta=1.0, w_line=-5, w_edge=0, gamma=0.1)
refx = [5, 9, 13, 17, 21, 25, 30, 34, 38, 42]
refy = [136, 135, 134, 133, 132, 131, 129, 128, 127, 125]
assert_equal(np.array(snake[:10, 0], dtype=np.int32), refx)
assert_equal(np.array(snake[:10, 1], dtype=np.int32), refy)
@skipif(not new_scipy)
def test_free_reference():
img = data.text()
x = np.linspace(5, 424, 100)
y = np.linspace(70, 40, 100)
init = np.array([x, y]).T
snake = active_contour(gaussian_filter(img, 3), init, bc='free',
alpha=0.1, beta=1.0, w_line=-5, w_edge=0, gamma=0.1)
refx = [10, 13, 16, 19, 23, 26, 29, 32, 36, 39]
refy = [76, 76, 75, 74, 73, 72, 71, 70, 69, 69]
assert_equal(np.array(snake[:10, 0], dtype=np.int32), refx)
assert_equal(np.array(snake[:10, 1], dtype=np.int32), refy)
@skipif(not new_scipy)
def test_RGB():
img = gaussian_filter(data.text(), 1)
imgR = np.zeros((img.shape[0], img.shape[1], 3))
imgG = np.zeros((img.shape[0], img.shape[1], 3))
imgRGB = np.zeros((img.shape[0], img.shape[1], 3))
imgR[:, :, 0] = img
imgG[:, :, 1] = img
imgRGB[:, :, :] = img[:, :, None]
x = np.linspace(5, 424, 100)
y = np.linspace(136, 50, 100)
init = np.array([x, y]).T
snake = active_contour(imgR, init, bc='fixed',
alpha=0.1, beta=1.0, w_line=-5, w_edge=0, gamma=0.1)
refx = [5, 9, 13, 17, 21, 25, 30, 34, 38, 42]
refy = [136, 135, 134, 133, 132, 131, 129, 128, 127, 125]
assert_equal(np.array(snake[:10, 0], dtype=np.int32), refx)
assert_equal(np.array(snake[:10, 1], dtype=np.int32), refy)
snake = active_contour(imgG, init, bc='fixed',
alpha=0.1, beta=1.0, w_line=-5, w_edge=0, gamma=0.1)
assert_equal(np.array(snake[:10, 0], dtype=np.int32), refx)
assert_equal(np.array(snake[:10, 1], dtype=np.int32), refy)
snake = active_contour(imgRGB, init, bc='fixed',
alpha=0.1, beta=1.0, w_line=-5/3., w_edge=0, gamma=0.1)
assert_equal(np.array(snake[:10, 0], dtype=np.int32), refx)
assert_equal(np.array(snake[:10, 1], dtype=np.int32), refy)
@skipif(not new_scipy)
def test_end_points():
img = data.astronaut()
img = rgb2gray(img)
s = np.linspace(0, 2*np.pi, 400)
x = 220 + 100*np.cos(s)
y = 100 + 100*np.sin(s)
init = np.array([x, y]).T
snake = active_contour(gaussian_filter(img, 3), init,
bc='periodic', alpha=0.015, beta=10, w_line=0, w_edge=1,
gamma=0.001, max_iterations=100)
assert np.sum(np.abs(snake[0, :]-snake[-1, :])) < 2
snake = active_contour(gaussian_filter(img, 3), init,
bc='free', alpha=0.015, beta=10, w_line=0, w_edge=1,
gamma=0.001, max_iterations=100)
assert np.sum(np.abs(snake[0, :]-snake[-1, :])) > 2
snake = active_contour(gaussian_filter(img, 3), init,
bc='fixed', alpha=0.015, beta=10, w_line=0, w_edge=1,
gamma=0.001, max_iterations=100)
assert_allclose(snake[0, :], [x[0], y[0]], atol=1e-5)
@skipif(not new_scipy)
def test_bad_input():
img = np.zeros((10, 10))
x = np.linspace(5, 424, 100)
y = np.linspace(136, 50, 100)
init = np.array([x, y]).T
assert_raises(ValueError, active_contour, img, init,
bc='wrong')
assert_raises(ValueError, active_contour, img, init,
max_iterations=-15)
if __name__ == "__main__":
np.testing.run_module_suite()