Active Contours using Parameteric Curves¶

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}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$

Important: Please read the installation page for details about how to install the toolboxes.

In [2]:

using PyPlot
using NtToolBox
# using Autoreload
# arequire("NtToolBox")

Parameteric 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.

In [4]:

gamma0 = [.78, .14, .42, .18, .32, .16, .75, .83, .57, .68, .46, .40, .72, .79, .91, .90] + 1im*[.87, .82, .75, .63, .34, .17, .08, .46, .50, .25, .27, .57, .73, .57, .75, .79];

Display the initial curve.

In [5]:

periodize = gamma -> [gamma; [gamma[1]]]
function cplot(gamma, s = "b", lw = 1) 
    plot(real(periodize(gamma)), imag(periodize(gamma)), s, linewidth = lw)
    axis("equal")
    axis("off")
end
cplot(gamma0, "b.-")

Out[5]:

(0.1015,0.9485,0.040499999999999994,0.9095)

Number of points of the discrete curve.

In [6]:

p = 256;

Shortcut to re-sample a curve according to arc length.

In [7]:

using Dierckx
curvabs = gamma -> [0; cumsum(1e-5 .+ abs(gamma[1:end - 1] - gamma[2:end]) )]
resample11 = (gamma, d) -> evaluate(Spline1D(d./d[end], real(gamma), k=1), (0:p-1)./p)
resample12 = (gamma, d) -> evaluate(Spline1D(d./d[end], imag(gamma), k=1), (0:p-1)./p)
resample = gamma -> resample11( [gamma; gamma[1]], curvabs( [gamma; gamma[1]] ) ) + 1im*resample12( [gamma; gamma[1]], curvabs( [gamma; gamma[1]] ) );

ArgumentError: Module Dierckx not found in current path.
Run `Pkg.add("Dierckx")` to install the Dierckx package.

 in require(::Symbol) at ./loading.jl:365
 in require(::Symbol) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?

Initial curve $ \ga_1(t)$.

In [6]:

gamma1 = resample(gamma0);

Out[6]:

256-element Array{Complex{Float64},1}:
     0.78+0.87im    
 0.761553+0.868559im
 0.743106+0.867118im
 0.724659+0.865676im
 0.706212+0.864235im
 0.687764+0.862794im
 0.669317+0.861353im
  0.65087+0.859912im
 0.632423+0.858471im
 0.613976+0.857029im
 0.595529+0.855588im
 0.577082+0.854147im
 0.558635+0.852706im
         ⋮          
 0.889715+0.719573im
 0.899979+0.734968im
 0.909894+0.750424im
 0.905407+0.768371im
 0.900921+0.786317im
 0.887764+0.798157im
 0.872369+0.808421im
 0.856974+0.818684im
 0.841579+0.828947im
 0.826185+0.83921im 
  0.81079+0.849474im
 0.795395+0.859737im

Display the initial curve.

In [7]:

cplot(gamma1, "k");

Out[7]:

(0.10727700934498144,0.9481138973125205,0.040705710689687556,0.909490204252872)

Shortcut for forward and backward finite differences.

In [8]:

BwdDiff = c -> c - [c[end]; c[1:end - 1]]
FwdDiff = c -> [c[2:end]; c[1]] - c
dotp = (c1, c2) -> real(c1.*conj(c2));

Out[8]:

(::#15) (generic function with 1 method)

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.

In [9]:

normalize = v -> v./max(abs(v), 1e-10)
tangent = gamma -> normalize( FwdDiff(gamma) )
normal = gamma -> -1im*tangent(gamma);

Out[9]:

(::#21) (generic function with 1 method)

Move the curve in the normal direction, by computing $ \ga_1(t) \pm \delta n_{\ga_1}(t) $.

In [10]:

delta = .03
gamma2 = gamma1 .+ delta .* normal(gamma1)
gamma3 = gamma1 .- delta .* normal(gamma1);

Out[10]:

256-element Array{Complex{Float64},1}:
 0.782337+0.840091im
  0.76389+0.83865im 
 0.745442+0.837209im
 0.726995+0.835768im
 0.708548+0.834326im
 0.690101+0.832885im
 0.671654+0.831444im
 0.653207+0.830003im
  0.63476+0.828562im
 0.616313+0.827121im
 0.597866+0.825679im
 0.579418+0.824238im
 0.560971+0.822797im
         ⋮          
 0.864754+0.736214im
 0.874728+0.751167im
  0.88079+0.743148im
 0.876303+0.761095im
 0.880853+0.764018im
 0.871123+0.773196im
 0.855728+0.783459im
 0.840333+0.793722im
 0.824938+0.803986im
 0.809544+0.814249im
 0.794149+0.824512im
 0.778754+0.834775im

Display the curves.

In [11]:

cplot(gamma1, "k")
cplot(gamma2, "r--")
cplot(gamma3, "b--")
axis("tight") 
axis("off")

Out[11]:

(0.09818192340611892,0.979037189504653,0.011000256259768873,0.942237845572892)

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) $.

In [12]:

normalC = gamma -> BwdDiff(tangent(gamma)) ./ abs( FwdDiff(gamma) );

Out[12]:

(::#23) (generic function with 1 method)

Time step for the evolution. It should be very small because we use an explicit time stepping and the curve has strong curvature.

In [13]:

dt = 0.001 / 100;

Out[13]:

1.0e-5

Number of iterations.

In [14]:

Tmax = 3.0 / 100
niter = round(Tmax/dt);

Out[14]:

3000.0

Initialize the curve for $s=0$.

In [15]:

gamma = gamma1;

Out[15]:

256-element Array{Complex{Float64},1}:
     0.78+0.87im    
 0.761553+0.868559im
 0.743106+0.867118im
 0.724659+0.865676im
 0.706212+0.864235im
 0.687764+0.862794im
 0.669317+0.861353im
  0.65087+0.859912im
 0.632423+0.858471im
 0.613976+0.857029im
 0.595529+0.855588im
 0.577082+0.854147im
 0.558635+0.852706im
         ⋮          
 0.889715+0.719573im
 0.899979+0.734968im
 0.909894+0.750424im
 0.905407+0.768371im
 0.900921+0.786317im
 0.887764+0.798157im
 0.872369+0.808421im
 0.856974+0.818684im
 0.841579+0.828947im
 0.826185+0.83921im 
  0.81079+0.849474im
 0.795395+0.859737im

Evolution of the curve.

In [16]:

gamma = gamma .+ dt .* normalC(gamma);

Out[16]:

256-element Array{Complex{Float64},1}:
 0.779911+0.869658im
 0.761553+0.868559im
 0.743106+0.867118im
 0.724659+0.865676im
 0.706212+0.864235im
 0.687764+0.862794im
 0.669317+0.861353im
  0.65087+0.859912im
 0.632423+0.858471im
 0.613976+0.857029im
 0.595529+0.855588im
 0.577082+0.854147im
 0.558635+0.852706im
         ⋮          
 0.889715+0.719573im
 0.899971+0.734974im
 0.909471+0.750493im
 0.905407+0.768371im
 0.900638+0.786147im
 0.887716+0.798096im
 0.872369+0.808421im
 0.856974+0.818684im
 0.841579+0.828947im
 0.826185+0.83921im 
  0.81079+0.849474im
 0.795395+0.859737im

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).

In [17]:

gamma = resample(gamma);

Out[17]:

256-element Array{Complex{Float64},1}:
 0.779911+0.869658im
  0.76173+0.868569im
 0.743572+0.867154im
 0.725414+0.865735im
 0.707255+0.864317im
 0.689097+0.862898im
 0.670939+0.86148im 
 0.652781+0.860061im
 0.634623+0.858642im
 0.616465+0.857224im
 0.598307+0.855805im
 0.580149+0.854387im
  0.56199+0.852968im
         ⋮          
 0.890755+0.721134im
 0.900799+0.736326im
 0.909116+0.752056im
 0.905023+0.769802im
  0.89969+0.787023im
 0.886136+0.799159im
 0.871021+0.80932im 
 0.855866+0.819423im
 0.840712+0.829526im
 0.825557+0.839629im
 0.810403+0.849732im
 0.795246+0.859832im

Exercise 1: Perform the curve evolution. You need to resample it a few times.

In [18]:

include("NtSolutions/segmentation_2_snakes_param/exo1.jl")

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$.

In [19]:

n = 200;

Out[19]:

Create a synthetic weight $W(x)$.

In [20]:

nbumps = 40
theta = rand(nbumps,1).*(2*pi)
r = .6*n/2
a = [.62*n,.6*n]
x = round( a[1] + r*cos(theta) )
y = round( a[2] + r*sin(theta) )
W = zeros(n, n)
for i in 1:nbumps
    W[Int(x[i]), Int(y[i])] = 1
end
W = gaussian_blur(W, 6.0)
W = rescale( -min(W, .05), .3,1);

Out[20]:

200×200 Array{Float64,2}:
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  …  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  …  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  …  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 ⋮                        ⋮              ⋱            ⋮                      
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  …  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  …  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0

Display the metric $W$.

In [21]:

imageplot(W)

Pre-compute the gradient $\nabla W(x)$ of the metric.

In [22]:

G = grad(W)
G = G[:,:,1] .+ 1im*G[:, :, 2];

Out[22]:

200×200 Array{Complex{Float64},2}:
                  0.0+0.0im  …                   0.0+0.0im
                  0.0+0.0im                      0.0+0.0im
                  0.0+0.0im                      0.0+0.0im
                  0.0+0.0im                      0.0+0.0im
                  0.0+0.0im                      0.0+0.0im
                  0.0+0.0im  …                   0.0+0.0im
                  0.0+0.0im                      0.0+0.0im
                  0.0+0.0im                      0.0+0.0im
                  0.0+0.0im                      0.0+0.0im
                  0.0+0.0im                      0.0+0.0im
                  0.0+0.0im  …                   0.0+0.0im
                  0.0+0.0im                      0.0+0.0im
                  0.0+0.0im                      0.0+0.0im
                     ⋮       ⋱                            
 5.73319e-13+6.05738e-13im      1.39155e-12+1.44751e-12im 
 2.52687e-13+2.62457e-13im      6.15064e-13+6.29274e-13im 
  1.0858e-13+1.11022e-13im   …  2.64899e-13+2.66898e-13im 
 4.57412e-14+4.59632e-14im      1.11466e-13+1.10578e-13im 
 1.84297e-14+1.82077e-14im       4.55191e-14+4.4853e-14im 
 7.32747e-15+7.32747e-15im      1.82077e-14+1.77636e-14im 
 2.88658e-15+2.88658e-15im      7.10543e-15+6.88338e-15im 
 1.11022e-15+1.11022e-15im   …  2.88658e-15+2.66454e-15im 
 4.44089e-16+4.44089e-16im      8.88178e-16+8.88178e-16im 
 2.22045e-16+2.22045e-16im      4.44089e-16+4.44089e-16im 
                  0.0+0.0im     2.22045e-16+2.22045e-16im 
                  0.0+0.0im                      0.0+0.0im

Display the image of the magnitude $\norm{\nabla W(x)}$ of the gradient.

In [23]:

imageplot(abs(G))

Shortcut to evaluate the gradient and the potential along a curve.

In [24]:

EvalG = gamma -> NtToolBox.bilinear_interpolate(G, imag(gamma), real(gamma))
EvalW = gamma -> NtToolBox.bilinear_interpolate(W, imag(gamma), real(gamma));

Out[24]:

(::#27) (generic function with 1 method)

Create a circular curve $\ga_0$.

In [153]:

r = .98*n/2 # radius
p = 128 # number of points on the curve
theta = transpose( linspace(0, 2*pi, p + 1) )
theta = theta[1 : end - 1]
gamma0 = n/2 * (1 + 1im) +  r*(cos(theta) + 1im*sin(theta));

Initialize the curve at time $t=0$ with a circle.

In [26]:

gamma = gamma0;

For this experiment, the time step should be larger, because the curve is in $[0,n-1] \times [0,n-1]$.

In [27]:

dt = 1;

Number of iterations.

In [28]:

Tmax = 5000
niter = round(Tmax/ dt);

Display the curve on the background.

In [29]:

lw = 2
clf
imageplot(transpose(W))
cplot(gamma, "r", lw)

Out[29]:

(-0.5,199.5,199.5,-0.5)

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.

In [30]:

dotp = (c1,c2) -> real(c1.*conj(c2));

Out[30]:

(::#29) (generic function with 1 method)

Evolution of the curve according to this gradient.

In [31]:

N = normal(gamma)
g = - EvalW(gamma) .* normalC(gamma) + dotp(EvalG(gamma), N) .* N
gamma = (gamma - dt*g);

Out[31]:

128-element Array{Complex{Float64},1}:
  197.99+100.0im  
 197.872+104.808im
 197.518+109.605im
 196.929+114.378im
 196.107+119.117im
 195.053+123.81im 
  193.77+128.445im
 192.262+133.012im
 190.531+137.499im
 188.582+141.896im
 186.419+146.192im
 184.049+150.377im
 181.475+154.44im 
        ⋮         
 181.476+45.5598im
 184.049+49.6232im
 186.419+53.8079im
 188.582+58.104im 
 190.531+62.5009im
 192.262+66.9882im
  193.77+71.5551im
 195.053+76.1904im
 196.107+80.8831im
 196.929+85.6219im
 197.518+90.3953im
 197.872+95.1919im

To avoid the curve from being poorly sampled, it is important to re-sample it evenly.

In [32]:

gamma = resample(gamma);

Out[32]:

128-element Array{Complex{Float64},1}:
  197.99+100.0im  
 197.872+104.808im
 197.518+109.605im
 196.929+114.378im
 196.107+119.117im
 195.053+123.81im 
  193.77+128.445im
 192.262+133.012im
 190.531+137.499im
 188.582+141.896im
 186.419+146.192im
 184.049+150.377im
 181.475+154.44im 
        ⋮         
 181.476+45.5598im
 184.049+49.6232im
 186.419+53.8079im
 188.582+58.104im 
 190.531+62.5009im
 192.262+66.9882im
  193.77+71.5551im
 195.053+76.1904im
 196.107+80.8831im
 196.929+85.6219im
 197.518+90.3953im
 197.872+95.1919im

Exercise 2: Perform the curve evolution.

In [33]:

include("NtSolutions/segmentation_2_snakes_param/exo2.jl")


# gamma = copy(gamma0);
# displist = round(linspace(0, niter, 10))
# k = 1;
# clf; 
# imageplot(transpose(W));
# track  = gamma
# for i in 0 : niter
#     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) | (i == niter)
#             lw = 4
#         end
#         cplot(gamma, "r", lw);
#         k = k + 1;
#         axis("equal"); axis("off");
#     end
# end

Medical Image Segmentation¶

One can use a gradient-based metric to perform edge