Active Contours using Level Sets¶

Important: Please read the installation page for details about how to install the toolboxes. $\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\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}$

This tour explores image segementation using level set methods.

In [2]:

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

Managing level set functions¶

In the level set formalism, the evolution of some curve $ (\ga(t))_{t=0}^1 $ is computed by evolving the zero level of a function $\phi : \RR^2 \rightarrow \RR $ $$ \enscond{\ga(s)}{ s \in [0,1] } = \enscond{x \in \RR^2}{\phi(x)=0}. $$ This corresponds to replacing the parameteric representation $\ga$ of the curve by an implicit representation. This requires an additional dimension (and hence more storage) but ease the handling of topological change of the curve during the evolution.

Discretazion size $n \times n$ of the domain $[0,1]^2$.

In [3]:

# include("NtToolBox/src/ndgrid.jl")
n = 200
(Y, X) = meshgrid(collect(1:n), collect(1:n));

One can create a circular shape by using the signed distance function to a circle $$ \phi_1(x) = \sqrt{ (x_1-c_1)^2 + (x_2-c_2)^2 } - r $$ where $r>0$ is the radius and $c \in \RR^2$ the center.

Radius $r$.

In [4]:

r = n/3.;

Center $c$.

In [5]:

c = [r, r] .+ 10;

Distance function $\phi_1$.

In [6]:

phi1 = sqrt((X .- c[1]).^2 + (Y .- c[2]).^2) - r;

Exercise 1

Load a square shape $\phi_2$ at a different position for the center.

In [7]:

include("NtSolutions/segmentation_3_snakes_levelset/exo1.jl");

In [7]:

## Insert your code here.

Display the curves associated to $\phi_1$ and $\phi_2$.

In [7]:

figure(figsize = (10,5))

subplot(1,2,1)
NtToolBox.plot_levelset(phi1)

subplot(1,2,2)
NtToolBox.plot_levelset(phi2)

Exercise 2

Compute the intersection and the union of the two shapes. Store the union in $\phi_0$ (phi0) that we will use in the remaining part of the tour.

In [8]:

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

In [10]:

## Insert your code here.

Mean Curvature Motion.¶

The mean curvature motion corresponds to the minimizing flow of the length of the curve $$ \int_0^1 \norm{\ga'(s)} d s. $$

It is implemeted in a level set formalism by a familly $\phi_t$ of level set function parameterized by an artificial time $t \geq 0$, that satisfies the following PDE $$ \pd{\phi_t}{t} = -G(\phi_t) \qwhereq G(\phi) = -\norm{\nabla \phi} \text{div} \pa{ \frac{\nabla \phi}{\norm{\nabla \phi}} } $$ and where $\nabla \phi_t(x) \in \RR^2$ is the spacial gradient.

This flow is computed using a gradient descent $\phi^{(0)} = \phi_0$ and $$ \phi^{(\ell+1)} = \phi^{(\ell)} - \tau G(\phi^{(\ell)}), $$ where $\tau>0$ is small enough time step.

Maximum time of the evolution $0 \leq t \leq t_{\max}$.

In [9]:

Tmax = 200;

Time step $\tau>0$ (should be small).

In [10]:

tau = .5;

Number of iterations.

In [11]:

niter = Int(Tmax/tau);

Initial shape $\phi^{(0)}$ at $t=0$.

In [12]:

phi = copy(phi0);

We now compute the right hand side of the evolution equation.

Compute the gradient $\nabla \phi$. We use centered differences for the discretization of the gradient.

In [13]:

g0 = NtToolBox.Grad(phi, "sym", 2);

Norm $\norm{\nabla \phi}$ of the gradient.

In [14]:

d = max(eps().*ones(n, n), sqrt(sum(g0.^2, 3)));

Normalized gradient.

In [15]:

g = g0./cat(3, d[:, :], d[:, :]);

The curvature term.

In [16]:

K = - d.*NtToolBox.Div(g[:, :, 1], g[:, :, 2], "sym", 2);

Perform one step of the gradient descent.

In [17]:

phi = phi - tau.*K;

Exercise 3

Implement the mean curvature motion.

In [18]:

include("NtSolutions/segmentation_3_snakes_levelset/exo3.jl")

In [21]:

## Insert your code here.

Levelset Re-distancing¶

During PDE resolution, a level set function $\phi$ might become ill-conditionned, so that the zero crossing is not sharp enough. The quality of the level set function is restored by computing the signed distance function to the zero level set.

This corresponds to first extracting the zero level set $$ \Cc = \enscond{x \in \RR^2 }{\phi(x)=0}, $$ and then solving the following eikonal equation PDE on $\tilde \phi$ (in viscosity sense) $$ \norm{\nabla \tilde \phi(x)} = 1 \qandq \forall y \in \Cc, \tilde\phi(y)=0. $$ The one can replace $\phi$ by $\text{sign}(\phi(x))\tilde \phi(x)$ which is the signed distance function to $\Cc$.

We set $\phi=\phi_0^3$ so that they are both valid level set function of the same curve, but $\phi$ is not the signed distance function.

In [19]:

phi = phi0.^3;

Solve the eikonal PDE using the Fast Marching algorithm. You have to install a C++ compiler (https://wiki.python.org/moin/WindowsCompilers#Microsoft_Visual_C.2B-.2B-_14.0_standalone:_Visual_C.2B-.2B-_Build_Tools_2015_.28x86.2C_x64.2C_ARM.29) and the package scikit-fmm (skfmm) to run this function (pip install scikit_fmm in the console). You also need to have PyCall package installed in Julia.

In [20]:

include("NtToolBox/src/perform_fast_marching.jl")
include("NtToolBox/src/perform_redistancing.jl")

phi1 = perform_redistancing(phi0);

Display the level sets.

In [21]:

figure(figsize = (10, 5))

subplot(1, 2, 1)
plot_levelset(phi)
title("Before redistancing")

subplot(1, 2, 2)
plot_levelset(phi1)
title("After redistancing")

show()

Edge-based Segmentation with Geodesic Active Contour¶

Geodesic active contours compute loval minimum of a weighted geodesic distance that attract the curve toward the features of the background image.

Note: these active contours should not be confounded with the geodesic shortest paths, that are globally minimizing geodesics between two points. Here the active contour is a close curve progressively decreasing a weighted geodesic length that is only a local minimum (the global minimum would be a single point).

Size of the image.

In [22]:

n = 200;

First we load an image $f_0 \in \RR^{n \times n}$ to segment.

In [23]:

f0 = NtToolBox.rescale(NtToolBox.load_image("NtToolBox/src/data/cortex.png", n));

Given a background image $f_0$ to segment, one needs to compute an edge-stopping function $W$. It should be small in area of high gradient, and high in area of large gradient.

We use here $$ W(x) = \al + \frac{\be}{\epsilon + d(x) } \qwhereq d = \norm{\nabla f_0} \star h_a, $$ and where $h_a$ is a blurring kernel of size $a>0$.

Compute the magnitude of the gradient $d_0(x) = \norm{\nabla f_0(x)}$.

In [24]:

g = NtToolBox.Grad(f0, "sym", 2)
d0 = sqrt(sum(g.^2, 3));

Blur size $a$.

In [25]:

a = 5;

Compute the blurring $d = d_0 \star h_a$.

In [26]:

d = NtToolBox.perform_blurring(d0[:, :, 1], [a], "per");

Parameter $\epsilon>0$.

In [27]:

epsilon = 1e-1;

We set the $\al$ and $\be$ parameters to adjust the overall values of $W$ (equivalently we use the function rescale).

In [28]:

W = 1./(epsilon + d)
W = NtToolBox.rescale(-d, 0.1, 1);

Display it.

In [29]:

figure(figsize = (10, 5))
NtToolBox.imageplot(f0, "Image to segment", [1, 2, 1])
NtToolBox.imageplot(W, "Weight", [1, 2, 2])