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]:
200

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 detection in medical images.

Load an image $f$.

In [34]:
n = 256
name = "NtToolBox/src/data/cortex.png"
f = load_image(name, n);
Out[34]:
256×256 Array{Float32,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 ⋮                        ⋮              ⋱       ⋮                        ⋮  
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0

Display.

In [35]:
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.

In [36]:
G = grad(f)
d0 = sqrt(sum(G.^2, 3))
imageplot(d0[:, :, 1])

Blur it by $h_a$.

In [37]:
a = 2
d = gaussian_blur(d0, a)
imageplot(d[:, :, 1])
In [95]:
f = copy(d0)

n = maximum(size(f))
# t = [collect(0:n/2); collect(-n/2:-2)]
# (Y, X) = meshgrid(t, t)
# h = exp( -(X.^2 + Y.^2)./(2.0*float(sigma)^2) )
# h = h./sum(h)
# real(plan_ifft((plan_fft(f)*f).*(plan_fft(h)*h))*((plan_fft(f)*f).*(plan_fft(h)*h)))
Out[95]:
128.0

Compute a decreasing function of the gradient to define $W$.

In [38]:
d = min(d, .4)
W = rescale(-d, .8, 1);
Out[38]:
256×256×1 Array{Float64,3}:
[:, :, 1] =
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  …  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  …  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  …  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 ⋮                        ⋮              ⋱       ⋮                        ⋮  
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  …  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  …  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  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 it.

In [39]:
imageplot(W[:, :, 1])

Number of points.

In [40]:
p = 128;

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

In [57]:
include("NtSolutions/segmentation_2_snakes_param/exo3.jl")
Out[57]:
(-0.5,255.5,255.5,-0.5)

Step size.

In [58]:
dt = 2;

Number of iterations.

In [59]:
Tmax = 9000
niter = round(Tmax/ dt);
Out[59]:
4500.0

Exercise 4: Perform the curve evolution.

In [61]:
include("NtSolutions/segmentation_2_snakes_param/exo4.jl")

# G = grad(W);
# G = G[:, :, 1] + 1im*G[:, :, 2]
# EvalG = gamma -> NtToolBox.bilinear_interpolate(G, imag(gamma), real(gamma))
# EvalW = gamma -> NtToolBox.bilinear_interpolate(W, imag(gamma), real(gamma))
# #
# gamma = gamma0
# displist = round(linspace(0, niter, 10))
# k = 1;
# clf
# imageplot(transpose(f))
# for i in 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) | (i == niter)
#             lw = 4;
#         end
#         cplot(gamma, "r", lw);
#         k = k + 1;
#         axis("equal"); axis("off");
#     end
# end

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

In [125]:
n = 256
f = load_image(name, n)
f = f[45:105, 60:120]
n = size(f)[1];
Out[125]:
61

Display.

In [126]:
imageplot(f)

Exercise 5: Compute an edge attracting criterion $W(x)>0$, that is small in area of strong gradient.

In [128]:
include("NtSolutions/segmentation_2_snakes_param/exo5.jl")

Start and end points $x_0$ and $x_1$.

In [129]:
x0 = 4 + 55*1im
x1 = 53 + 4*1im;
Out[129]:
53 + 4im

Initial curve $\ga_0$.

In [154]:
p = 128
t = transpose(linspace(0, 1, p))
gamma0 = t*x1 + (1 - t)*x0;
Out[154]:
1×128 Array{Complex{Float64},2}:
 4.0+55.0im  4.38583+54.5984im  …  52.6142+4.40157im  53.0+4.0im

Initialize the evolution.

In [138]:
gamma = gamma0;
Out[138]:
128-element Array{Complex{Float64},1}:
     4.0+55.0im   
 4.38583+54.5984im
 4.77165+54.1969im
 5.15748+53.7953im
 5.54331+53.3937im
 5.92913+52.9921im
 6.31496+52.5906im
 6.70079+52.189im 
 7.08661+51.7874im
 7.47244+51.3858im
 7.85827+50.9843im
 8.24409+50.5827im
 8.62992+50.1811im
        ⋮         
 48.7559+8.41732im
 49.1417+8.01575im
 49.5276+7.61417im
 49.9134+7.2126im 
 50.2992+6.81102im
  50.685+6.40945im
 51.0709+6.00787im
 51.4567+5.6063im 
 51.8425+5.20472im
 52.2283+4.80315im
 52.6142+4.40157im
    53.0+4.0im    

Display.

In [140]:
clf
imageplot(transpose(W))
cplot(gamma[:],"r", 2)
plot(real(gamma[1]), imag(gamma[1]), "b.", markersize = 20)
plot(real(gamma[end]), imag(gamma[end]), "b.", markersize = 20);

Re-sampling for non-periodic curves.

In [141]:
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]] ) )
Out[141]:
(::#57) (generic function with 1 method)

Time step.

In [142]:
dt = 1/10;
Out[142]:
0.1

Number of iterations.

In [143]:
Tmax = 2000*4/ 7
niter = round(Tmax/ dt);
Out[143]:
11429.0

Exercise 6: Perform the curve evolution. Be careful to impose the boundary conditions at each step.

In [156]:
include("NtSolutions/segmentation_2_snakes_param/exo6.jl")