Regularization with Singular Energies

Martin Burger

Institute for Computational and Applied Mathematics European Institute for Molecular Imaging (EIMI) Center for Nonlinear Science (CeNoS)

Westfälische Wilhelms-Universität Münster

Introduction

 In many inverse and imaging problems, a priori information on the solution (or its interesting parts) needs to introduced in an effective way, e.g.

smoothness away from edges (

images

- (almost) piecewise constant (

material

)

densities with interfaces

) sparsity in some basis / frame (

MRI, ..

), or in space (

deconvolution, EEG/MEG, ..

) 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07 ¸ 2 kAu ¡ f k 2 + 1 2 kL uk 2 !

min u

Examples

 Nanoscopy with optical (4pi) techniques, cell imaging: sparsity in space 4.6.2007

©

Andreas Schönle, MPI Göttingen

Regularization with Singular Energies AIP 2007, Vancouver, June 07 ¸ 2 kAu ¡ f k 2 + 1 2 kL uk 2 !

min u

Examples

 PET / CT imaging of mouse heart coronal sagittal transverse 4.6.2007

©

Dept of Nuclear Medicine / SFB 658, WWU Münster

Regularization with Singular Energies AIP 2007, Vancouver, June 07 ¸ 2 kAu ¡ f k 2 + 1 2 kL uk 2 !

min u

Examples

 MRI images for EEG/MEG modelling T1 PD ©

Institute for BioSystemAnalysis, WWU Münster

 Anisotropic Structures essential for field simulations and dipole source reconstructions in MEG 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07 ¸ 2 kAu ¡ f k 2 + 1 2 kL uk 2 !

min u

Examples

 MRI images for EEG/MEG modelling

Baillet, Mosher, Leahy, IEEE Signal Processing Magazine, 2001, 18(6), pp. 14-30.

4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07 ¸ 2 kAu ¡ f k 2 + 1 2 kL uk 2 !

min u

Examples

 Aerial images of buildings: strong anisotropy ©

Aerowest GmbH Münster Visualization Project, Dept. of Computer Science, WWU

4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07 ¸ 2 kAu ¡ f k 2 + 1 2 kL uk 2 !

min u

Introduction

 A-priori information or desired structures have to be incorporated into reconstruction methods and regularization schemes  One possibility are singular energies (not differentiable and not strictly convex) , see also various speakers at AIP 07:

Daubechies, Candes, Fornasier, Saab, Leitao, Mizera, Kindermann, Lorenz, Ramlau, Rauhut, Zhariy, Ring, Villegas, Klann, …

4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07 ¸ 2 kAu ¡ f k 2 + 1 2 kL uk 2 !

min u

Linear Regularization

 Classical regularization schemes for inverse problems and imaging are based on linear smoothing = quadratic energy functionals  Example: Tikhonov regularization for linear

¸ 2 kAu ¡ f k

2

+ 1 2 kL uk

2

!

min

u 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07 ¸ 2 kAu ¡ f k 2 + 1 2 kL uk 2 !

min u

Linear Tikhonov

 These energy functionals are strictly convex and differentiable – standard tools from analysis and computation (Newton methods etc.) can be used  Disadvantage: possible oversmoothing , seen from first-order optimality condition 

L L u = ¡ ¸ A

¤

(Auf )

Hence

u

is in the range of

(L*L) -1 A*

4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07 ¸ 2 kAu ¡ f k 2 + 1 2 kL uk 2 !

min u

Deblurring

 Classical example: integral equation of the first kind, regularization in

L 2

(

L = Id

),

A

= Fredholm integral operator with kernel

k

u(x) = ¸ k(y; x)( ¡ k(y; z)u(z) + f (z))dy dz  Smoothness of regularized solution is determined by smoothness of kernel  For typical convolution kernels like Gaussians,

u

is analytic !

4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07 ¸ 2 kAu ¡ f k 2 + 1 2 kL uk 2 !

min u

Image Smoothing

 Classical image smoothing: data in

L 2 Id

),

L

= gradient (

H 1

-Seminorm)

¡ ¢ u + ¸ u = ¸ f

(

A =

 On a reasonable domain, standard elliptic regularity implies

u 2 H

2

( ) ,!

C( )

 Reconstruction contains no edges, blurs the image (with Green kernel) 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07 ¸ 2 kAu ¡ f k 2 + 1 2 kL uk 2 !

min u

Sparse Reconstructions ?

 Let

A

`

2

(Z)

be an operator on (basis repre sentation of a Hilbert space operator, wavelet )  Penalization by squared norm (

L = Id

)  Optimality condition for components of

u

u

k

= ¸ (A

¤

( ¡ Au + f ))

k 4.6.2007

 Decay of components determined by

A*.

Even if data are generated by sparse signal (finite number of nonzeros), reconstruction is not sparse ! Regularization with Singular Energies AIP 2007, Vancouver, June 07 ¸ 2 kAu ¡ f k 2 + 1 2 kL uk 2 !

min u

Error Estimates

 Error estimates for ill-posed problems can be obtained only under stronger conditions (source conditions)

9w : u = A

¤

w

cf. Groetsch, Engl-Hanke-Neubauer, Colton-Kress, Natterer. Nonlinear: Engl-Kunisch-Neubauer.

 Equivalent to u being minimizer of Tikhonov functional with some data  For many inverse problems unrealistic due to extreme smoothness assumptions 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07 ¸ 2 kAu ¡ f k 2 + 1 2 kL uk 2 !

min u

Error Estimates

 Condition can be weakened to

9v : u = ¾(A

¤

A)v

cf. Neubauer et al (algebraic), Hohage (logarithmic), Mathe-Pereverzyev (general).

 Advantage: more realistic conditions  Disadvantage: Estimates get worse with s 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07 ¸ 2 kAu ¡ f k 2 + 1 2 kL uk 2 !

min u

Nonlinear Problems

 Analogous (local) theory for nonlinear inverse problems (as long as forward operator is Frechet differentiable, and additional technical conditions satisfied) 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07 ¸ 2 kAu ¡ f k 2 + 1 2 kL uk 2 !

min u



Singular Energies

R r (r u)

for penalization:  

r

is strictly convex and smooth : usual smoothing behaviour / elliptic regularity 

r

is not convex : problem not well-posed  Try borderline case: singular energy 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07 ¸ 2 kAu ¡ f k 2 + 1 2 kL uk 2 !

min u

Total Variation Methods r (p) = jpj

 Simplest choice yields total variation method (

Rudin-Osher-Fatemi 89, Acar Vogel 93, Chambolle Lions 96, …)

 Singular energy: nondifferentiable, not strictly convex

Z juj

B V

„ “

Z =

g2 C 1 0

sup

;kgk 1 · 1

u(r ¢g) dx

4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Sparsity



A

operator on ` 2 (Z) \ ` 1 (Z)  Singular energy J (u) = kuk 1 = P ju k j 

¸ 2 kAu ¡ f k

2

+ J (u) !

min

u

Daubechies et al 04-07, Loubes 06/07, Ramlau, Maass, Klann 07, …

4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07 ¸ 2 kAu ¡ f k 2 + 1 2 kL uk 2 !

min u

Sparsity

 Optimality condition for components of

u

sign (u

k

) = ¸ A

¤

(f ¡ Au)

k 

A*

`

2

(Z) sign (u

k

) 2 `

2

(Z)

 Implies sparsity of u, since ksign(u)k 2 = number of nonzero component s 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07 ¸ 2 kAu ¡ f k 2 + 1 2 kL uk 2 !

min u

Sparsity / Thresholding

 For

A

being the identity , the result is well known soft-thresholding :

Donoho, Johnstone 94, Chambolle, DeVore, Lee, Lucier 98 , …

sign (u

k

) + ¸ u

k

= ¸ f

k implies

u

k

= : 8 < f f 0

k k

¡ +

¸ 1 ¸ 1

f

k

> f <

k

else ¡

¸ 1 ¸ 1 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07 ¸ 2 kAu ¡ f k 2 + 1 2 kL uk 2 !

min u



ROF Model

Z

ROF model for denoising

2 (u ¡ f )

2

+ juj

B V

!

min

u 2 B V

Rudin-Osher Fatemi 89/92, Chambolle-Lions 96, Scherzer Dobson 96, Meyer 01,…

4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

ROF Model



p + ¸ u = ¸ f ; p 2 @juj

B V  Dual variable p enters – related to mean curvature of edges for total variation 

@J (u) = f p 2 X j 8v 2 X : J (u) + hp; v ¡ ui · J (v)g

4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

ROF Model clean Reconstruction

(code by Jinjun Xu)

noisy ROF 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

ROF Model

 ROF model denoises cartoon images resp. computes the cartoon of an arbitrary image, natural spatial multi-scale decomposition by varying l 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Numerical Differentiation with TV

Bachmayr, 2007

4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Singular Energies

 Methods with singular energies have great potential, but still some problems: difficult to analyze estimates and to obtain error systematic errors (like loss of contrast) computational challenges - strong bias – how to incorporate uncertain a-priori structures (adaptively) ?

4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Systematic Errors and Bias

Berkels, mb, Droske, Nemitz, Rumpf 06

4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Loss of Contrast

 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

4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Iterative Refinement

 Idea: add the residual („noise“) back to the image to pronounce the features decreased procedure u k = arg min u · ¸ 2 Z (u ¡ f ¡ v k Iterative ¡ 1 ) 2 + juj B V ¸ v k = v k ¡ 1 + (f ¡ u k ); v 0 = 0

Osher-mb-Goldfarb-Xu-Yin 04

4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Iterative Refinement & ISS

 Improves reconstructions significantly 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Iterative Refinement & ISS

4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Iterative Refinement

 u Works for inverse problems k = arg min u · ¸ 2 kAu in a similar way ¡ f ¡ v k ¡ 1 k 2 + J (u) ¸ v k = v k ¡ 1 + (f Au k ); v 0 = 0

Osher-mb-Goldfarb-Xu-Yin 04

4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Iterative Refinement & ISS  Observation from optimality condition p k = ¸ A ¤ ( ¡ Au k + f + v k ¡ 1 ) 2 @J (u k )  residual and subgradient p k = ¸ A ¤ v k 4.6.2007

 Iterates determined by equivalent ¸ = minimization of 2 kAu ¡ f k ¸ 2 kAu 2 + J (u) J (u ¡ f k 2 k ¡ 1 + J (u) J (u k ) ¡ 1 ¡ ¸ hv ) k ¡ hp ¡ 1 k ; Au Au ¡ 1 ; u u k k ¡ 1 i ¡ 1 i Regularization with Singular Energies AIP 2007, Vancouver, June 07

Iterative Refinement & ISS  Dual update formula p k = p k ¡ 1 + ¸ A ¤ ( ¡ Au k + f )  Iterative refinement = dual proximal method = Bregman iteration. Minimization in each step of 2 kAu ¡ f k 2 + D p k J ¡ 1 (u; u k ¡ 1 )  D J (u; v) = J (u) J (v) q 2 @J (v) ¡ hq; u vi 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Iterative Refinement & ISS

 Choice of parameter l less important, can be kept small (oversmoothing). Regularizing effect comes from appropriate stopping.  Quantitative stopping rules available, or „stop when you are happy“ in imaging – S.O.

 Limit l to zero can be studied. Yields gradient flow for the dual variable („inverse @ t p(t) = A ¤ ( ¡ Au(t) + f ); p 2 @J (u)

mb-Gilboa-Osher-Xu 06, mb-Frick-Osher-Scherzer 06

4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Iterative Refinement & ISS

 Efficient numerical schemes for flow

mb-Gilboa-Osher-Xu 06

 Analysis of iteration and partly of the flow

Osher-mb-Goldfarb-Xu-Yin 04, mb-Frick-Osher Scherzer 06

 Error estimates in Bregman distance

mb-He-Resmerita 07

4.6.2007

 Non-quadratic fidelity is possible, some caution needed for

L 1

fidelity

He-mb-Osher 05, mb-Frick-Osher-Scherzer 06

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Iterative Refinement & ISS

 Application to other energies, e.g. Besov norms (wavelets), is straight-forward  Starting from soft shrinkage, iterated refinement yields firm shrinkage , inverse scale space becomes hard shrinkage

Osher-Xu 06

 Bregman is distance natural sparsity measure , number of nonzero components is constant in error estimates 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Surface Smoothing

 Smoothing of surfaces in level set represenation  3D Ultrasound, Kretz / GE Med.

4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Inverse Scale Space

4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

MRI Reconstruction

MRI Data Siemens Magnetom Avanto 1.5 T Scanner

He, Chang, Osher, Fang, Speier 06

Penalization TV + Besov 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

MRI Reconstruction

MRI Data Siemens Magnetom Avanto 1.5 T Scanner

He, Chang, Osher, Fang, Speier 06

Penalization TV + Besov 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Iterative Refinement & ISS

 Generalization to nonlinear inverse problems possible

Bachmayr, Thesis 07

 Different ways of approximating nonlinearity lead to different iterative schemes (similar to iterated Tikhonov / Landweber / Levenberg-Marquardt)  Example: parameter identification (diffusivity) in elliptic PDE 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Iterative Refinement & ISS

Exact Solution Reconstructions at 1 % noise Iterative 1 Iterative 2 Standard TV 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Current / Future Work: EM-TV

 Generalization to combination with EM / Richardson-Lucy method in progress A.Sawatzky, C.Brune

 Application 1: 4pi / STED nanoscopy  Application 2: PET/CT imaging with O 15 isotopes (fast decay, hence bad statistics) 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Adaptive Bias / Parametrization

 Bias of one functional often too strong J (u; ®)  parametrized by  Example: adaptive anisotropy in total variation methods 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Adaptive Anisotropy

 In aerial images the typical anisotropy is rectangular, houses have 90 ° angles  But not all of them have the same orientation 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Adaptive Anisotropy

 Bias for edges with 90 ° angles from J (u; ®) = (jv 1 j + jv 2 j) dx; v = R ® r u 

R

a is rotation matrix for angle the orientation a to capture  Since orientation is not constant over the image, a has to vary and to be found adaptively by minimization 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Adaptive Anisotropy

 To avoid microstructure, variation of a has to be regularized, too  ¸ 2 Possible functional to be minimized (u ¡ f ) 2 + J (u; ®) + ¹ 2 1 Z jr ®j 2 + ¹ 2 2 Z j ¢ ®j 2 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Adaptive Anisotropy

 Improves angles, still loses contrast 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Adaptive Anisotropy

 Contrast correction again by iterative refinement  Angle parameter provides classification of orientations in the image 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Adaptive Anisotropy

 Cartoon reconstruction and orientational classification of aerial images 4.6.2007

Berkels, mb, Droske, Nemitz, Rumpf 06

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Adaptive Anisotropy

 Cartoon reconstruction and orientational classification of aerial images

Berkels, mb, Droske, Nemitz, Rumpf 06

4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Adaptive Anisotropy

 Analogous problem in segmentation of MRI brain images for EEG/MEG  Regularization by total variation (= length of curves in segmentation) kills noise, but also elongated structures  Adapt anisotropy (locally like sharp ellipse) to find sulci accurately and provide classification of normals (for dipole fitting, source reconstruction) 4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Adaptive Anisotropy

Denis Neiter, results of internship 2007

4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07

Download and Contact

Papers and talks at www.math.uni-muenster.de/u/burger or by email [email protected] Thanks for input and suggestions to:

S.Osher, J.Xu, G.Gilboa, D.Goldfarb, W.Yin, L.He, E.Resmerita, M.Bachmayr, B.Berkels, M.Droske, O.Nemitz, M.Rumpf, K.Frick, O.Scherzer, A.Schönle, T.Hohage, C.Wolters, T.Kösters, K.Schäfers, F.Wübbeling, A.Sawatzky, D.Neiter

4.6.2007

Regularization with Singular Energies AIP 2007, Vancouver, June 07