Reconstruction de volumes ` a partir de coupes Simon Masnou - PowerPoint PPT Presentation

Reconstruction de volumes ` a partir de coupes Simon Masnou Institut Camille Jordan Universit´ e Lyon 1 S´ eminaire Parisien des Math´ ematiques Appliqu´ ees ` a l’Imagerie Institut Henri Poincar´ e 2 f´ evrier 2017 en collaboration avec Elie Bretin (INSA Lyon) et Fran¸ cois Dayrens (ENS Lyon)

Motivation Frequent problem in medical imaging (MRI, CT) and computational geometry: how to reconstruct a volume from a few slices (or more generally from partial data)?

Motivation

Formulation with inner / outer constraints 0.5 0.4 0.3 0.2 0.1 0 −0.1 −0.2 −0.3 −0.4 −0.5 −0.5 0 0.5

Formulation with inner / outer constraints What is a ”good shape” satisfying the constraints?

Geometric optimization in real life

Modeling Let ω int , ω ext ⊂ R N Geometric optimization problem J ( E ) | ω int ⊂ E ⊂ R N � ω ext � � inf where J is a geometric energy ◮ Natural choice: J =perimeter or Willmore energy ◮ A natural topology is the L 1 topology of characteristic functions of sets ◮ The problem is however ill-posed (at least for the perimeter) when | ω int | = | ω ext | = 0.

Perimeter in the BV sense Perimeter E has finite perimeter if its characteristic function ✶ E ∈ BV Denote P ( E ) = TV ( ✶ E ) its perimeter. The perimeter functional is lower semicontinuous for the L 1 topology. P ( E ) = length ( ∂ E ) P ( E ) = area ( ∂ E )

Perimeter in the BV sense Perimeter E has finite perimeter if its characteristic function ✶ E ∈ BV Denote P ( E ) = TV ( ✶ E ) its perimeter. The perimeter functional is lower semicontinuous for the L 1 topology. P ( E ) � = length ( ∂ E ) P ( E ) � = length ( ∂ E )

Natural formulation of the reconstruction problem for the perimeter P ( E ) | ω int ⊂ E 1 , ω ext ⊂ E 0 � � inf 0.5 0.4 0.3 0.2 0.1 0 −0.1 −0.2 −0.3 −0.4 −0.5 −0.5 0 0.5

Bernoulli-Euler elastic energy in R 2 Curvature × × Let γ be a C 2 curve in R 2 , κ = det( γ ′′ , γ ′ ) | γ ′ | 3 1 /κ × × Bernoulli-Euler energy Let E be a set with C 2 boundary, � κ 2 d H 1 . W ( E ) = ∂ E γ

Willmore energy (in R 3 ) Mean curvature Let M = C 2 surface in R 3 , H = 1 2( κ 1 + κ 2 ) κ 1 , κ 2 : principal curvatures Willmore energy If E has C 2 boundary ∂ E , � H 2 d H 2 . W ( E ) = image credits: Wikipedia ∂ E

Natural formulation of the reconstruction problem for the Willmore energy The Willmore energy is not lower semicontinuous in L 1 . For minimization purposes, use its relaxation W (i.e. its lower semicontinuous envelope). We address the following problem: W ( E ) | ω int ⊂ E 1 , ω ext ⊂ E 0 � � inf

Approximation of the problem I: Perimeter approximation � ε |∇ u ε | 2 d x ≈ 1 ε Area ≈ 1 Thus, εε P ( E ) = P ( E ) as ε → 0. However, any constant function has zero energy! How to force u ε to be close to a characteristic function, i.e. a binary function?

Perimeter approximation Use a double-well potential, for instance G ( s ) = 1 2 s 2 (1 − s ) 2 . G 0 1 �� 1 � If sup ε G ( u ε ) d x < + ∞ then u ε → 0 or 1 a.e. as ε → 0. ε Therefore, u ε approximates a characteristic function.

The Cahn-Hilliard functional (Van der Waals)-Cahn-Hilliard energy The phase-field approximation of perimeter is given by � � ε 2 |∇ u | 2 + 1 � P ε ( u ) = ε G ( u ) d x ε u ε ✶ E E E where G is a double-well potential. G e.g., G ( s ) = 1 2 s 2 (1 − s ) 2 0 1

Phase-field approximation of perimeter Convergence of P ε (Modica, Mortola - 1977) P ε converges to � λ P ( E ) si u = ✶ E ∈ BV P ( u ) = + ∞ otherwise in the sense of Γ-convergence where λ is a constant depending only on potential G . Property of Γ-convergence Let X be a metric space and ( F ε ) a sequence of equicoercive functionals converging to F in the sense of Γ-convergence in X . If u ε is a minimizer of F ε , then there exists a minimizer u of F , s.t. u ε → u .

Optimal profile One can define the phase-field optimal profile associated with E : � 1 � q ( s ) = 1 2(1 − tanh( s u ε ( x ) = q ε d s ( x , E ) with 2)) Signed distance q 1 d s ( x , E ) = d ( x , E ) − d ( x , R N � E ) 0 Convergences For a bounded set E ◮ u ε → ✶ E ◮ P ε ( u ε ) → λ P ( E ) if E has finite perimeter as ε → 0.

Phase field approximation of the Willmore energy The L 2 -gradient of P ε satisfies −∇ L 2 P ε ( u ) = ε ∆ u − 1 ε G ′ ( u ) . The gradient flow of perimeter is the mean curvature flow and −∇ L 2 P ε ( u ε ) approximates the mean curvature of ∂ E in the transition zone of u ε when u ε ≈ ✶ E . Approximation of the Willmore energy In R 2 and R 3 , the energy � 2 1 � ε ∆ u − 1 � ε G ′ ( u ) u �→ P ε ( u ) + W ε ( u ) = P ε ( u ) + d x 2 ε Γ-converges to E �→ λ ( P ( E ) + W ( E )) if E is C 2 and compact ◮ De Giorgi + Bellettini, Paolini (1993) + R¨ oger, Sch¨ atzle (2006)

Optimal profile With the same phase-field profile associated with E � 1 � u ε ( x ) = q ε d s ( x , E ) one has Convergences For a bounded set E ◮ u ε → ✶ E ◮ P ε ( u ε ) → λ P ( E ) if E has finite perimeter ◮ W ε ( u ε ) → λ W ( E ) if ∂ E is C 2 as ε → 0.

Inclusion-exclusion constraints Let ω int , ω ext ⊂ R N Geometric optimization problem inf { J ( E ) | ω int ⊂ E 1 , ω ext ⊂ E 0 } where J is either P , or W One defines obstacle constraints: � 1 � � 1 � u int ε d s ( x , ω int ) u ext ε d s ( x , ω ext ) ε ( x ) = q and ε ( x ) = 1 − q Key property ω int ⊂ E ⊂ R N � ω ext u int � u ε � u ext ⇐ ⇒ ε ε In the phase field approximation, constraints can be interpreted as a linear obstacle problem!

Numerical scheme for perimeter Approximating a solution to min { P ε ( u ) | u int � u � u ext ε } ε ◮ Initialize u 0 ; ◮ At step n , given u n , use a splitting method: ◮ u n +1 / 2 is obtained by one step of an implicit discrete L 2 gradient flow for P ε , i.e. u n +1 / 2 − u n = δ t ( ε ∆ u n +1 / 2 − 1 ε G ′ ( u n +1 / 2 ) (discrete Allen-Cahn equation) ◮ Get u n +1 from u n +1 / 2 by projecting onto the constraints u int � u � u ext ε ε

Implicit discrete gradient flow Finding u n +1 / 2 is equivalent to finding a fixed point of the map: �� �� u n + δ t v �→ ( I d − δ t ε ∆) − 1 ε G ′ ( v ) . Picard iterations give a stable scheme, and solving in Fourier domain provides an excellent spatial accuracy

Matlab code (projection is embedded into the fixed point scheme)

Numerical scheme for Willmore Same principle, but now the L 2 flow is: � ∂ t v = ∆ µ − 1 ε 2 G ′′ ( v ) µ, µ = 1 ε G ′ ( v ) − ε ∆ v , It can be discretized at step n as u n +1 / 2 = u n + δ t ∆ µ n +1 / 2 − 1 � � ε 2 G ′′ ( u n +1 / 2 ) µ n +1 / 2 � µ n +1 / 2 = 1 ε G ′ ( u n +1 / 2 ) − ε ∆ u n +1 / 2 . whose solution ( u n +1 / 2 , µ n +1 / 2 ) is a fixed point of the map: � I d � − 1 � u n − δ t � − δ t ∆ ε G ′′ ( u ) µ v �→ , 1 ε ∆ ε 2 G ′ ( u ) I d Again, an efficient and accurate scheme can be designed using Fourier transform.

First experiments 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 −0.1 −0.1 −0.2 −0.2 −0.3 −0.3 −0.4 −0.4 −0.5 −0.5 −0.5 0 0.5 −0.5 0 0.5 Willmore energy Perimeter

Interpretation In some cases, the energy � if u int � u ε � u ext P ε ( u ) ε ε P 1 ,ε ( u ε ) = + ∞ otherwise converges to F 1 ( u ) = λ ( P ( E ) + H ( E )) if u = ✶ E , and u ε → u as ε → 0 ◮ H ( E ) = length (in 2D) or area (in 3D) of the set ( E 0 ∩ ω int ) � ( E 1 ∩ ω out ) However, the term λ H ( E ) may favor constraints violation

Constraints violation A situation where violating the outer constraint is more favorable for F 1 : F 1 (left configuration) < F 1 (right configuration) In contrast, defining F 2 ( E ) = λ ( P ( E ) + 2 H ( E )) one has F 2 (left configuration) > F 2 (right configuration)

Remedy: use fat constraints Thicken the constraints to give them volume √ ε Ω int ω int ε � 1 The function U int ε d s ( x , Ω int � ε ( x ) = q ε ) takes values in [0 , 1]. ω int U int = 1 ε

Convergence result Theorem (Bretin, Dayrens, M.) The energy � si U int � u � U ext P ε ( u ) ε ε P 2 ,ε ( u ) = + ∞ sinon Γ-converges to F 2 ( u ) = λ ( P ( E ) + 2 H ( E )) if u = ✶ E , as ε → 0 ◮ Minimal sets for F 2 satisfy inclusion-exclusion constraints in reasonable cases. ◮ Characterizing the Γ-limit for the Willmore energy is delicate due to the non locality and ghost parts.

Numerical experiments 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 −0.1 −0.1 −0.2 −0.2 −0.3 −0.3 −0.4 −0.4 −0.5 −0.5 −0.5 0 0.5 −0.5 0 0.5 ”Thin” perimeter ”Fat” perimeter

3D reconstruction with Willmore energy Reconstruction of a 3D brain image from real MRI slices

Confined elastica (i.e. a minimizer of the constrained Bernoulli-Euler energy) An elastica in a fox head

Can be used for smoothing pixellized surfaces cf Bretin, Lachaud, Oudet, 2011 where was used a penalization of the constraints violation set volume (acting as a repulsion force)

Other ”slices”