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{\Si}{\Sigma} \newcommand{\be}{\beta} \newcommand{\de}{\delta} \newcommand{\De}{\Delta} \renewcommand{\phi}{\varphi} \renewcommand{\th}{\theta} \newcommand{\om}{\omega} \newcommand{\Om}{\Omega} $
Important: You need to download the file nt_toolbox.py
from the
root of the github repository.
from nt_toolbox import *
%load_ext autoreload
%autoreload 2
The autoreload extension is already loaded. To reload it, use: %reload_ext autoreload
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([.78, .14, .42, .18, .32, .16, .75, .83, .57, .68, .46, .40, .72, .79, .91, .90]) + 1j*array([.87, .82, .75, .63, .34, .17, .08, .46, .50, .25, .27, .57, .73, .57, .75, .79])
Display the initial curve.
def periodize(gamma): return 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.
def interpc(x,xf,yf): return interp(x,xf,real(yf)) + 1j * interp(x,xf,imag(yf));
def curvabs(gamma): return concatenate( ([0], cumsum( 1e-5 + abs(gamma[:-1:]-gamma[1::]) ) ) );
def resample1(gamma,d): return interpc(arange(0,p)/float(p), d/d[-1],gamma);
def resample(gamma): return 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.
def shiftR(c): return concatenate( ([c[-1]],c[:-1:]) );
def shiftL(c): return concatenate( (c[1::],[c[0]]) );
def BwdDiff(c): return c - shiftR(c);
def FwdDiff(c): return shiftL(c) - c;
def dotp(c1,c2): return dot(real(c1),real(c2)) + dot(imag(c1),imag(c2));
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.
def normalize(v): return v/maximum(abs(v),1e-10)
def tangent(gamma): return normalize( FwdDiff(gamma) )
def normal(gamma): return -1j*tangent(gamma)
Move the curve in the normal direction, by computing $ \ga_1(t) \pm \delta n_{\ga_1}(t) $.
delta = .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');
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) $.
def normalC(gamma): return 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);
Perform the curve evolution. You need to 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 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.
Create a synthetic weight $W(x)$.
n = 200;
nbumps = 40;
theta = random.rand(nbumps,1)*2*pi;
r = .6*n/2;
a = array([.62*n,.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]),int(y[i])] = 1;
W = gaussian_blur(W,6.0);
W = rescale( -minimum(W,.05), .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))
def EvalG(gamma): interp(1:n,1:n, G, imag(gamma), real(gamma)); def EvalW(gamma): interp(1:n,1:n, W, imag(gamma), real(gamma));
x = linspace(0,10,5)
y = linspace(0,6,5)
x = [0.5, 2]
y = [0.5, 1]
interp2([],[],f,x,y)
print x,y
print [x,y]
--------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) <ipython-input-201-b8baedde3359> in <module>() 3 x = [0.5, 2] 4 y = [0.5, 1] ----> 5 interp2([],[],f,x,y) 6 print x,y 7 print [x,y] /Users/gpeyre/Dropbox/github/numerical-tours/python/nt_toolbox.py in interp2(x, y, f, xi, yi, k) 275 interp2(1:n,1:n, G, imag(gamma), real(gamma)); 276 """ --> 277 return ndimage.map_coordinates(f, [x, y], order=k); 278 279 def grad(f): /Users/gpeyre/anaconda/lib/python2.7/site-packages/scipy/ndimage/interpolation.pyc in map_coordinates(input, coordinates, output, order, mode, cval, prefilter) 295 output_shape = coordinates.shape[1:] 296 if input.ndim < 1 or len(output_shape) < 1: --> 297 raise RuntimeError('input and output rank must be > 0') 298 if coordinates.shape[0] != input.ndim: 299 raise RuntimeError('invalid shape for coordinate array') RuntimeError: input and output rank must be > 0