From 320977c8b9d8d1d78852a6f8f2986e06fb0acf68 Mon Sep 17 00:00:00 2001 From: Julius Bier Kirkegaard Date: Tue, 1 Dec 2015 12:42:12 +0000 Subject: [PATCH] Fixed doctest problem --- skimage/segmentation/active_contour_model.py | 466 +++++++++---------- 1 file changed, 233 insertions(+), 233 deletions(-) diff --git a/skimage/segmentation/active_contour_model.py b/skimage/segmentation/active_contour_model.py index 1f5d88ce..43904306 100644 --- a/skimage/segmentation/active_contour_model.py +++ b/skimage/segmentation/active_contour_model.py @@ -1,233 +1,233 @@ -import warnings -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 - >>> int(np.mean(np.sqrt((45-snake[:, 0])**2 + - (35-snake[:, 1])**2))) #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 +import warnings +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