diff --git a/doc/examples/plot_active_contours.py b/doc/examples/plot_active_contours.py
new file mode 100644
index 00000000..a8cf7432
--- /dev/null
+++ b/doc/examples/plot_active_contours.py
@@ -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
+"""
diff --git a/skimage/segmentation/__init__.py b/skimage/segmentation/__init__.py
index f79fb482..63f7d66d 100644
--- a/skimage/segmentation/__init__.py
+++ b/skimage/segmentation/__init__.py
@@ -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',
diff --git a/skimage/segmentation/active_contour_model.py b/skimage/segmentation/active_contour_model.py
new file mode 100644
index 00000000..b2f89b3b
--- /dev/null
+++ b/skimage/segmentation/active_contour_model.py
@@ -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
diff --git a/skimage/segmentation/tests/test_active_contour_model.py b/skimage/segmentation/tests/test_active_contour_model.py
new file mode 100644
index 00000000..48a86d45
--- /dev/null
+++ b/skimage/segmentation/tests/test_active_contour_model.py
@@ -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()