Bregman-EM-TV Methods for Deconvolution Problems with Poisson Noise - - PowerPoint PPT Presentation

bregman em tv methods for deconvolution problems with
SMART_READER_LITE
LIVE PREVIEW

Bregman-EM-TV Methods for Deconvolution Problems with Poisson Noise - - PowerPoint PPT Presentation

Introduction EM-TV Methods Numerical Results Bregman-EM-TV Methods for Deconvolution Problems with Poisson Noise Christoph Brune, Alex Sawatzky Prof. Martin Burger Institute of Computational and Applied Mathematics University of M unster,


slide-1
SLIDE 1

Introduction EM-TV Methods Numerical Results

Bregman-EM-TV Methods for Deconvolution Problems with Poisson Noise

Christoph Brune, Alex Sawatzky

  • Prof. Martin Burger

Institute of Computational and Applied Mathematics University of M¨ unster, Germany Conference on Applied Inverse Problems, Vienna 2009

1 Christoph Brune University of M¨ unster

slide-2
SLIDE 2

Introduction EM-TV Methods Numerical Results

Outline

1

Introduction Inverse Problems in Microscopy Goal / Motivation

2

EM-TV Methods EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

3

Numerical Results

2 Christoph Brune University of M¨ unster

slide-3
SLIDE 3

Introduction EM-TV Methods Numerical Results Inverse Problems in Microscopy Goal / Motivation

STED-Microscope

Stimulated Emission Depletion (Prof. Hell, MPI G¨

  • ttingen)

Fluorescence by focused laser beam; Stimulation + Depletion

3 Christoph Brune University of M¨ unster

slide-4
SLIDE 4

Introduction EM-TV Methods Numerical Results Inverse Problems in Microscopy Goal / Motivation

Optical Nanoscopy in so far unknown sharpness

STED microscope overcomes physical limits imposed by diffraction Correlation between laser intensity, convolution kernel, granularity Industry Partner: Leica Microsystems CMS GmbH

4 Christoph Brune University of M¨ unster

slide-5
SLIDE 5

Introduction EM-TV Methods Numerical Results Inverse Problems in Microscopy Goal / Motivation

Image Reconstruction

Basic problem: We consider the inverse problem, a linear equation K˜ u + b = g K : L1(Ω) → L1(Σ) linear and compact operator with non-closed range and obtaining positivity ill-posed problem (Hadamard) ˜ u exact image (desired solution) b assumed constant background model g exact data (values before detection) Typical examples are Fredholm integral operators of the first kind (Ku)(x) =

k(x, y)u(y) dy , x ∈ Σ , k non-negative .

5 Christoph Brune University of M¨ unster

slide-6
SLIDE 6

Introduction EM-TV Methods Numerical Results Inverse Problems in Microscopy Goal / Motivation

Reconstruction in Optical Nanoscopy

Deconvolution problems with Poisson noise Ku + b = f (Ku)(x) =

h(x − τ)u(τ)dτ Nanoscopy provides blurred images Statistical noise model, described by likelihood p(f |u) → fi realization of random var. Fi Poisson-noise resulting from laser sampling (photon counts)

exact data ˜ u kernel h (psf) convolved data convolved & noisy data f

6 Christoph Brune University of M¨ unster

slide-7
SLIDE 7

Introduction EM-TV Methods Numerical Results Inverse Problems in Microscopy Goal / Motivation

Reconstruction in Optical Nanoscopy

Deconvolution problems with Poisson noise Ku + b = f (Ku)(x) =

h(x − τ)u(τ)dτ Nanoscopy provides blurred images Statistical noise model, described by likelihood p(f |u) → fi realization of random var. Fi Poisson-noise resulting from laser sampling (photon counts)

exact data ˜ u kernel h (psf) convolved data convolved & noisy data f

6 Christoph Brune University of M¨ unster

slide-8
SLIDE 8

Introduction EM-TV Methods Numerical Results Inverse Problems in Microscopy Goal / Motivation

Goal / Motivation

3D image reconstruction method Accurate, fast and robust deblurring and noise removal Cartoon reconstructions → segmentation, quantification Well-posedness of a continuous variational model Positivity preserving algorithm Proved convergence of algorithm Simultaneous contrast enhancement

7 Christoph Brune University of M¨ unster

slide-9
SLIDE 9

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

Expectation Maximization (EM)

Statistical modeling, Bayes theory Formulation as a random experiment Fi = Poiss{(Ku)i + b} EM-Algorithm by Dempster, Laird, Rubin (E)step: Averaging complete data log-likelihood over its cond. distr. (M)step: Maximization of conditional expectation Use of maximum a-posteriori (MAP) estimation wanted: ˆ u = arg max

u

p(u|f )

8 Christoph Brune University of M¨ unster

slide-10
SLIDE 10

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

MAP Likelihood Approach

Consider the a-posteriori likelihood: p(u | f ) = p(f | u) p(u) p(f ) , p(u) ∼ e−αR(u) (Gibbs)

i.i.d.

=

  • x∈Σ

pλ(F = f ) · p(u) , λ = E [F] = Ku =

  • x∈Σ

Ku(x)f (x) f (x)! e−Ku(x) · e−αR(u) = ⇒ ˆ u = arg max

u

p(u | f ) = arg min

u {− log p(u | f )}

= arg min

u {

  • Σ

Ku − f log Ku + f log f − f

  • independent of u

+ α R(u)} Kullback-Leibler divergence Regularization

9 Christoph Brune University of M¨ unster

slide-11
SLIDE 11

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

Variational Formulation for R(u) = 0

min

u

  • Σ

Ku − f log Ku Optimality condition K ∗1 − K ∗ f Ku

  • = 0,

(Scal. K ∗1 = 1) Leads to fixed point equation u = u K ∗ f Ku

  • Iteration scheme

uk+ 1

2 = ukK ∗

f Kuk

  • (EM step)

10 Christoph Brune University of M¨ unster

slide-12
SLIDE 12

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

Variational Formulation for R(u) = 0

min

u ≥ 0

  • Σ

Ku − f log Ku Optimality condition / KKT : ∃λ ≥ 0 λ u = 0, K ∗1 − K ∗ f Ku

  • −λ = 0,

(Scal. K ∗1 = 1) Leads to fixed point equation u = u K ∗ f Ku

  • Iteration scheme

uk+ 1

2 = ukK ∗

f Kuk

  • (EM step)

10 Christoph Brune University of M¨ unster

slide-13
SLIDE 13

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

Variational Formulation for R(u) = 0

min

u ≥ 0

  • Σ

Ku − f log Ku Optimality condition / KKT : ∃λ ≥ 0 λ u = 0, K ∗1 − K ∗ f Ku

  • −λ = 0,

(Scal. K ∗1 = 1) Leads to fixed point equation u = u K ∗ f Ku

  • Iteration scheme

uk+ 1

2 = ukK ∗

f Kuk

  • (EM step)

Convergence analysis for exact data: descent in objective functional

10 Christoph Brune University of M¨ unster

slide-14
SLIDE 14

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

Example: Nanoscopy

Convolved + Poisson noise EM, KL-Distance: 3.20

11 Christoph Brune University of M¨ unster

slide-15
SLIDE 15

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

Motivation for TV Regularization

Need additional (a-priori) information about:

Known structures in the image Desired structures to be investigated

EM algorithm uses uniform prior probability distribution Prior probability can be related to regularization functional p(u) ∼ e−α|u|BV TV-Regularization to preserve edges R(u) := |u|BV = sup

g∈C ∞

0 (Ω;Rd)

||g||∞≤1

u divg ( =

|∇u| )

12 Christoph Brune University of M¨ unster

slide-16
SLIDE 16

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

EM-Algorithm and TV-Regularization

min

u∈BV (Ω)

  • Σ

Ku − f log Ku + α |u|BV , α > 0 Combining nonlocal part (incl. K) and local regularization functional Convex, but not differentiable Optimality condition for u > 0 K ∗1 − K ∗( f Ku ) + α p = 0, p ∈ ∂|u|BV K ∗1 − K ∗( f Ku ) − α div( ∇u |∇u|) = 0 formal Assume K ∗1 = 1 in the following (operator scaling)

13 Christoph Brune University of M¨ unster

slide-17
SLIDE 17

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

EM-Algorithm and TV-Regularization

min

u∈BV (Ω) u≥0

  • Σ

Ku − f log Ku + α |u|BV , α > 0 Combining nonlocal part (incl. K) and local regularization functional Convex, but not differentiable Optimality condition for u > 0 u

  • 1 − K ∗( f

Ku ) + α p

  • = 0,

p ∈ ∂|u|BV

13 Christoph Brune University of M¨ unster

slide-18
SLIDE 18

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

EM-Algorithm and TV-Regularization

Evaluate nonlocal term at uk and subgradient at uk+1 1 − K ∗( f Kuk ) + α pk+1 = 0, pk+1 ∈ ∂|uk+1|BV Two-step method uk+ 1

2 = ukK ∗( f

Kuk ) EM step uk+1 = uk+ 1

2 − α ukpk+1

TV step Alternating EM and TV

14 Christoph Brune University of M¨ unster

slide-19
SLIDE 19

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

EM-Algorithm and TV-Regularization

Evaluate nonlocal term at uk and subgradient at uk+1 better: uk+1 − ukK ∗( f Kuk ) + α ukpk+1 = 0, pk+1 ∈ ∂|uk+1|BV Two-step method uk+ 1

2 = ukK ∗( f

Kuk ) EM step uk+1 = uk+ 1

2 − α ukpk+1

TV step Alternating EM and TV

14 Christoph Brune University of M¨ unster

slide-20
SLIDE 20

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

EM-Algorithm and TV-Regularization

Iteration step uk+ 1

2 uk+1 can be realized by:

Weighted ROF-Problem uk+1 = arg min

u∈BV (Ω)

{ 1 2

(u − uk+ 1

2 )2

uk + α |u|BV } Optimality condition for uk+1 > 0 uk+1 = uk+ 1

2 − α ukpk+1,

pk+1 ∈ ∂|uk+1|BV Weighting data-fidelity term by solution uk of last TV problem

15 Christoph Brune University of M¨ unster

slide-21
SLIDE 21

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

EM-Algorithm and TV-Regularization

Iteration step uk+ 1

2 uk+1 can be realized by:

Weighted ROF-Problem uk+1 = arg min

u∈BV (Ω)

{ 1 2

(u − uk+ 1

2 )2

uk + α |u|BV } Optimality condition for uk+1 > 0 uk+1 = uk+ 1

2 − α ukpk+1,

pk+1 ∈ ∂|uk+1|BV Weighting data-fidelity term by solution uk of last TV problem Adaptive Regularization Numerical realization: → dual [Chambolle 04,05], [Aujol 08],... → primal-dual

15 Christoph Brune University of M¨ unster

slide-22
SLIDE 22

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

Example: Nanoscopy

EM KL-distance: 3.21 EM-TV, α = 0.04 KL-distance: 2.43

16 Christoph Brune University of M¨ unster

slide-23
SLIDE 23

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

Analysis EM-TV

Iterates well-defined in BV

existence (semicontinuity, coercivity) uniqueness (if K is injective) stability (regarding disturbances of K)

Nonnegativity: due to state constraints and fixed point scheme Positivity preservation for f > 0, u0 > 0: EM-step: uk+ 1

2 = ukK ∗( f

Kuk )

TV-step: Maximum principle for weighted ROF Descent of the objective functional with damping (yields uniform bound in BV hence stability)

[BSB, Preprint, 09]

17 Christoph Brune University of M¨ unster

slide-24
SLIDE 24

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

EM-TV and Forward-Backward Splitting

Decomposition in convex optimization Find u satisfying 0 ∈ 1 − K ∗ f Ku

  • A(u)

+ α ∂|u|BV

  • B(u)

Interpreting two-step method as forward-backward splitting uk+ 1

2 = ukK ∗

f Kuk

  • uk+1 = uk+ 1

2 − α ukpk+1 ,

pk+1 ∈ ∂|uk+1|BV

18 Christoph Brune University of M¨ unster

slide-25
SLIDE 25

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

EM-TV and Forward-Backward Splitting

Decomposition in convex optimization Find u satisfying 0 ∈ 1 − K ∗ f Ku

  • A(u)

+ α ∂|u|BV

  • B(u)

Interpreting two-step method as forward-backward splitting uk+ 1

2 − uk

uk + Auk = 0 , forward step for A uk+1 − uk+ 1

2

uk + Buk+1 = 0 , backward step for B

18 Christoph Brune University of M¨ unster

slide-26
SLIDE 26

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

EM-TV and Forward-Backward Splitting

Decomposition in convex optimization Find u satisfying 0 ∈ 1 − K ∗ f Ku

  • A(u)

+ α ∂|u|BV

  • B(u)

Interpreting two-step method as forward-backward splitting uk+ 1

2 − uk

uk + Auk = 0 , forward step for A uk+1 − uk+ 1

2

uk + Buk+1 = 0 , backward step for B = ⇒ uk+1 = (I + ukB)−1 (I − ukA) uk

[Douglas 56], [Lions 79], [Tseng 91], [Combettes 07,08], [Setzer 08]

18 Christoph Brune University of M¨ unster

slide-27
SLIDE 27

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

Convergence Results EM-TV

TV-step with damping, ωk ∈ [0, 1] uk+1 = arg min

u∈BV (Ω)

1 2

(u − (ωk uk+ 1

2 + (1 − ωk)uk))2

uk + ωk α|u|BV Forward-Backward Splitting = ⇒ uk+1 = (I + ωkukB)−1 (I − ωkukA) uk Convergence for general K if ǫ ≤ ωkuk ≤

4m 1−ǫ,

ǫ ∈ (0, 2m] Modulus m describes strength of convexity of KL(f , K·)

[BSB, Preprint, 09]

19 Christoph Brune University of M¨ unster

slide-28
SLIDE 28

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

Convergence Results EM-TV

TV-step with damping, ωk ∈ [0, 1] uk+1 = arg min

u∈BV (Ω)

1 2

(u − (ωk uk+ 1

2 + (1 − ωk)uk))2

uk + ωk α|u|BV Forward-Backward Splitting = ⇒ uk+1 = (I + ωkukB)−1 (I − ωkukA) uk Theorem

∃ ǫ > 0 0 < c1 ≤ ωk ≤

(uk+1 − uk)2 uk dλ

  • (1 − ǫ)

sup

v∈[uk ,uk+1]

1 2

  • Σ

f (Kuk+1 − Kuk)2 (Kv)2 dµ [BSB, Preprint, 09]

19 Christoph Brune University of M¨ unster

slide-29
SLIDE 29

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

Convergence Results EM-TV

TV-step with damping, ωk ∈ [0, 1] uk+1 = arg min

u∈BV (Ω)

1 2

(u − (ωk uk+ 1

2 + (1 − ωk)uk))2

uk + ωk α|u|BV Forward-Backward Splitting = ⇒ uk+1 = (I + ωkukB)−1 (I − ωkukA) uk Corollary K = I

∃ ǫ > 0 0 < c1 ≤ ωk+1 ≤ 2 (1 − ǫ) (min {(1 − ωk)f + ωkuk})2 max f max {(1 − ωk)f + ωkuk} [BSB, Preprint, 09]

19 Christoph Brune University of M¨ unster

slide-30
SLIDE 30

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

Convergence Results EM-TV

TV-step with damping, ωk ∈ [0, 1] uk+1 = arg min

u∈BV (Ω)

1 2

(u − (ωk uk+ 1

2 + (1 − ωk)uk))2

uk + ωk α|u|BV Forward-Backward Splitting = ⇒ uk+1 = (I + ωkukB)−1 (I − ωkukA) uk α = 10 ω ≡ 1 ω ≡ 0.2 ω ≡ 0.07

[BSB, Preprint, 09]

19 Christoph Brune University of M¨ unster

slide-31
SLIDE 31

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

Systematic Error of TV

Variational problem so far: min

u∈BV (Ω) u≥0

  • Σ

Ku − f log Ku + α |u|BV Experimental data: Nanoscopy

Protein Syntaxin (1024x1024) EM-TV Reconstruction

Systematic error of TV - contrast reduction

[ Meyer 01 ], [ Osher-Burger-Goldfarb-Xu-Yin 05 ]

20 Christoph Brune University of M¨ unster

slide-32
SLIDE 32

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

Contrast Enhancement by Iterative Regularization

Compute u1 via proposed EM-TV Bregman-EM-TV (Inverse Scale Space Method) pl ∈ ∂|ul|BV ul+1 = arg min

u∈BV (Ω) u≥0

{

  • Σ

(Ku − f log Ku) + α (|u|BV − |ul|BV − pl, u − ul

  • Dpl

|·|BV (u,ul) Bregman distance

) } Dpl

|·|BV (u, ul) ≥ 0

= 0 if pl ∈ ∂|u|BV particularly if u ≡ c ul ∀c ∈ R+ → edges at same position enhanced contrast

21 Christoph Brune University of M¨ unster

slide-33
SLIDE 33

Introduction EM-TV Methods Numerical Results EM / Richardson-Lucy Nested EM-TV and Analysis Contrast Enhancement by Bregman-EM-TV

Bregman-EM-TV Method

Realization via EM-TV splitting framework (for each l)          uk+ 1

2

= uk K ∗ f Kuk

  • (EM step)

uk+1 = arg min

u∈BV (Ω)

  • 1

2

(u − (uk+ 1

2 + αukpl)2

uk + α |u|BV

  • (TV step)

         Update: pl+1 = pl − 1 αul+1

  • 1 − K ∗
  • f

Kul+1

  • Dual Bregman-EM-TV

Contrast enhancement for general data fidelities Interpretation: residual ↔ dynamic background model Error estimates / convergence rates for dual ISS method [BSB - Primal and Dual Bregman Methods 09, UCLA CAM # 09-47]

22 Christoph Brune University of M¨ unster

slide-34
SLIDE 34

Introduction EM-TV Methods Numerical Results

Nanoscopy: Testobject

23 Christoph Brune University of M¨ unster

slide-35
SLIDE 35

Introduction EM-TV Methods Numerical Results

Nanoscopy: Testobject, EM-TV

Figure: KL-distance: 2.43

24 Christoph Brune University of M¨ unster

slide-36
SLIDE 36

Introduction EM-TV Methods Numerical Results

Nanoscopy: Testobject, Bregman-EM-TV

Figure: KL-distance: 1.43

25 Christoph Brune University of M¨ unster

slide-37
SLIDE 37

Introduction EM-TV Methods Numerical Results

Bregman-EM-TV Methods, Optimality

Figure: left: Functional values, right: KL(u, ˜ u)

26 Christoph Brune University of M¨ unster

slide-38
SLIDE 38

Introduction EM-TV Methods Numerical Results

Experimental Results 2D Nanoscopy

Figure: Protein Syntaxin in cell membrane, 53nm, MPI G¨

  • ttingen

27 Christoph Brune University of M¨ unster

slide-39
SLIDE 39

Introduction EM-TV Methods Numerical Results

Syntaxin, EM-TV Reconstruction

Figure: EM-TV Reconstruction

28 Christoph Brune University of M¨ unster

slide-40
SLIDE 40

Introduction EM-TV Methods Numerical Results

Syntaxin, Bregman-EM-TV Reconstruction

Figure: Bregman-EM-TV Reconstruction

29 Christoph Brune University of M¨ unster

slide-41
SLIDE 41

Introduction EM-TV Methods Numerical Results

3D Reconstruction Bregman-EM-TV

Colloidal crystal structure

(Harke et. al, MPI G¨

  • ttingen)

3D Bregman-EM-TV

30 Christoph Brune University of M¨ unster

slide-42
SLIDE 42

Introduction EM-TV Methods Numerical Results

Conclusions and Open Questions

Conclusions: Statistical model for deconvolution problems with Poisson noise EM-TV: accurate, fast and robust reconstruction algorithm Analytical results: well-posedness, positivity preservation, convergence Simultaneous contrast enhancement via Bregman Methods 2D and 3D results in optical nanoscopy Open Questions / Future Work: Appropriate stopping rules for Bregman-EM-TV methods Joint 4D reconstruction and motion estimation/compensation Nonlocal regularization techniques

31 Christoph Brune University of M¨ unster

slide-43
SLIDE 43

Introduction EM-TV Methods Numerical Results

Thank you for your attention! Questions? http://imaging.uni-muenster.de Acknowledgements:

32 Christoph Brune University of M¨ unster