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.
using PyPlot
using NtToolBox
# using Autoreload
# arequire("NtToolBox")
WARNING: Method definition ndgrid(AbstractArray{T<:Any, 1}) in module NtToolBox at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:3 overwritten at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:3. WARNING: Method definition ndgrid(AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}) in module NtToolBox at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:6 overwritten at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:6. WARNING: Method definition ndgrid_fill(Any, Any, Any, Any) in module NtToolBox at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:13 overwritten at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:13. WARNING: Method definition ndgrid(AbstractArray{#T<:Any, 1}...) in module NtToolBox at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:19 overwritten at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:19. WARNING: Method definition meshgrid(AbstractArray{T<:Any, 1}) in module NtToolBox at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:33 overwritten at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:33. WARNING: Method definition meshgrid(AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}) in module NtToolBox at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:36 overwritten at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:36. WARNING: Method definition meshgrid(AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}) in module NtToolBox at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:44 overwritten at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:44.
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 = [.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.
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.-")
(0.1015,0.9485,0.040499999999999994,0.9095)
Number of points of the discrete curve.
p = 256;
Shortcut to re-sample a curve according to arc length.
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)$.
gamma1 = resample(gamma0);
UndefVarError: resample not defined
Display the initial curve.
cplot(gamma1, "k");
UndefVarError: gamma1 not defined
Shortcut for forward and backward finite differences.
BwdDiff = c -> c - [c[end]; c[1:end - 1]]
FwdDiff = c -> [c[2:end]; c[1]] - c
dotp = (c1, c2) -> real(c1.*conj(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.
normalize = v -> v./max(abs(v), 1e-10)
tangent = gamma -> normalize( FwdDiff(gamma) )
normal = gamma -> -1im*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);
UndefVarError: gamma1 not defined
Display the curves.
cplot(gamma1, "k")
cplot(gamma2, "r--")
cplot(gamma3, "b--")
axis("tight")
axis("off")
UndefVarError: gamma1 not defined
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 = 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;
UndefVarError: gamma1 not defined
Evolution of the curve.
gamma = gamma .+ dt .* normalC(gamma);
MethodError: no method matching endof(::Base.Math.#gamma) Closest candidates are: endof(::SimpleVector) at essentials.jl:169 endof(::Char) at char.jl:18 endof(::String) at strings/string.jl:37 ... in (::##5#6)(::Function) at ./In[8]:2 in (::##11#12)(::Function) at ./In[9]:2 in (::##15#16)(::Function) at ./In[12]:1
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);
UndefVarError: resample not defined
Exercise 1: Perform the curve evolution. You need to resample it a few times.
include("NtSolutions/segmentation_2_snakes_param/exo1.jl")
LoadError: UndefVarError: gamma1 not defined while loading /Users/gpeyre/Dropbox/github/numerical-tours/julia/NtSolutions/segmentation_2_snakes_param/exo1.jl, in expression starting on line 1 in include_from_node1(::String) at ./loading.jl:488 in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
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 = 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);
Display the metric $W$.
imageplot(W)
Pre-compute the gradient $\nabla W(x)$ of the metric.
G = grad(W)
G = G[:,:,1] .+ 1im*G[:, :, 2];
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 = gamma -> NtToolBox.bilinear_interpolate(G, imag(gamma), real(gamma))
EvalW = gamma -> NtToolBox.bilinear_interpolate(W, imag(gamma), real(gamma));
Create a circular curve $\ga_0$.
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.
gamma = gamma0;
WARNING: imported binding for gamma overwritten in module Main
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)
(-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.
dotp = (c1,c2) -> real(c1.*conj(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);
MethodError: Cannot `convert` an object of type Array{Float64,1} to an object of type Int64 This may have arisen from a call to the constructor Int64(...), since type constructors fall back to convert methods. in bilinear_interpolate(::Array{Float64,2}, ::Array{Float64,1}, ::Array{Float64,1}) at /Users/gpeyre/.julia/v0.5/NtToolBox/src/signal.jl:414 in (::##19#20)(::Array{Complex{Float64},1}) at ./In[24]:2
To avoid the curve from being poorly sampled, it is important to re-sample it evenly.
gamma = resample(gamma);
UndefVarError: resample not defined
Exercise 2: Perform the curve evolution.
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
LoadError: MethodError: Cannot `convert` an object of type Array{Float64,1} to an object of type Int64 This may have arisen from a call to the constructor Int64(...), since type constructors fall back to convert methods. while loading /Users/gpeyre/Dropbox/github/numerical-tours/julia/NtSolutions/segmentation_2_snakes_param/exo2.jl, in expression starting on line 7 in bilinear_interpolate(::Array{Float64,2}, ::Array{Float64,1}, ::Array{Float64,1}) at /Users/gpeyre/.julia/v0.5/NtToolBox/src/signal.jl:414 in (::##19#20)(::Array{Complex{Float64},1}) at ./In[24]:2 in macro expansion; at /Users/gpeyre/Dropbox/github/numerical-tours/julia/NtSolutions/segmentation_2_snakes_param/exo2.jl:9 [inlined] in anonymous at ./<missing>:? in include_from_node1(::String) at ./loading.jl:488 in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
One can use a gradient-based metric to perform edge detection in medical images.
Load an image $f$.
n = 256
name = "NtToolBox/src/data/cortex.png"
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, 3))
imageplot(d0[:, :, 1])
Blur it by $h_a$.
a = 2
d = gaussian_blur(d0, a)
imageplot(d[:, :, 1])
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)))
256
Compute a decreasing function of the gradient to define $W$.
d = min(d, .4)
W = rescale(-d, .8, 1);
Display it.
imageplot(W[:, :, 1])
Number of points.
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.
include("NtSolutions/segmentation_2_snakes_param/exo3.jl")
LoadError: transpose not implemented for Array{Float32,3}. Consider adding parentheses, e.g. A*(B*C.') instead of A*B*C' to avoid explicit calculation of the transposed matrix. while loading /Users/gpeyre/Dropbox/github/numerical-tours/julia/NtSolutions/segmentation_2_snakes_param/exo3.jl, in expression starting on line 8 in transpose(::Array{Float32,3}) at ./abstractarraymath.jl:10 in include_from_node1(::String) at ./loading.jl:488 in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
Step size.
dt = 2;
Number of iterations.
Tmax = 9000
niter = round(Tmax/ dt);
Exercise 4: Perform the curve evolution.
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
LoadError: transpose not implemented for Array{Float32,3}. Consider adding parentheses, e.g. A*(B*C.') instead of A*B*C' to avoid explicit calculation of the transposed matrix. while loading /Users/gpeyre/Dropbox/github/numerical-tours/julia/NtSolutions/segmentation_2_snakes_param/exo4.jl, in expression starting on line 10 in transpose(::Array{Float32,3}) at ./abstractarraymath.jl:10 in include_from_node1(::String) at ./loading.jl:488 in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
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 = size(f)[1];
Display.
imageplot(f)
Exercise 5: Compute an edge attracting criterion $W(x)>0$, that is small in area of strong gradient.
include("NtSolutions/segmentation_2_snakes_param/exo5.jl")
LoadError: UndefVarError: ff not defined while loading /Users/gpeyre/Dropbox/github/numerical-tours/julia/NtSolutions/segmentation_2_snakes_param/exo5.jl, in expression starting on line 1 in include_from_node1(::String) at ./loading.jl:488 in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
Start and end points $x_0$ and $x_1$.
x0 = 4 + 55*1im
x1 = 53 + 4*1im;
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[1]), imag(gamma[1]), "b.", markersize = 20)
plot(real(gamma[end]), imag(gamma[end]), "b.", markersize = 20);
transpose not implemented for Array{Float64,3}. Consider adding parentheses, e.g. A*(B*C.') instead of A*B*C' to avoid explicit calculation of the transposed matrix. in transpose(::Array{Float64,3}) at ./abstractarraymath.jl:10
Re-sampling for non-periodic curves.
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:?
Time step.
dt = 1/10;
Number of iterations.
Tmax = 2000*4/ 7
niter = round(Tmax/ dt);
Exercise 6: Perform the curve evolution. Be careful to impose the boundary conditions at each step.
include("NtSolutions/segmentation_2_snakes_param/exo6.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
# 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[1] = x0
# gamma[end] = x1
# 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");
# plot(real(gamma[1]), imag(gamma[1]), "b.", markersize = 20)
# plot(real(gamma[end]), imag(gamma[end]), "b.", markersize = 20);
# end
# end
LoadError: MethodError: Cannot `convert` an object of type Array{Float64,1} to an object of type Int64 This may have arisen from a call to the constructor Int64(...), since type constructors fall back to convert methods. while loading /Users/gpeyre/Dropbox/github/numerical-tours/julia/NtSolutions/segmentation_2_snakes_param/exo6.jl, in expression starting on line 11 in bilinear_interpolate(::Array{Float64,3}, ::Array{Float64,1}, ::Array{Float64,1}) at /Users/gpeyre/.julia/v0.5/NtToolBox/src/signal.jl:414 in (::##29#30)(::Array{Complex{Float64},1}) at /Users/gpeyre/Dropbox/github/numerical-tours/julia/NtSolutions/segmentation_2_snakes_param/exo6.jl:4 in macro expansion; at /Users/gpeyre/Dropbox/github/numerical-tours/julia/NtSolutions/segmentation_2_snakes_param/exo6.jl:13 [inlined] in anonymous at ./<missing>:? in include_from_node1(::String) at ./loading.jl:488 in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?