mirror of
https://github.com/wassname/scikit-image.git
synced 2026-06-27 23:39:49 +08:00
Kyle Jackson
235 lines
7.9 KiB
Python
235 lines
7.9 KiB
Python
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
|
|
|
|
Create and smooth image:
|
|
|
|
>>> img = np.zeros((100, 100))
|
|
>>> rr, cc = circle_perimeter(35, 45, 25)
|
|
>>> img[rr, cc] = 1
|
|
>>> img = gaussian(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].astype(np.float), snake[:, 1].astype(np.float)
|
|
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
|