Variational image segmentation Nicolas Rougon Institut Mines-Tlcom - - PowerPoint PPT Presentation

Foundations Parameterized active contours Level set methods

Variational image segmentation

Nicolas Rougon

Institut Mines-Télécom / Télécom SudParis ARTEMIS Department; CNRS UMR 8145 nicolas.rougon@telecom-sudparis.eu

May 13, 2020

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods

Outline

1

Foundations

2

Parameterized active contours

3

Level set methods

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods

From edge detection...

◮ Edge detectors are Local: deliver point sets, not curves / surfaces (i.e. manifolds) → post-processing is needed to build segments / patches Blind: detect all image edges → post-processing is needed to filter out irrelevant elements Luminance-based: hence, noise-sensitive and geometry-free → cannot integrate prior shape information

Topology

connectivity (# of parts) homotopy (# of holes) boundary (open / close)

Geometry

smoothness subspace (functions / template, atlas)

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods

...To object segmentation

◮ In many applications, the focus is on object segmentation rather than scene segmentation Expected solutions for object boundaries are manifolds The aim is to detect specific image structures Prior information on object shape (topology / geometry) can be used to constrain the solution space

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods

Active contours

◮ Active contours are deterministic model-based approaches which tackle object segmentation as a shape optimization problem ◮ Solutions are sought as manifolds of a given shape space which achieve the best trade-off between image features (observed data) known object properties (prior model) ◮ This trade-off is expressed via a segmentation criterion combining a data consistency term a prior term ◮ This framework allows for integrating arbitrary features / priors shape color / texture motion

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods

Active contours - A long history

◮ Key milestones Terzopoulos, 1984 parameterized active contours (“snakes”) Caselles et al., 1997 level set-based active contours Chan & Vese, 2001 level set-based active regions Paragios & Deriche, 2002 statistical active regions Bresson et al., 2007 primal-dual level set methods still developing

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods

Active contours - Key ideas

◮ Physically-based modeling An object boundary is modeled as a deformable manifold of the image domain with state described by an energy Prior elastic properties are defined by an internal energy → elastic forces (model cohesion) ≡ shape constraints Interactions with image data are defined by an external energy → image forces (model attraction) ≡ data consistency Starting from an initialization, the model evolves until reaching an equilibrium state with minimal (total) energy ◮ Segmentation is achieved by designing image energy / forces so as to attract the model on relevant image structures energy ≡ segmentation criterion

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods

Geometric modeling

◮ In the sequel, we deal with codimension-1 simply connected, closed submanifolds C over the image domain Ω ⊂ Rn curves in R2 with generic point x = (x, y)T surfaces in R3 with generic point x = (x, y, z)T partitioning Ω into an inner region Rin (≡ object) and an outer region Rout (≡ background) ◮ Active contour models divide into 2 families depending on the geometric representation of C explicit representations lead to parameterized active contours implicit representations lead to level set models

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods

Explicit representations

◮ Explicit representations involve a parameterization (embedding) from an abstract parameter space Λ ⊂ Rn−1 to Rn n = 2 C =

• x(u) =

x(u), y(u) T u ∈ Λ

• n = 3

C =

• x(u) =

x(u), y(u), z(u) T u = (u, v) ∈ Λ

• ◮ The parameterization is obviously not unique

the derivatives Dp

ux convey both metric (extrinsic) and shape

(intrinsic) information choosing an arclength (normalized metrics) results into an intrinsic representation

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods

Implicit representations

◮ Implicit representations define the manifold as the zero level set

• f an implicit function over Ω ⊂ Rn

C =

• x
• ϕ(x) = 0
• ◮ The implicit function is obviously not unique

the (Cartesian) derivatives Dpϕ convey shape information only the representation is intrinsic

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Outline

1

Foundations

2

Parameterized active contours

3

Level set methods

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Internal energies

◮ Prior knowledge on C is specified via an internal energy Eint(C) which sets constraints on evolution and feasible segmentations energy additivity provides a flexible mechanism for integrating multiple priors Eint(C) =

• E i

int(C)

active contours are physically-based models implementing linear elasticity priors acting as shape constraints

– membrane prior → stretching / shrinking – thin-plate prior → bending – balloon prior → inflation / deflation

These physical priors are interpreted geometrically as smoothness constraints

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Membranes

◮ A membrane is a codimension-1 elastic material which can only deform by stretching / contracting along its tangent space Tx n = 2 n = 3 xu Tx = span(xu) xu xv Tx = span(xu, xv)

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

2D membrane

2D membrane energy E(C) =

• α(u) |xu|2 du

Hyperparameter: nonnegative function α over Λ ◮ Elastic properties → deformation α is a stiffness density defining local resistance to stretching

setting α(u0) = 0 enables free deformation at x(u0)

the standard model is the homogeneous membrane for which α(u) = α > 0 for all u ∈ Λ

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

2D membrane

2D membrane energy E(C) =

• α(u) |xu|2 du

◮ Elastic properties → deformation a membrane is the continuum mechanics analog of a spring

x1 x2 x0 xi xi-1 αi

Proof: spring network N E(N) =

• αi |xi − xi−1|2

Letting i → ∞ and xi−1 → xi xi − xi−1 = xi − xi−1 δu δu → xu du αi → α(u)

• Nicolas Rougon

IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

2D membrane

2D membrane energy E(C) =

• α(u) |xu|2 du

◮ Smoothness properties → regularization E is a 1st-order Tikhonov stabilizer

minimizing E over some shape space yields a.e. C1-continuous solutions

α is a regularization function controlling C1-continuity locally

smoothness at x(u) increases with α(u) setting α(u0) = 0 allows C to be C1-discontinuous at x(u0)

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

2D membrane

2D membrane energy E(C) =

• α(u) |xu|2 du

◮ Global geometric properties → regularization For a homogeneous membrane (α(u) = α) parameterized by an arclength s (|xs| = 1) E(C) = α

• ds = α Length (C)

A membrane is a minimal length segmenting curve (geodesic)

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

3D membrane

3D membrane energy E(C) = w10(u) |xu|2 + w01(u) |xv|2 du Hyperparameters: nonnegative functions (wij)i+j=1 over Λ ◮ Elastic properties → deformation (wij) are directional stiffness densities defining local resistance to stretching along coordinate axes

setting w10(u0) = 0 (w01(u0) = 0 ) enables free deformation along u (v) at x(u0)

the standard model is the homogeneous isotropic membrane for which wij(u) = α > 0 for all u ∈ Λ

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

3D membrane

3D membrane energy E(C) = w10(u) |xu|2 + w01(u) |xv|2 du ◮ Smoothness properties → regularization E is a 1st-order Tikhonov stabilizer

minimizing E over some shape space yields a.e. C1-continuous solutions

(wij) are regularization functions controlling C1-continuity locally

smoothness at x(u) increases with wij(u) setting wij(u0) = 0 allows C to be C1-discontinuous at x(u0)

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Thin-plates

◮ A thin-plate (also called shell) is a codimension-1 elastic material which can only deform by bending along its normal n n = 2 n = 3 n n xss n

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

2D thin-plate

2D thin-plate energy E(C) =

• β(u) |xuu|2 du

Hyperparameter: nonnegative function β over Λ ◮ Elastic properties → deformation β is a bending stiffness density defining local resistance to flexion

setting β(u0) = 0 enables free deformation at x(u0)

the standard model is the homogeneous thin-plate for which β(u) = β > 0 for all u ∈ Λ

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

2D thin-plate

2D thin-plate energy E(C) =

• β(u) |xuu|2 du

◮ Smoothness properties → regularization E is a 2nd-order Tikhonov stabilizer

minimizing E over some shape space yields a.e. C2-continuous solutions

β is a regularization function controlling C2-continuity locally

smoothness at x(u) increases with β(u) setting β(u0) = 0 allows C to be C2-discontinuous at x(u0)

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

3D thin-plate

3D thin-plate energy E(C) = w20(u) |xuu|2 + 2w11(u) |xuv|2 + w02(u) |xvv|2 du Hyperparameters: nonnegative functions (wij)i+j=2 over Λ ◮ Elastic properties → deformation (wij) are directional bending stiffness densities defining local resistance to flexion

setting wij(u0) enables free deformation at x(u0)

the standard model is the homogeneous isotropic thin-plate for which wij(u) = β > 0 for all u ∈ Λ

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

3D thin-plate

3D thin-plate energy E(C) = w20(u) |xuu|2 + 2w11(u) |xuv|2 + w02(u) |xvv|2 du ◮ Smoothness properties → regularization E is a 2nd-order Tikhonov stabilizer

minimizing E over some shape space yields a.e. C2-continuous solutions

(wij) are regularization functions controlling C2-continuity locally

smoothness at x(u) increases with wij(u) setting wij(u0) = 0 allows C to be C2-discontinuous at x(u0)

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Balloons

◮ A balloon is a codimension-1 elastic material which can only deform by inflating / deflating along its normal n n = 2 n = 3 n n ⊥ xu n n xu ∧ xv

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

2D balloon

2D pressure energy E(C) = ±

• det

p(x), xu du

Hyperparameter: irrotational vector field p over Ω ⊂ R2 ◮ Elastic properties → deformation p is a pressure density defining local resistance to inflation

setting p(x(u0)) = 0 enables free deformation at x(u0)

the standard model is the homogeneous isotropic pressure for which p(x) = p x (p > 0) for all x ∈ Ω

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

2D balloon

2D pressure energy E(C) = ±

• det

p(x), xu du

◮ Global geometric properties → regularization For a homogeneous isotropic pressure (p(x) = p x) E(C) = ±p

• det

x, xu du = ±p · Area (Rin)

A 2D balloon is a minimal / maximal area segmenting curve

Proof: use Green-Stokes theorem to transform the boundary integral into a domain integral

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

2D balloon

Theorem (Green-Stokes) Given a C1 vector field V : Ω → R2 and a compact domain R ⊂ Ω with smooth boundary ∂R (unit normal n)

• R

∇ · V(x) dx =

• ∂R

(V · n) (x(s)) ds Letting V =

P, Q T, this rewrites as

• R

(Px + Qy) dxdy =

• ∂R

(Pyu − Qxu) du =

• ∂R

det (V, xu) du Setting V(x) = x, we get 2 Area (R) =

• det (x, xu) du
• Nicolas Rougon

IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

3D balloon

3D pressure energy E(C) = ±

• det

p(x), xu, xv du

Hyperparameter: irrotational vector field p over Ω ⊂ R3 ◮ Elastic properties → deformation p is a pressure density defining local resistance to inflation

setting p(x(u0)) = 0 enables free deformation at x(u0)

the standard model is the homogeneous isotropic pressure for which p(x) = p x (p > 0) for all x ∈ Ω

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

3D balloon

3D pressure energy E(C) = ±

• det

p(x), xu, xv du

◮ Global geometric properties → regularization For a homogeneous isotropic pressure (p(x) = p x) E(C) = ±p

• det

x, xu, xv du = ±p · Volume (Rin)

A 3D balloon is a minimal/maximal volume segmenting surface

Proof: use Green-Ostrogradski theorem to transform the boundary integral into a domain integral

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

3D balloon

Theorem (Green-Ostrogradski) Given a C1 vector field V : Ω → R3 and a compact domain R ⊂ Ω with smooth boundary ∂R (unit normal n =

xu∧xv |xu∧xv|)

• R

∇ · V(x) dx =

• ∂R

V · dS =

• ∂R

(V · n) (x(u)) |xu ∧ xv| du Letting V =

P, Q, R T and using a · (b ∧ c) = det (a, b, c)

this rewrites as

• R

(Px + Qy + Rz) dx =

• ∂R

det (V, xu, xv) du Setting V(x) = x, we get 3 Volume (R) =

• det (x, xu, xv) du
• Nicolas Rougon

IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Internal energies

Eint(C) =

• Vint

x(u) du

where: Vint = Vmembrane + Vthin−plate + Vballoon 2D active contour Vint(x) = α(u) |xu|2 + β(u) |xuu|2 ± det

p(x), xu)

3D active contour Vint(x) = w10(u) |xu|2 + w01(u) |xv|2 + w20(u) |xuu|2 + 2w11(u) |xuv|2 + w02(u) |xvv|2 ± det

p(x), xu, xv)

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

External energies

◮ The interactions of C with image content are specified via an external energy Eext(C, L) energy additivity provides a flexible mechanism for integrating multiple relevant image information for the problem at hand Eext(C, L) =

• λi E i

ext(C, L)

the influence of image energy E i

ext in the overall segmentation

criterion is controlled by the hyperparameter λi

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

External energies

◮ Image energies E i

ext are derived from potential maps V i ext(L)

computed from image data over Ω contour-based E i

ext(C, L) =

• Λ

V i

ext(L)(x(u)) du

region-based E i

ext(C, L) =

• Rin

V i

ext(L)(x) dx

◮ Active contour modeling relies on designing relevant image potentials for the application at hand. V i

ext(L) are expected to

depend on discriminative image features be minimal along C (if contour-based) or over Rin (if region-based)

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

External energies

◮ Luminance potentials deterministic target V i

ext(L)

properties Dark/bright object dark bright L −L contour / region-based Low-texture object region variance (L − µin)2 µin : mean luminance over Rin ≡ constant image model (a.k.a. cartoon model)

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

External energies

◮ Luminance potentials statistical target V i

ext(L)(x)

properties High-texture

• bject

region log-likelihood − log pRinL(x)

• pRin: pdf of luminance over Rin

parametric estimator pRin(L | Θ)

e.g. Gaussian Mixture Model

Gaussian pdf ≡ cartoon model non-parametric estimator

e.g. kernel density estimator

supervised (pRin learnt from a training set) | unsupervised setting

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

External energies

◮ Edge potentials target V i

ext(L)

properties High-contrast boundary polynomial bounded decreasing −|∇L|p ϕ(|∇L|) p ≥ 1 ϕ(0) = 1 , lim

x→+∞ ϕ(x) = 0

Laplacian zero-crossing |∆L|p p ≥ 1 Binary edge map C ◮ medium | large capture range distance-based convolutive dC −Kσ ⋆ C distance function to C Kσ : kernel with width σ e.g. Gaussian

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

External energies

◮ Example application Hand tracking using active contours

1 Feature extraction

Canny-Deriche edge detection edge map post-processing

raw edge map hand silhouette C 2 Potential map generation

chamfer distance dC to C

Vext(L) Vext(L) level lines Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Total energy

◮ Summing elastic and image energies yields a total energy E(C) ≡ segmentation criterion E(C) = Eint(C) + Eext(C, L) =

• V

x(u) du

where: V = Vint + λi V i

ext(L)

◮ Admissible segmentations C∗ are minimizers of E(C) over the shape space S C∗ = arg min

C ∈ S E(C)

C∗ is an equilibrium state between internal / external actions, derived from elastic / image energies ◮ Modeling issues shape space S

• ptimization technique

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Shape spaces

◮ Finite-dimensional shape spaces → parametric estimation The coordinates of x ∈ C can be expressed in some basis (φi)i x(u) =

• θi φi(u)

C is completely defined by Θ = (θi)i ≡ model parameters E(Θ) is computed in closed-form and minimized w.r.t. Θ standard bases: B-spline B-snakes | eiu Fourier-snakes | wavelets wavelet-snakes | ... ◮ Infinite-dimensional shape spaces → nonparametric estimation the location of every point x ∈ C must be estimated E(C) is minimized w.r.t. x

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Energy minimization

◮ Minimization techniques iteratively decrease E(C) by deforming C until a minimum is reached Direct (a.k.a. gradient-free) methods operate on E(C) Variational methods use the derivative of E(C)

dim S < ∞ ∂ΘE(Θ) =

• ∂θiE(Θ)
• i is a standard (vector) derivative

dim S = ∞ ∂xE(C) is a functional derivative ≡ force F(x) acting on x ∈ C

F(x) = ∂xEint(C) + ∂xEext(C, L) = Fint(x) + Fext(x)

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Gradient-free minimization

◮ Gradient-free optimizers are search methods applicable to discrete variables and feasible spaces C is sampled → point set / mesh a discrete expression of E(C) is derived finite difference | finite element methods candidate deformations of C are hypothetized from a finite set of possible variations for x or Θ for each of them, energy variation is computed the deformation maximizing energy decrease is applied this process is iterared until E(C) cannot be ց

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Gradient-free minimization

◮ Efficient algorithms are available dynamic programming Nelder-Mead (a.k.a. nonlinear) simplex . . . ◮ Performances fast robust to noise global | “good” local minimizer accuracy depends on sampling density model point clustering around attractors of Vext → fixed by periodic resampling / remeshing during iterations

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Variational calculus

◮ Let E(f ) =

V (f ) be a functional over some Hilbert function

space with standard dot product < f , g > =

fg

Gâteaux derivative When it exists, the limit DηE(f ) = lim

ε → 0

E(f + εη) − E(f ) ε defines the Gâteaux derivative of E at f in the direction η The Gâteaux derivative verifies DηE(f ) = < ∂f V , η > ∂f V is referred to as the variational derivative of E(f )

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Variational calculus

◮ The extrema of E(f ) verify DηE(f ) = < ∂f V , η > = 0 ∀η Euler-Lagrange equations The minimizers of E(f ) are solutions of the PDE ∂f V = 0 known as Euler-Lagrange equations ◮ Hereafter, we derive general closed-form expressions for ∂f V when f : Ω ⊂ Rp → Rq Key results for solving a variety of image understanding issues modeled as deterministic optimization problems

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Variational calculus

◮ We first consider the case of a functional E(f ) =

V (f , fu) du

• ver a space of univariate scalar functions f (u)

Performing Taylor expansion E(f + εη) =

• V (f + εη, fu + εηu) du

≈ V (f , fu) + ε ∂V ∂f η + ε ∂V ∂fu ηu

• du

= E(f ) + ε

∂V

∂f η + ∂V ∂fu ηu

• du

Integrating by parts (

b

a g′h = [gh]b a −

gh′)

E(f +εη) ≈ E(f ) + ε

∂V

∂f − ∂ ∂u

∂V

∂fu

• η du + ε

∂V

∂fu η

b

a

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Variational calculus

Enforcing zero boundary conditions (BCs) for η η(a) = η(b) = 0 E(f + εη) ≈ E(f ) + ε

∂V

∂f − ∂ ∂u

∂V

∂fu

• η du

Taking the limit lim

ε → 0

E(f + εη) − E(f ) ε =

∂V

∂f − ∂ ∂u

∂V

∂fu

• , η
• Hence the variational derivative ∂f V

∂f V = ∂V ∂f − ∂ ∂u

∂V

∂fu

• Nicolas Rougon

IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Variational calculus

◮ Consider an arbitrary functional E(f ) =

V (. . . , fu(i), . . . ) du

involving the ith derivative fu(i) of f Performing Taylor expansion E(f + εη) =

• V (. . . , fu(i) + εηu(i), . . . ) du

≈ E(f ) + ε . . . + ∂V ∂fu(i) ηu(i) + . . .

• du

Integrating by parts i times with zero BCs for all ηu(k) (k ≤ i) E(f +εη) ≈ E(f ) + ε . . . (−1)i ∂i ∂ui

• ∂V

∂fu(i)

• + . . .
• η du

Taking the limit lim

ε → 0

E(f + εη) − E(f ) ε =

• . . . (−1)i ∂i

∂ui

• ∂V

∂fu(i)

• + . . . , η
• Nicolas Rougon

IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Variational calculus

This proves that ∂f V has a generating term (−1)i ∂i ∂ui

• ∂V

∂fu(i)

• riginating from the occurrence of the derivative fu(i) in V

(using the convention fu(0) = f ) Variational derivative (univariate scalar problems) Given a functional E(f ) =

V (f , fu, . . . , fu(n)) du with f (u) ∈ R

∂f V =

n

• i=0

(−1)i ∂i ∂ui

• ∂V

∂fu(i)

• Nicolas Rougon

IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Variational calculus

◮ Example: Let E(f ) = ff 2

u + f 2 uu

• du

E(f ) depends on a potential of the form V (f , fu, fuu). Hence ∂f V = ∂V ∂f − ∂ ∂u

∂V

∂fu

• + ∂2

∂u2

∂V

∂fuu

• We get

∂f V = f 2

u − ∂

∂u (2ffu) + ∂2 ∂u2 (2fuu) = f 2

u − 2f 2 u − 2ffuu + 2fuuuu

= −f 2

u − 2ffuu + 2fuuuu

• Nicolas Rougon

IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Variational calculus

◮ We now consider a functional E(f ) =

V (f , fu) du over a

space of univariate vector functions f (u) =

f j(u) ∈ Rq

V (f , fu) can be viewed as a function V (f 1, . . . , f q, f 1

u , . . . , f q u )

• f 2q independent scalar variables

Performing similar computations coordinatewise readily yields ∂f jV = ∂V ∂f j − ∂ ∂u

• ∂V

∂f j

u

• Letting

∂ ∂f =

• ∂

∂f j

• and

∂ ∂fu =

• ∂

∂f j

u

• leads to the vector form

∂f V = ∂V ∂f − ∂ ∂u

∂V

∂fu

• ∈ Rq

→ Vector generalization of the scalar case

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Variational calculus

◮ These conclusions extend to the case of an arbitrary functional E(f )=

V (. . . fu(i) . . . ) du involving higher-order derivatives fu(i)

Variational derivative (univariate vector problems) Given a functional E(f ) =

V (f , fu, . . . , fu(n)) du with f (u) ∈ Rq

∂f V =

n

• i=0

(−1)i ∂i ∂ui

• ∂V

∂fu(i)

• The generating term of ∂f V originates from the occurrence of

the derivative fu(i) in V ◮ Use cases 2D parameterized active contours | minimal paths

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Variational calculus

◮ We next consider a functional E(f ) =

V (f , fu1, fu2) du1du2

• ver a space of bivariate scalar functions f (u1, u2)

Performing Taylor expansion E(f + εη) =

• V (f + εη, fu1 + εηu1, fu2 + εηu2) du1du2

≈ E(f ) + ε

∂V

∂f η + ∂V ∂fu1 ηu1 + ∂V ∂fu2 ηu2

• du1du2

Using Green’s theorem (

• Ω (Pu1 + Qu2) =
• ∂Ω (Pys − Qxs) )

E(f + εη) ≈ E(f ) + ε

∂V

∂f − ∂ ∂u1

∂V

∂fu1

• −

∂ ∂u2

∂V

∂fu2

• η du1du2

+ ε

• ∂Ω

η

∂V

∂fu1 ys − ∂V ∂fu2 xs

• ds

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Variational calculus

Enforcing zero BCs for η along ∂Ω E(f +εη) ≈ E(f ) + ε

∂V

∂f − ∂ ∂u1

∂V

∂fu1

• − ∂

∂u2

∂V

∂fu2

• η du1dvu2

Hence the variational derivative ∂f V ∂f V = ∂V ∂f − ∂ ∂u1

∂V

∂fu1

• −

∂ ∂u2

∂V

∂fu2

• → Bivariate extension of the univariate case

◮ Extension to functionals E(f ) =

V (f , fu1, . . . , fup) du1 . . . dup

• ver a space of multivariate scalar functions f (u1, . . . , up)

∂f V = ∂V ∂f −

p

• k=1

∂ ∂uk

∂V

∂fuk

• Nicolas Rougon

IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Variational calculus

◮ Given a multi-index i = (i1, . . . , ip) ∈ Np and u = (u1, . . . , up) denote |i| = i1 + · · · + ip denote

∂i ∂ui = ∂ i1+···+ip ∂ui1

1 ...∂u ip p

and fu(i) = fu(i1)

1

...u

(ip) p

◮ Consider an arbitrary functional E(f ) =

V (. . . fu(i) . . . ) du

involving the partial derivative fu(i) of total order |i| of f Performing Taylor expansion and applying Green’s theorem |i| times using zero BCs for all ηu(k) (|k| ≤ |i|) E(f +εη) ≈ E(f ) + ε . . . (−1)|i| ∂i ∂ui

• ∂V

∂fu(i)

• + . . .
• η du

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Variational calculus

Variational derivative (multivariate scalar problems) Given a functional E(f ) =

V (f , fu . . . fu(n)) du with f (u) ∈ R

∂f V =

|n|

• |i|=0

(−1)|i| ∂i ∂ui

• ∂V

∂fu(i)

• ◮ Use cases

level set active contours gray level image restoration (denoising, deblurring...) shape from X

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Example 1 | Variational 2D image denoising

E(L) =

• Ω

L(x) − Ln(x) 2dx + λ

• Ω

L2

x + L2 y

dx

Gaussian noise model + 1st-order quadratic stabilizer ∇L2 Since E(L) =

• V (L, Lx, Ly) dxdy

∂LV = ∂V ∂L − ∂ ∂x

∂V

∂Lx

• − ∂

∂y

• ∂V

∂Ly

• We get

∂f V = 2(L − Ln) − λ

∂

∂x (2Lx) + ∂ ∂y (2Ly)

• = 2

L − Ln − λ∆L

• Nicolas Rougon

IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Example 2 | Variational 2D image denoising

E(L) =

• Ω

L(x) − Ln(x) 2dx + λ

• Ω

L2

xx + 2L2 xy + L2 yy

dx

2nd-order quadratic stabilizer D2L2 = tr

• (D2L)2

Since E(L) =

• V (L, Lxx, Lxy, Lyy) dxdy

∂LV = ∂V ∂L + ∂2 ∂x2

∂V

∂Lxx

• +

∂2 ∂x∂y

• ∂V

∂Lxy

• + ∂2

∂y2

• ∂V

∂Lyy

• We get

∂f V = 2(L − Ln) + λ

• ∂2

∂x2 (2Lxx) + ∂2 ∂x∂y (4Lxy) + ∂2 ∂y2 (2Lyy)

• = 2
• L − Ln + λ (Lxxxx + 2Lxxyy + 2Lyyyy)
• = 2
• L − Ln + λ ∆∆L
• ∆∆ is the biharmonic operator

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Variational calculus

◮ The previous result extends to multivariate vector functions f (u) =

f j(u) ∈ Rq in a coordinatewise fashion

Variational derivative (multivariate vector problems) Given a functional E(f ) =

V (f , fu . . . fu(n)) du with f (u) ∈ Rq

∂f V =

|n|

• |i|=0

(−1)|i| ∂i ∂ui

• ∂V

∂fu(i)

• ◮ Use cases

parameterized active surfaces | surface reconstruction

• ptical flow estimation | image registration

color image / vector field restoration

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Elastic forces

Eint(C) =

• Vint

x(u) du

where: Vint = Vmembrane + Vthin−plate + Vballoon

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

2D internal forces

type Vint Fint(x) membrane α |xu|2 − ∂ ∂u

α xu

• thin-plate

β |xuu|2 ∂2 ∂u2

β xuu

• balloon

det

p(x), xu

• −

∇ · p |xu| n

quadratic potentials yield quasi-linear elastic forces pressure forces induce coordinate coupling in Euler equations

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

2D balloon forces

◮ Proof: Let p = (p1, p2)T Vballoon =

• p1(x)

xu p2(x) yu

• = p1(x) yu − p2(x) xu = V (x, xu)

Computing ∂xV = ∂V ∂x − ∂ ∂u

∂V

∂xu

• yields

Fballoon(x) =

   

∂p1 ∂x yu − ∂p2 ∂x xu + ∂p2 ∂u ∂p1 ∂y yu − ∂p2 ∂y xu − ∂p1 ∂u

   

Using chain rule

∂

∂u = xu ∂ ∂x + yu ∂ ∂y = xu · ∇

• Fballoon(x) =

∂p1

∂x + ∂p2 ∂y yu −xu

• = −

∇·p |xu| n

• Nicolas Rougon

IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

3D internal forces

type Vint Fint(x) membrane w10 |xu|2 + w01 |xv|2 − ∂ ∂u

w10 xu − ∂

∂v

w01 xv

• thin-plate

w20 |xuu|2 + w02 |xvv|2 + 2w11 |xuv|2 ∂2 ∂u2

w20 xuu + ∂2

∂v2

w02 xvv

• + 2

∂2 ∂u∂v

w11 xuv

• balloon

det

p(x), xu, xv

• −

∇ · p |xu ∧ xv| n

quadratic potentials yield quasi-linear elastic forces pressure forces induce coordinate coupling in Euler equations

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Elastic forces

◮ Owing to their linearity, elastic forces are compactly expressed via a linear differential operator K, known as the stiffness tensor Felastic(x) = Fmembrane(x) + Fthin−plate(x) = Kx K n = 2 ∂2 ∂u2

• β ∂2

∂u2

• − ∂

∂u

• α ∂

∂u

• n = 3

∂2 ∂u2

• w20

∂2 ∂u2

• + 2

∂2 ∂u∂v

• w11

∂2 ∂u∂v

• + ∂2

∂v2

• w02

∂2 ∂v2

• − ∂

∂u

• w10

∂ ∂u

• − ∂

∂v

• w01

∂ ∂v

• Nicolas Rougon

IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Image forces

type Eext Fext(x) contour-based

• Λ

Vext(L)(x(u)) du ∇Vext(L)(x) region-based

• Rin

Vin(L)(x) dx −Vin(L)(x) n image forces are nonlinear

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Region-based image forces

◮ Proof: Assume Vin(L) is smooth and does not depend on Rin There exists a smooth vector field V : Ω → Rn such that Vin(L) = ∇ · V Using Green’s theorem (|xs| = |xs1 ∧ xs2| = 1)

n = 2

• Rin

Vin(L)(x) dx =

• ∂Rin

det (V(x), xs) ds n = 3

• Rin

Vin(L)(x) dx =

• ∂Rin

det (V(x), xs1, xs2) ds

→ Same expressions as the balloon energy Using the expression of the balloon force Fin(x) = −

∇ · V n = −Vin(L)(x) n

• Nicolas Rougon

IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

2-phase active region segmentation

◮ A generic 2-phase active region segmentation model is obtained by integrating data consistency constraints over the object region Rin and the background Rout = Ω\Rin Eext(C, L) = λin

• Rin

Vin(L)(x) dx + λout

• Rout

Vout(L)(x) dx Since object (background) boundaries C = ∂Rin (∂Rout) have opposite normals, the resulting image force is Fext(x) =

• λoutVout(L)(x) − λinVin(L)(x)
• n

◮ The segmenting manifold evolves according to a region competition principle

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Variational optimization

◮ Since image forces are nonlinear, Euler equations cannot be solved in closed-form ◮ This issue is overcome by using numerical optimization methods Starting from an initialization, these techniques iteratively update the variable x to be estimated (i.e. deform C) so as to decrease E(C) until a minimum is reached A time variable t is introduced C(t) =

• x(u, t)
• (u, t) ∈ Λ × R+

The evolution law for x is driven by ∂xV , and can be mathematically- or physically-based Relevant schemes must ensure fast and stable convergence, and avoid non significant local minima

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Gradient-descent methods

◮ Gradient-descent methods are 1st-order schemes in which the variable x to be optimized is updated along the steepest slope −∂xV of the energy surface

C E deformatjon initjalizatjon segmentatjon

• δxV

They involve a speed law with prototype xt = −∂xV (x) Many schemes are available

– fixed | adaptive time step – conjugate gradient – Levenberg-Marquardt – Newton | quasi-Newton | BFGS – stochastic gradient – . . .

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Gradient-descent methods

◮ Bibliography A practical overview of gradient-free & gradient-descent methods with C++ implementations can be found in

• W. Press - S. Teukolsky - W. Vetterling - B. Flannery

Numerical Recipes - The Art of Scientific Computing Cambridge University Press, 2007 (3rd Edition) ◮ For the sake of simplicty, hereafter we focus our attention on the basic fixed time step gradient-descent

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Physically-based methods

◮ Physically-based methods rely on Lagrangian mechanics, leading to 2nd-order schemes derived from the least action principle

∂ ∂t (M xt)

+ Cxt + F(x) =

inertia dissipation applied forces

Hyperparameters: mass tensor M, damping tensor C the standard model involve homogeneous isotropic mass and damping distributions (M = m I, C = γ I) (m, γ > 0) m xtt + γ xt + F(x) = 0

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Discretization

◮ Discretizing the model requires Sampling C(t) (or, equivalently, the parameter space Λ)

– this yields a point set xt

i = x(ui, t) ∈ Ω ⊆ Rn ( i ∈ [1..N])

– C(t) is then represented as a vector of points Xt = (xt

i )

Deriving discrete expressions for the forces F(xt

i )

Discretizing the evolution law for C(t) at samples xt

i

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Discretizing forces

◮ Image forces An image force map Fext over the pixel grid is readily derived from the potential map Vext

– for contour-based image energies defined by smooth potentials, Fext = ∇Vext is estimated using a simple Sobel filter

Image forces Fext(xi) at contour points (recall xi ∈ Rn) are then computed by interpolation from the map Fext

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Discretizing forces

◮ Internal forces involve spatial derivatives along C The latter can be discretized using finite differences (FD) xu(ui) ≈ D−

u xi = xi − xi−1

δu xuu(ui) ≈ D2

u xi = xi+1 − 2xi + xi−1

δu2 xuuuu(ui) ≈ D2

u

D2

u xi

= xi+2 − 4xi+1 + 6xi − 4xi−1 + xi−2

δu4 where δu is the spatial sampling step The stiffness tensor K is discretized into a stiffness matrix K acting on the discrete shape vector X

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Stiffness matrix

◮ Example: 2D closed contour with homogeneous elastic properties K =

                 

β β β −α−4β β β β −α−4β −α−4β 2α+6β β −α−4β −α−4β β 2α+6β −α−4β

                 

◮ Properties of K arbitrary n / material functions positive symmetric band, hence sparse band size depending

• n FD scheme

cyclic if C closed singular letting 1 = (1)i

• ne has K1 = 0

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Discrete evolution law

◮ The evolution law involves time derivatives The latter are discretized using FD xt(ui, t) ≈ D−

t xt i = xt i − xt−1 i

δt xtt(ui, t) ≈ D−

t

D−

t xt i

= xt

i − 2xt−1 i

+ xt−2

i

δt2 where δt is the time step

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Discrete speed law

◮ For simplicity, we consider a (1st-order) speed law xt = −F(x) Dividing applied forces into linear (Felastic) and nonlinear (Fnl = Fballoon + Fext) parts F(x) = Kx + Fnl(x) the evolution law splits into a linear l.h.s. and a nonlinear r.h.s. xt + Kx = −Fnl(x) Several discretization options are available, depending on the time at which force terms are evaluated xt

i − xt−1 i

δt + Kx(?)

i

= −Fnl(x(?)

i

) leading to various numerical iterative schemes

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Discrete speed law

◮ Explicit scheme All force terms are evaluated at time t − 1 xt

i − xt−1 i

δt + Kxt−1

i

= −Fnl(xt−1

i

) This yields a computationally simple scheme where current estimates are computed directly from previous ones xt

i =

1 − δt K xt−1

i

− δt Fnl(xt−1

i

)

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Discrete speed law

◮ Implicit scheme All force terms are evaluated at time t xt

i − xt−1 i

δt + Kxt

i = −Fnl(xt i )

Computing xt

i requires inverting the update equation l.h.s

1 + δt K xt

i + δt Fnl(xt i ) = xt−1 i

Intractable due to the nonlinearity of image forces

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Discrete speed law

◮ Semi-implicit scheme Linear (resp. nonlinear) terms are evaluated at time t (resp. t −1) xt

i − xt−1 i

δt + Kxt

i = −Fnl(xt−1 i

) This yields a mildly demanding scheme xt

i =

1 + δt K −1

xt−1

i

− δt Fnl(xt−1

i

)

• where the time-independent (sparse, symmetric) matrix

1 + δt K can be pre-inverted efficiently

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Stable convergence

◮ The time step δt conditions the convergence speed C E

• δt F(xt)

xt⁺¹ xt⁺² xt

To avoid divergence / oscillations around minima, the time step must be upper-bounded δt ≤ (δt)max This constraint is known as the Courant-Friedrichs-Lewy (CFL) condition for the numerical scheme The CFL condition is (only) necessary for convergence The value of (δt)max depends on the FD scheme

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Stable convergence

◮ Implicit solvers are usually less prone to numerical instability (δt)explicit

max

≤ (δt)semi−implicit

max

≪ (δt)implicit

max

Explicit schemes have low computational cost per iteration but require more iterations to converge Semi-implicit schemes have higher computational cost per iteration but converge more rapidly The choice of a numerical scheme should be dictated by the (average) distance of initialization to solution

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Energies Energy minimization Model discretization

Discrete Lagrangian evolution law

◮ All previous observations hold for the 2nd-order evolution law of Lagragian mechanics m xtt + γ xt + F(x) = 0 which lead to more complex but more stable discrete schemes m xt

i − 2xt−1 i

+ xt−2

i

δt2 + γ xt

i − xt−1 i

δt + Kx(?)

i

= −Fnl(x(?)

i

)

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Outline

1

Foundations

2

Parameterized active contours

3

Level set methods

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Limitations of parameterized active contours

◮ Explicit representation ◮ implicit representation iterative reparameterization (≡ resampling / remeshing) → computational complexity fixed topology → topological rigidity ◮ Energy properties ◮ energy redefinition non-intrinsic → solutions depend on metric non-convex, hence non-significant local minima → sensitivity to initialization global smoothness → non-preservation of shape discontinuities

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Geodesic active contours

◮ Homogeneous membrane in a contrast potential g2|∇L|

• g bounded, decreasing s.t. g(x) +∞

− − → 0 E0(C, L) =

• α |xu|2du + λ
• g2|∇L|(x)

du |xu| g(|∇L|(x)) x xu

E0 can be viewed as the square on the hypotenuse of the right triangle T with edges {√α |xu|, √ λ g(x)} in phase space Span(x, xu) Maupertuis principle E0 has the same minima as the intrinsic energy E ≡ Area (T ) E(C, L) =

• g

|∇L|(x) |xu| du =

• g

|∇L|(x) ds

This model is known as geodesic active contours

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Geodesic active contours

E(C, L) =

• g

|∇L|(x) ds

◮ 2 interpretations by analogy with a homogeneous membrane Emembrane(C) =

• α ds

Energy of a membrane with image-dependent stiffness g(|∇L|)

– since g(u)

+∞

− − → 0, stiffness decreases along image edges. This allows to account for corner points

Length of an image curve using the image-dependent Riemannian metric dσ = g

|∇L| ds – segmentations are minimal length curves (geodesics) for dσ

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Geodesic active contours

◮ Variational derivative ∂xV = η(x) n where the force magnitude η(x) is η(x) = ∇g · n − gk ∂xV involves shape properties only (normal n, curvature k) (i.e. is metric-free). This is because E is intrinsic A more efficient model is obtained by adding a pressure term η(x) = ∇g · n − g (k ± p) ◮ Geometric speed law xt = −η(x) n

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Geodesic active contours

◮ Proof: Let V (x, xu) = g

|∇L|(x) |xu| and choose u ≡ s

Computing ∂xV = ∂V ∂x − ∂ ∂s

∂V

∂xs

• yields

∂xV = ∇g |xs| − ∂ ∂s

gt

• = ∇g − ∂g

∂s t − gk n Using chain rule

∂

∂s = xs ∂ ∂x + ys ∂ ∂y = t · ∇

and

expressing ∇g in frame (t, n)

∇g = (∇g · t) t + (∇g · n) n

• ∂xV =

gk − ∇g · n n

• Nicolas Rougon

IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Geodesic active surfaces

◮ Generalization to a surface C ⊆ R3 in a 3D contrast potential g

|∇L| is achieved by substituting the arclength ds = |xu| du

with the unit surface element da = |xu ∧ xv| du E(C, L) =

• g

|∇L|(x) |xu ∧ xv| du

◮ Variational derivative ∂xV = η(x) n η(x) = ∇g · n − gH E is intrinsic, so that ∂xV involves shape properties only (normal n, mean curvature H) A more efficient model is obtained by adding a pressure term η(x) = ∇g · n − g (H ± p)

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Level set representation

◮ The evolving manifold C(t) is represented as the zero level set of a smooth function ϕ over the image grid Ω C(t) =

• x ∈ Ω
• ϕ(x, t) = 0
• ◮ Shape / motion properties of C(t) relate to the derivatives of ϕ

unit normal n = − ∇ϕ |∇ϕ| (mean) curvature curv(ϕ) = ∇ · n speed Dϕ = ϕt + ∇ϕ · xt = 0 xt = −η(x) n ϕt = −η |∇ϕ|

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Level set representation

◮ Though ϕ is a smooth function, its zero level set can be non-smooth → image discontinuity preservation change its topolopy in time → topological adaptivity

C⁽t ⁾ φ = 0 φ (x,t) φ (x,t) φ = 0 C⁽t ⁾

fusion | split

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Geodesic active contours

◮ Level set-based evolution law ≡ reaction-diffusion PDE ϕt = ±p g |∇ϕ| + |∇ϕ| ∇ ·

• g

∇ϕ |∇ϕ|

• hyperbolic

parabolic

ϕ is initialized as the normalized signed distance to C(0) ϕ(x, 0) = 1 max

Ω dC(0)(x)

• −dC(0)(x)

if x ∈ Rin +dC(0)(x) if x ∈ Rout

−dC⁽⁰⁾(x) −dC⁽⁰⁾(x) dC⁽⁰⁾(x) C⁽⁰⁾ C⁽⁰⁾

From the binary image of C(0) dC(0) is computed using a discrete distance function sign (≡ orientation) is assigned from a connected component labeling

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Chamfer distances

◮ Estimation of geodesic distance map dij to object O ⊆ Ω ⊂ Z2 via sequential propagation of local integer approximations mij

• f Euclidean distance in pixel neighborhood N

d₁ d₁ d₂ d₂ d₁ d₁ d₂ d₂ m⁺ m-

N ⁺ N -

Initialization dij =

• if xij ∈ O

+∞

• therwise

Forward / backward scans using half-masks m± dij = min

(ı′′) ∈ N ±

• dij , di+ı′,j+′ + m±

ı′′

• Optimal coefficients mij

d1 = 3 d2 = 4 Normalization by

1 d1

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Chamfer distances

◮ Estimation of geodesic distance map dijk to object O ⊆ Ω ⊂ Z3 via sequential propagation of local integer approximations mijk

• f Euclidean distance in voxel neighborhood N

d₁ d₁ d₂ d₂ d₁ d₁ d₂ d₂ d₂ d₂ d₃ d₃ d₂ d₂ d₃ d₃ d₂ d₂ d₃ d₃ d₂ d₂ d₃ d₃ d₁ d₁

k’= 0 k’= -1 k’= 1

m⁺ m-

Initialization dijk =

• if xijk ∈ O

+∞

• therwise

Forward / backward scans using half-masks m± dijk = min

(ı′′k′) ∈ N ±

• dijk , di+ı′,j+′,k+k′ + m±

ı′′k′

• Optimal coefficients mijk

d1 = 3 d2 = 4 d3 = 5 Normalization by

1 d1

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Discrete evolution law

ϕt = ±p g |∇ϕ| + |∇ϕ| ∇ ·

• g ∇ϕ

|∇ϕ|

• ◮ The potential map g

|∇L| is computed analytically from the

(smooth) contrast map |∇L| ← Canny-Deriche gradient filter g

|∇L|

• ij = g

|∇L|ij

• ◮ The evolution law is discretized using FD

distinct schemes are used for the hyperbolic / parabolic terms explicit iterative scheme ϕ(t) = ϕ(t−1) + δt

• ±p g |∇ϕ| + |∇ϕ| ∇ ·
• g ∇ϕ

|∇ϕ|

(t−1)

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Discrete evolution law

◮ Parabolic term ≡ smoothing effect → standard FD schemes apply |∇ϕ| ∇ ·

• g ∇ϕ

|∇ϕ|

• = g |∇ϕ| curv(ϕ) + ∇g · ∇ϕ

∇g is computed using a Sobel gradient filter (g is smooth) both terms are computed analytically by plugging standard FD estimates of 1st / 2nd-order derivatives of ϕ ϕx(xij) ≈ Dxϕij = ϕi+1,j − ϕi−1,j 2δx ϕxx(xij) ≈ D2

xϕij = ϕi+1,j − 2ϕij + ϕi−1,j

δ2

x

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Discrete evolution law

◮ Hyperbolic term ≡ transport effect (speed η(x) = ±p g)

x y

η(x) C⁽t ⁾ C⁽t ⁺¹ ⁾

n- n⁺

Singularities arise along non-convex parts of level sets due to self-collisions ≡ shocks At shocks, ϕ is not differentiable. Yet, sided derivatives of ϕ exist and allow to define a proper normal motion direction n− or n+ Building C(t+1) geometrically by moving x ∈ C(t) according to the speed law does not yield a simply connected curve The proper way for building C(t+1) is Huygens’ construction → this is performed by monotonic conservative FD schemes

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Discrete evolution law

◮ Hyperbolic term ≡ transport effect Monotonic conservative FD schemes

x y

η(x) C⁽t ⁾ C⁽t ⁺¹ ⁾

n- n⁺

– lateral FD approximations of sided derivatives ϕ+

x (xij) ≈ D+ x ϕij = ϕi+1,j − ϕij

δx ϕ−

x (xij) ≈ D− x ϕij = ϕij − ϕi−1,j

δx – automated selection of relevant sided derivatives based on their sign using Boolean switch functions

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Discrete evolution law

◮ Hyperbolic term ≡ transport effect Monotonic conservative FD schemes

– forward motion (p > 0)

p g |∇ϕij| ≈ p g

ξ=x,y

min(D+

ξ ϕij, 0)2 + max(D− ξ ϕij, 0)2

– backward motion (p < 0)

−p g |∇ϕij| ≈ −p g

ξ=x,y

min(D−

ξ ϕij, 0)2 + max(D+ ξ ϕij, 0)2

C⁽t ⁾ C⁽t ⁺¹ ⁾ n- n⁺ n

straightforward 3D extension

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Narrow band method

◮ The computational load is reduced by updating ϕ only in a narrow band (NB) around its zero level set created by thresholding the initial distance function dC(0) Neuman boundary conditions (i.e. zero derivatives) are used ◮ A peripheral test band (TB) allows to warn for near NB exit created by binary eroding the NB ϕ is reinitialized to dC(t) whenever C(t) enters the TB

NB TB

C ⁽⁰⁾ C ⁽t ⁾

NB half-width ≈ 3-5 pixels TB width ≈ 1-2 pixels

initialization reinitialization

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Geodesic active contour segmentation

◮ Example applications Vessel segmentation in retinal angiography

initialization evolution segmentation

Vertebra segmentation in CT

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Geodesic active surface segmentation

◮ Example application: cortex segmentation in brain MRI

• riginal dataset

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Geodesic active surface segmentation

◮ Example application: cortex segmentation in brain MRI

initialization iteration #10 iteration #20 iteration #30 iteration #40 iteration #80 iteration #120 iteration #150 Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Reminder | 2-phase active region segmentation

◮ A generic 2-phase active region segmentation model is obtained by integrating data consistency constraints over the object region Rin and the background Rout = Ω\Rin Eext(C, L) = λin

• Rin

Vin(L)(x) dx + λout

• Rout

Vout(L)(x) dx The resulting image force is Fext(x) =

• λoutVout(L)(x) − λinVin(L)(x)
• n

◮ C(t) evolves according to a region competition principle

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Level set-based active regions

◮ The Chan-Vese (CV) model is a 2-phase active region model which segments Ω into object / background regions r = {in, out} with minimum luminance variance, i.e. Vr(L)(x) =

L(x)−µr 2

E(C, µ) = Eext(C, L) + α Length(C) ± p Area(Rin) it is also known as active contours without edges solutions are cartoon segmentations with minimum length and

• ptimal area

Variables

– segmenting curve C – parameters µ = (µin, µout) ≡ piecewise constant (cartoon) image fitting model

Hyperparameters: (λin, λout, α, p) > 0

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Level set-based active regions

◮ Level set formulation Rin = { x ∈ Ω | ϕ(x) > 0 } Rout = { x ∈ Ω | ϕ(x) < 0 } image energies rewrite as fixed domain integrals Ein(ϕ, L) =

• Ω

Vin(L)(x) H

ϕ(x) dx

Eout(ϕ, L) =

• Ω

Vout(L)(x)

1 − H(ϕ(x)) dx

Hε 1

2ε

δε

a smooth appproximation Hǫ (with derivative δǫ) of the Heaviside function H is used so that E is differentiable Additional hyperparameter: ǫ > 0

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Level set-based active regions

◮ E(C, µ) is minimized using an alternating minimization scheme Repeat until convergence given ϕ, optimize µ ◮ closed-form solution µr = 1 Area(Rr)

• Rr

L(x) dx given µ, optimize ϕ (n iterations) ϕt = δǫ(ϕ)

• λoutVout(L)(x)−λinVin(L)(x) + α curv(ϕ) ± p
• convergence criterion: |ϕ(t) − ϕ(t−1)|∞ < cst

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Level set-based active regions

◮ The CV model extends to 2n phases ≡ n level set functions ◮ multi-object segmentation ◮ single object segmentation involving complex object / background structure higher-dimensional / nonparametric image fitting models

– piecewise smooth – statistical mixture

◮ high-texture scene segmentation

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

The Mumford-Shah model

◮ The CV model is a special case of a nonparametric region-based segmentation model which seeks to partition Ω into regions over which luminance is best approximated (in the least square sense) by some piecewise-smooth function u E(C, u) = λ

• Ω

L(x) − u(x) 2dx

+

• Ω\C

|∇u(x)|2 dx + α Length(C)

image fidelity u piecewise-smooth C smooth

E(C, u) is known as the Mumford-Shah functional Variables

– segmenting curve C – image fitting model u ≡ image structure (i.e. non-texture) component

Hyperparameters: (λ, α) > 0

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

The Mumford-Shah model

E(C, u) = λ

• Ω

L(x)−u(x) 2dx +

• Ω\C

|∇u(x)|2 dx + α Length(C) ◮ A generic but hard problem extensive mathematical studies → many results on solution properties available yet, no closed-form solutions / no algorithms to compute solutions in the general case ◮ Extension to high-texture scene segmentation by using statistical similarity metrics as image fidelity term

Nicolas Rougon IMA4509 | Variational image segmentation

Foundations Parameterized active contours Level set methods Geodesic active contours Level set-based active regions

Variational image segmentation

Nicolas Rougon

Institut Mines-Télécom / Télécom SudParis ARTEMIS Department; CNRS UMR 8145 nicolas.rougon@telecom-sudparis.eu

May 13, 2020

Nicolas Rougon IMA4509 | Variational image segmentation