Parametric Active Contours¶
Evolve a closed curve under smoothing and image-dependent forces. Track how curvature and external forces cooperate to attract the curve toward object boundaries, and inspect the role of initialization and step size.
Run this tour¶
Run the cells in order with a Python 3 kernel. The first cell locates the companion data and toolbox and installs missing dependencies when needed. All worked examples include their implementation directly in this notebook. Random seeds make comparisons reproducible; you can change them to explore other samples.
# Locate the companion toolbox locally, or fetch it for a standalone/Colab copy.
from pathlib import Path
import importlib.util
import os
import subprocess
import sys
working = Path.cwd()
candidates = [working, working / "python", working.parent / "python"]
python_dir = next((p for p in candidates if (p / "nt_toolbox").is_dir()), None)
if python_dir is None:
checkout = working / "numerical-tours-support"
if not checkout.exists():
subprocess.run(
[
"git",
"clone",
"--depth",
"1",
"--branch",
"master",
"https://github.com/gpeyre/numerical-tours.git",
str(checkout),
],
check=True,
)
python_dir = checkout / "python"
os.chdir(python_dir)
if str(python_dir) not in sys.path:
sys.path.insert(0, str(python_dir))
requirements = python_dir / "requirements.txt"
if any(
importlib.util.find_spec(name) is None
for name in [
"numpy",
"scipy",
"matplotlib",
"skimage",
"sklearn",
"pywt",
"ipywidgets",
"cvxpy",
"skfmm",
"autograd",
"progressbar",
"celer",
]
):
subprocess.run(
[sys.executable, "-m", "pip", "install", "-r", str(requirements)], check=True
)
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(0)
plt.rcParams.update(
{
"figure.figsize": (8, 4),
"figure.dpi": 100,
"axes.spines.top": False,
"axes.spines.right": False,
"font.size": 11,
"image.cmap": "gray",
}
)
%matplotlib inline
This tour explores image segmentation using parametric active contours. $\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\si}{\sigma}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$
from nt_toolbox.general import rescale
from nt_toolbox.signal import (
bilinear_interpolate,
gaussian_blur,
grad,
imageplot,
load_image,
)
import numpy as np
import matplotlib.pyplot as plt
from numpy import (
abs,
arange,
around,
array,
concatenate,
cos,
cumsum,
equal,
imag,
interp,
linspace,
maximum,
minimum,
number,
pi,
random,
real,
sin,
sqrt,
sum,
transpose,
zeros,
)
from matplotlib.pyplot import axis, clf, matplotlib, plot
%matplotlib inline
Parametric Curves¶
In this tours, the active contours are represented using parametric curve $ \ga : [0,1] \rightarrow \RR^2 $.
This curve is discretized using a piewise linear curve with $p$ segments, and is stored as a complex vector of points in the plane $\ga \in \CC^p$.
Initial polygon.
gamma0 = array(
[
0.78,
0.14,
0.42,
0.18,
0.32,
0.16,
0.75,
0.83,
0.57,
0.68,
0.46,
0.40,
0.72,
0.79,
0.91,
0.90,
]
) + 1j * array(
[
0.87,
0.82,
0.75,
0.63,
0.34,
0.17,
0.08,
0.46,
0.50,
0.25,
0.27,
0.57,
0.73,
0.57,
0.75,
0.79,
]
)
Display the initial curve.
periodize = lambda gamma: concatenate((gamma, [gamma[0]]))
def cplot(gamma, s="b", lw=1):
plot(real(periodize(gamma)), imag(periodize(gamma)), s, linewidth=lw)
axis("equal")
axis("off")
cplot(gamma0, "b.-");
Number of points of the discrete curve.
p = 256
Shortcut to re-sample a curve according to arc length.
interpc = lambda x, xf, yf: interp(x, xf, real(yf)) + 1j * interp(x, xf, imag(yf))
curvabs = lambda gamma: concatenate(([0], cumsum(1e-5 + abs(gamma[:-1:] - gamma[1::]))))
resample1 = lambda gamma, d: interpc(arange(0, p) / float(p), d / d[-1], gamma)
resample = lambda gamma: resample1(periodize(gamma), curvabs(periodize(gamma)))
Initial curve $ \ga_1(t)$.
gamma1 = resample(gamma0)
Display the initial curve.
cplot(gamma1, "k")
Shortcut for forward and backward finite differences.
shiftR = lambda c: concatenate(([c[-1]], c[:-1:]))
shiftL = lambda c: concatenate((c[1::], [c[0]]))
BwdDiff = lambda c: c - shiftR(c)
FwdDiff = lambda c: shiftL(c) - c
The tangent to the curve is computed as $$ t_\ga(s) = \frac{\ga'(t)}{\norm{\ga'(t)}} $$ and the normal is $ n_\ga(t) = t_\ga(t)^\bot. $
Shortcut to compute the tangent and the normal to a curve.
normalize = lambda v: v / maximum(abs(v), 1e-10)
tangent = lambda gamma: normalize(FwdDiff(gamma))
normal = lambda gamma: -1j * tangent(gamma)
Move the curve in the normal direction, by computing $ \ga_1(t) \pm \delta n_{\ga_1}(t) $.
delta = 0.03
gamma2 = gamma1 + delta * normal(gamma1)
gamma3 = gamma1 - delta * normal(gamma1)
Display the curves.
cplot(gamma1, "k")
cplot(gamma2, "r--")
cplot(gamma3, "b--")
axis("tight")
axis("off")
Evolution by Mean Curvature¶
A curve evolution is a series of curves $ s \mapsto \ga_s $ indexed by an evolution parameter $s \geq 0$. The intial curve $\ga_0$ for $s=0$ is evolved, usually by minizing some energy $E(\ga)$ in a gradient descent $$ \frac{\partial \ga_s}{\partial s} = \nabla E(\ga_s). $$
Note that the gradient of an energy is defined with respect to the curve-dependent inner product $$ \dotp{a}{b} = \int_0^1 \dotp{a(t)}{b(t)} \norm{\ga'(t)} d t. $$ The set of curves can thus be thought as being a Riemannian surface.
The simplest evolution is the mean curvature evolution. It corresponds to minimization of the curve length $$ E(\ga) = \int_0^1 \norm{\ga'(t)} d t $$
The gradient of the length is $$ \nabla E(\ga)(t) = -\kappa_\ga(t) n_\ga(t) $$ where $ \kappa_\ga $ is the curvature, defined as $$ \kappa_\ga(t) = \frac{1}{\norm{\ga'(t)}} \dotp{ t_\ga'(t) }{ n_\ga(t) } . $$
Shortcut for normal times curvature $ \kappa_\ga(t) n_\ga(t) $.
normalC = lambda gamma: BwdDiff(tangent(gamma)) / abs(FwdDiff(gamma))
Time step for the evolution. It should be very small because we use an explicit time stepping and the curve has strong curvature.
dt = 0.001 / 100
Number of iterations.
Tmax = 3.0 / 100
niter = round(Tmax / dt)
Initialize the curve for $s=0$.
gamma = gamma1
Evolution of the curve.
gamma = gamma + dt * normalC(gamma)
To stabilize the evolution, it is important to re-sample the curve so that it is unit-speed parametrized. You do not need to do it every time step though (to speed up).
gamma = resample(gamma)
Worked example 1: Perform the curve evolution. We now resample it a few times.
gamma = gamma1
displist = around(linspace(0, niter, 10))
k = 0
for i in arange(0, niter + 1):
gamma = resample(gamma + dt * normalC(gamma))
if i == displist[k]:
lw = 1
if i == 0 or i == niter:
lw = 4
cplot(gamma, "r", lw)
k = k + 1
axis("tight")
axis("off")
Geodesic Active Contours¶
Geodesic active contours minimize a weighted length $$ E(\ga) = \int_0^1 W(\ga(t)) \norm{\ga'(t)} d t, $$ where $W(x)>0$ is the geodesic metric, that should be small in areas where the image should be segmented.
Size of the image $n$.
n = 200
Create a synthetic weight $W(x)$.
nbumps = 40
theta = random.rand(nbumps, 1) * 2 * pi
r = 0.6 * n / 2
a = array([0.62 * n, 0.6 * n])
x = around(a[0] + r * cos(theta))
y = around(a[1] + r * sin(theta))
W = zeros([n, n])
for i in arange(0, nbumps):
W[int(x[i].item()), int(y[i].item())] = 1
W = gaussian_blur(W, 6.0)
W = rescale(-minimum(W, 0.05), 0.3, 1)
Display the metric $W$.
imageplot(W)
Pre-compute the gradient $\nabla W(x)$ of the metric.
G = grad(W)
G = G[:, :, 0] + 1j * G[:, :, 1]
Display the image of the magnitude $\norm{\nabla W(x)}$ of the gradient.
imageplot(abs(G))
Shortcut to evaluate the gradient and the potential along a curve.
EvalG = lambda gamma: bilinear_interpolate(G, imag(gamma), real(gamma))
EvalW = lambda gamma: bilinear_interpolate(W, imag(gamma), real(gamma))
Create a circular curve $\ga_0$.
r = 0.98 * n / 2 # radius
p = 128 # number of points on the curve
theta = transpose(linspace(0, 2 * pi, p + 1))
theta = theta[0:-1]
gamma0 = n / 2 * (1 + 1j) + r * (cos(theta) + 1j * sin(theta))
Initialize the curve at time $t=0$ with a circle.
gamma = gamma0
For this experiment, the time step should be larger, because the curve is in $[0,n-1] \times [0,n-1]$.
dt = 1
Number of iterations.
Tmax = 5000
niter = round(Tmax / dt)
Display the curve on the background.
lw = 2
clf
imageplot(transpose(W))
cplot(gamma, "r", lw)
The gradient of the energy is $$ \nabla E(\ga) = -W(\ga(t)) \kappa_\ga(t) n_\ga(t) + \dotp{\nabla W(\ga(t))}{ n_\ga(t) } n_\ga(t). $$
Pointwise innerproduct on the curve.
dotp = lambda c1, c2: real(c1) * real(c2) + imag(c1) * imag(c2)
Evolution of the curve according to this gradient.
N = normal(gamma)
g = -EvalW(gamma) * normalC(gamma) + dotp(EvalG(gamma), N) * N
gamma = gamma - dt * g
To avoid the curve from being poorly sampled, it is important to re-sample it evenly.
gamma = resample(gamma)
Worked example 2: Perform the curve evolution.
gamma = gamma0
displist = around(linspace(0, niter, 10))
k = 0
clf
imageplot(transpose(W))
for i in arange(0, niter + 1):
N = normal(gamma)
g = EvalW(gamma) * normalC(gamma) - dotp(EvalG(gamma), N) * N
gamma = resample(gamma + dt * g)
if i == displist[k]:
lw = 1
if i == 0 or i == niter:
lw = 4
cplot(gamma, "r", lw)
k = k + 1
axis("equal")
axis("off")
Medical Image Segmentation¶
One can use a gradient-based metric to perform edge detection in medical images.
Load an image $f$.
n = 256
name = "nt_toolbox/data/cortex.bmp"
f = load_image(name, n)
Display.
imageplot(f)
An edge detector metric can be defined as a decreasing function of the gradient magnitude. $$ W(x) = \psi( d \star h_a(x) ) \qwhereq d(x) = \norm{\nabla f(x)}. $$ where $h_a$ is a blurring kernel of width $a>0$.
Compute the magnitude of the gradient.
G = grad(f)
d0 = sqrt(sum(G**2, 2))
imageplot(d0)
Blur it by $h_a$.
a = 2
d = gaussian_blur(d0, a)
imageplot(d)
Compute a decreasing function of the gradient to define $W$.
d = minimum(d, 0.4)
W = rescale(-d, 0.8, 1)
Display it.
imageplot(W)
Number of points.
p = 128
Worked example 3: Create an initial circle $\gamma_0$ of $p$ points. When plotting the image, you need to transpose it to have axis coherent with the cplot.
r = 0.95 * n / 2
p = 128 # number of points on the curve
theta = transpose(linspace(0, 2 * pi, p + 1))
theta = theta[0:-1]
gamma0 = n / 2 * (1 + 1j) + r * (cos(theta) + 1j * sin(theta))
gamma = gamma0
clf
imageplot(transpose(f))
cplot(gamma, "r", 2)
Step size.
dt = 2
Number of iterations.
Tmax = 9000
niter = round(Tmax / dt)
Worked example 4: Perform the curve evolution.
G = grad(W)
G = G[:, :, 0] + 1j * G[:, :, 1]
EvalG = lambda gamma: bilinear_interpolate(G, imag(gamma), real(gamma))
EvalW = lambda gamma: bilinear_interpolate(W, imag(gamma), real(gamma))
#
gamma = gamma0
displist = around(linspace(0, niter, 10))
k = 0
clf
imageplot(transpose(f))
for i in arange(0, niter + 1):
n = normal(gamma)
g = EvalW(gamma) * normalC(gamma) - dotp(EvalG(gamma), n) * n
gamma = resample(gamma + dt * g)
if i == displist[k]:
lw = 1
if i == 0 or i == niter:
lw = 4
cplot(gamma, "r", lw)
k = k + 1
axis("equal")
axis("off")
Evolution of a Non-closed Curve¶
It is possible to perform the evolution of a non-closed curve by adding boundary constraint $$ \ga(0)=x_0 \qandq \ga(1)=x_1. $$
In this case, the algorithm find a local minimizer of the geodesic distance between the two points.
Note that a much more efficient way to solve this problem is to use the Fast Marching algorithm to find the global minimizer of the geodesic length.
Load an image $f$.
n = 256
f = load_image(name, n)
f = f[45:105, 60:120]
n = f.shape[0]
Display.
imageplot(f)
Worked example 5: Compute an edge attracting criterion $W(x)>0$, that is small in area of strong gradient.
G = grad(f)
G = sqrt(sum(G**2, 2))
sigma = 1.5
G = gaussian_blur(G, sigma)
G = minimum(G, 0.4)
W = rescale(-G, 0.4, 1)
clf
imageplot(W)
Start and end points $x_0$ and $x_1$.
x0 = 4 + 55j
x1 = 53 + 4j
Initial curve $\ga_0$.
p = 128
t = transpose(linspace(0, 1, p))
gamma0 = t * x1 + (1 - t) * x0
Initialize the evolution.
gamma = gamma0
Display.
clf
imageplot(transpose(W))
cplot(gamma, "r", 2)
plot(real(gamma[0]), imag(gamma[0]), "b.", markersize=20)
plot(real(gamma[-1]), imag(gamma[-1]), "b.", markersize=20);
Re-sampling for non-periodic curves.
curvabs = lambda gamma: concatenate(([0], cumsum(1e-5 + abs(gamma[:-1:] - gamma[1::]))))
resample1 = lambda gamma, d: interpc(arange(0, p) / float(p - 1), d / d[-1], gamma)
resample = lambda gamma: resample1(gamma, curvabs(gamma))
Time step.
dt = 1 / 10
Number of iterations.
Tmax = 2000 * 4 / 7
niter = round(Tmax / dt)
Worked example 6: Perform the curve evolution. Be careful to impose the boundary conditions at each step.
G = grad(W)
G = G[:, :, 0] + 1j * G[:, :, 1]
EvalG = lambda gamma: bilinear_interpolate(G, imag(gamma), real(gamma))
EvalW = lambda gamma: bilinear_interpolate(W, imag(gamma), real(gamma))
#
gamma = gamma0
displist = around(linspace(0, niter, 10))
k = 0
clf
imageplot(transpose(f))
for i in arange(0, niter + 1):
N = normal(gamma)
g = EvalW(gamma) * normalC(gamma) - dotp(EvalG(gamma), N) * N
gamma = gamma + dt * g
gamma = resample(gamma)
# impose start/end point
gamma[0] = x0
gamma[-1] = x1
if i == displist[k]:
lw = 1
if i == 0 or i == niter:
lw = 4
cplot(gamma, "r", lw)
k = k + 1
axis("equal")
axis("off")
plot(real(gamma[0]), imag(gamma[0]), "b.", markersize=20)
plot(real(gamma[-1]), imag(gamma[-1]), "b.", markersize=20)
References and further reading¶
John Canny. A Computational Approach to Edge Detection. 1986, IEEE Transactions on Pattern Analysis and Machine Intelligence 8(6), 679–698. Detection, localization, and nonmaximum suppression of edges.
Michael Kass, Andrew Witkin, and Demetri Terzopoulos. Snakes: Active Contour Models. 1988, International Journal of Computer Vision 1, 321–331. Parametric curves driven by regularization and image forces.
Stanley Osher and James A. Sethian. Fronts Propagating with Curvature-Dependent Speed: Algorithms Based on Hamilton–Jacobi Formulations. 1988, Journal of Computational Physics 79(1), 12–49. Implicit front propagation using level-set functions.
Tony F. Chan and Luminita A. Vese. Active Contours Without Edges. 2001, IEEE Transactions on Image Processing 10(2), 266–277. Region-based segmentation when image gradients are weak.
Gabriel Peyré. Advanced Signal, Image and Surface Processing. 2010, course notes. A mathematical companion to the Numerical Tours.