Martin Burger Institute for Computational and Applied Mathematics - - PowerPoint PPT Presentation

Martin Burger

Institute for Computational and Applied Mathematics

(4D) Variational Models Preserving Sharp Edges

Martin Burger

Mathematical Imaging Workgroup @WWU

2

Linz, 2011

Martin Burger

3

Some Philosophy

„No matter what question, L1 is the answer“ Stanley O. Regularization in data assimilation is at the same state it was 10 years ago in biomedical imaging The understanding and methods we gained in medical imaging can hopefully be useful in geosciences and data assimilation

Martin Burger

4

Biomedical Imaging: 2000 vs 2010

Modality State of the art 2000 State of the art 2010

Full CT Filtered Backprojection Exact Reconstruction PET/SPECT Filtered Backprojection /EM EM-TV / Dynamic Sparse PET-CT

• EM-AnatomicalTV

Acousto-Opt.

• Wavelet Sparse / TV

EEG/MEG LORETA Sparsity / Bayesian ECG-BSPM Least Norm L1 of normal derivative Microscopy None, linear Filter Poisson-TV / Shearlet-L1

Martin Burger

Based on joint work with

Martin Benning, Michael Möller, Felix Lucka, Jahn Müller (Münster) Stanley Osher (UCLA) Christoph Brune (Münster / UCLA / Vancouver) Fabian Lenz (Münster), Silvia Comelli (Milano/Münster) Eldad Haber (Vancouver) Mohammad Dawood, Klaus Schäfers (NucMed/EIMI Münster)

5

Linz, 2011

SFB 656

Martin Burger

6

Regularization of Inverse Problems

We want to solve Forward operator between Banach spaces with finite dimensional approximation (sampling, averaging)

Martin Burger

Dynamic Biomedical Imaging

7

Saarbrücken, 9.7.10

Maximum Likelihood / Bayes

Reconstruct maximum-likelihood estimate Model of posterior probability (Bayes) Yields regularized variational problem for finite m

Martin Burger

8

Minimization of penalized log-likelihood

General variational approach Combines nonlocal part (including K ) with local regularization functional Gaussian noise (note: covariance hidden in output norm)

Martin Burger

9

Example Gauss:

Additive noise, i.i.d. on each pixel, mean zero, variance s Minimization of negative posterior log-likelihood yields Asymptotic variational model

Martin Burger

10

Optimality

Existence and uniqueness by variational methods General case: optimality condition is a nonlinear integro-differential equation / inclusion (integral operator K, differential operator in J ) Gauss:

Martin Burger

11

Robustness

Due to noisy data robustness of with respect to errors in f is important Problem is robust for large a, but data are only reproduced for small a Convergence of solutions as f converges or as a to zero in weak* topology

Martin Burger

12

Structure of Solutions

Analysis by convex optimization techniques, duality Structure of subgradients important Possible solution satisfy source condition Allows to gain information about regularity (e.g. of edges)

Martin Burger

13

Structure of Solutions

Optimality condition for Structure of u determined completely by properties of uB and K* For smoothing operators K, singularity not present in uB cannot be detected Model error goes into K resp. K* and directly modifies u

Martin Burger

4D VAR

Given time dynamics starting from unknown initial value Variational Problem to estimate initial state for further prediction

14

Linz, 2011

Martin Burger

4D VAR = 3D Variational Problem

Elimination of further states from dynamics Effective Variational Problem for initial value in 3D

15

Linz, 2011

Martin Burger

Example: Linear Advection

Minimize quadratic fidelity + TV of initial value subject to Upwind discretization

16

Linz, 2011

Martin Burger

4D VAR for Linear Advection

Gibbs phenomenon as usual

17

Linz, 2011

Martin Burger

4D VAR for Linear Advection

Full observations (black), noisy(blue), 40 noisy samples (red)

18

Linz, 2011

Martin Burger

4D VAR for Linear Advection

Different noise variances

19

Linz, 2011

Martin Burger

Analysis of Model Error

Optimality Exact Operator for linear advection is almost unitary Hence

20

Linz, 2011

Martin Burger

21

Beyond Gaussian Priors

Again: optimality condition for MAP estimate If J is strictly convex and smooth, subdifferential is a singleton containing only the gradient of J, which can be inverted to

• btain a similar relation. Again operator determines structure

Only chance to obtain full robustness: multivalued

• subdifferential. Singular regularization

Martin Burger

22

Singular Regularization

Construct J such that the subdifferential at points you want to be robust is large Example: l1 sparsity Zeros are robust

Martin Burger

23

TV-Methods: Structural Prior (Cartooning)

Penalization of total Variation Formal Exact ROF-Model for denoising g : minimize total variation subject to

Rudin-Osher-Fatemi 89,92

Martin Burger

24

Why TV-Methods ?

Cartooning Linear Filter TV-Method

Martin Burger

ROF Model

clean noisy ROF

Martin Burger

26

H2O15 PET – Left Ventricular Time Frame

EM EM-Gauss EM-TV

Martin Burger

Dynamic Biomedical Imaging

27

Saarbrücken, 9.7.10

H2O15 PET – Right Ventricular Time Frame

EM EM-Gauss EM-TV

Martin Burger

4D VAR for Linear Advection

Gibbs phenomenon as usual

28

Linz, 2011

Martin Burger

4D VAR for Linear Advection

Full observations (black), noisy(blue), 40 noisy samples (red)

29

Linz, 2011

Martin Burger

4D VAR TV for Linear Advection

Comparison for full observations

30

Linz, 2011

Martin Burger

4D VAR TV for Linear Advection

Comparison for observed samples

31

Linz, 2011

Martin Burger

4D VAR TV for Linear Advection

Comparison for observed samples with noise

32

Linz, 2011

Martin Burger

Analysis of Model Error

Variational problem as before, add Optimality condition As before

33

Linz, 2011

Martin Burger

Analysis of Model Error

Structures are robust: apply T in region where If we find s solving Poisson equation with then

34

Linz, 2011

Martin Burger

Numerical Solution: Splitting or ALM

Operator Splitting into standard problem (dependent on code) and simple denoising-type problem Example: Peaceman Rachford-Splitting for

35

Linz, 2011

Martin Burger

36

Bayes and Uncertainty

Natural prior probabilities for singular regularizations can be constructed even in a Gaussian framework Interpret J(u) as a random variable with variance s2 Prior probability density MAP estimate minimizes

Martin Burger

37

Bayes and Uncertainty

Equivalence to original form via constraint regularization For appropriate choice of a and g, minimization of and is equivalent to subject to

Martin Burger

38

Uncertainty Quantification

Sampling with standard MCMC schemes difficult Novel Gibbs sampler by F.Lucka based on analytical integration of posterior distribution function in 1D Theoretical Insight: MSc Thesis Silvia Comelli CM Estimate for TV prior

Martin Burger

39

Uncertainty Quantification II

Error estimates in dependence on the noise, using source conditions Error estimates need appropriate distance measure,generalized Bregman-distance

mb-Osher 04, Resmerita 05, mb-Resmerita-He 07, Benning-mb 09

Estimates for Bayesian distributions in Bregman transport distances (w. H.Pikkarainen) = 2 Wasserstein distance in the Gaussian case

Martin Burger

40

Uncertainty Quantification III

Idea: construct linear functionals from nonlinear eigenvectors We have For TV-denoising (also for linear advection example), Estimate of maximal error for mean value on balls For l1-sparsity estimate of error in single components

Benning PhD 11, Benning-mb 11

Martin Burger

ROF minimization loses contrast, total variation of the reconstruction is smaller than total variation of clean image. Image features left in residual f-u

g, clean

f, noisy u, ROF f-u

mb-Gilboa-Osher-Xu 06

Loss of Contrast

Martin Burger

42

Loss of Contrast = Systematic Bias of TV

Becomes more severe in ill-posed problems with operator K Not just simple vision effect to be corrected, but loss of information Simple idea for Least-Squares: add back the noise to amplify = Augmented Lagrangian Osher-mb-Goldfarb-Xu-Yin 2005

Martin Burger

43

Bregman Iteration

Can be shown to be equivalent to Bregman iteration Immediate generalization to convex fidelities and regularizers Generalization to Gauss-Newton type Methods for nonlinear K: use linearization of K around last iterate ul Bachmayr-mb 2009

Martin Burger

44

Bregman Iteration

Properties like iterative regularization method Regularizing effect from appropriate termination of the iteration Better performance for oversmoothing single steps, i.e. regularization parameter a very large Limit: Inverse Scale Space Method

mb-Gilboa.Osher-Xu 2006

Martin Burger

45

Why does Inverse Scale Space work ?

Singular value decomposition in fully quadratic case Eigenfunctions: yields Convergence faster in small frequencies (large eigenvalues)

Martin Burger

46

Why does Inverse Scale Space work ?

Convex one-homogeneous regularization J (TV, l1, …) Eigenfunctions: yields Again large frequencies appear later. Not at all for small t ! Eigenvalues in TV indeed related to jump measures PhD-Thesis Benning, 2011

Martin Burger

47

Why does Inverse Scale Space work ?

Multiple frequencies not simple for nonlinear case However, various theoretical and computational results confirming exact scale decomposition PhD-Thesis Benning, 2011 / mb-Frick-Scherzer-Osher 2007 Complete characterization of inverse scale space for discrete l1-functionals, yields jump dynamics in time, adaptive basis pursuit method with guaranteed convergence mb-Möller-Benning-Osher, 2011

Martin Burger

48

Saarbrücken, 9.7.10

18F-FDG

PET

EM, 20 min EM-TV, 5s EM, 5s BREG, 5s

Jahn Müller, 2011 Data from Nuclear Medicine Department, UKM

Martin Burger

49

STED Microscopy

Christoph Brune, 2009 Data from MPI for

• Biophys. Chem.

Göttingen (K.Willig, A.Schönle, Hell)

Martin Burger

4D Reconstruction

4D imaging of transport with penalization of large velocities: Minimize subject to

50

Linz, 2011

Martin Burger

51

Linz, 2011

Analysis of Motion Model

Functional related to Benamou-Brenier formulation of optimal transport . Analysis different from optimal transport, since usually no initial and final densities are given (more related to mean-field games, Lasry-Lions 07) Existence by transformation to

• A-priori estimate for w in L2. Weak compactness
• A-priori estimates for u in Lp(0,T;BV) and for time derivative in

Lq(0,T;W-1,s)

• Adaptation of Aubin-Lions gives strong compactness of u in

Lr(0,T; Lr), and thus of the square-root in L2r(0,T; L2r)

Martin Burger

52

Linz, 2011

4D TV Model

Analysis relies on superlinear growth of F, although F=Identity seems a very reasonable choice Choosing F equal to the identity would imply we seek a minimal L1 norm of the vector of total variations. Favours sparsity, i.e. solutions with very large total variation at some time step allowed if small else. This does not correspond to a smooth motion model, hence superlinear choices preferable Some indications of this effect in numerical results

Martin Burger

53

Linz, 2011

Numerical solution

Complicated 4D variational problem combining various integral and differential operators + nonlinearity. Convexity achieved by formulation in momentum variable m = u V Efficient GPU implementation by Christoph Brune on CUDA with specially designed algorithms. All subproblems solvable by FFT or shrinkage Realized by introducing new variables and inexact Uzawa Augmented Lagrangian approach

Martin Burger

54

Linz, 2011

Augmented Lagrangian

Martin Burger

55

Linz, 2011

Inexact Uzawa Augmented Lagrangian

Martin Burger

56

Linz, 2011

Update of Primal Variables

Martin Burger

57

Linz, 2011

Results: Deblurring, Synthetic Data

Exact solution Blurred Data

Martin Burger

58

Linz, 2011

Results: Deblurring, Synthetic Data

Exact solution Reconstruction

Martin Burger

Results: Cardiac 18F-FDG PET (Eulerian)

PET Reconstruction (Data) Registration to Diastole Registration to Systole

59

Linz, 2011

Martin Burger

Info

http://imaging.uni-muenster.de http://www.cells-in-motion.de http://www.herzforscher.de

60

Linz, 2011