[PPT] - Blind deconvolution of 3D data in wide field fluorescence microscopy PowerPoint Presentation

SLIDE 1

Blind deconvolution of 3D data in wide field fluorescence microscopy

Ferréol Soulez1,2 Loïc Denis2,3 Yves Tourneur1 Eric Thiébaut2

1Centre Commun de Quantimétrie

Lyon I, France

2Centre de Recherche Astrophysique de Lyon

Lyon I, France

3Laboratoire Hubert Curien

St Etienne, France

ISBI 2012

1 / 1

SLIDE 2

Wide Field Fluorescence Microscopy

Uniform illumination of the whole specimen, Imaging at the emission wavelenght, Moving the focal plane produces a 3D representation

f the specimen.

Coarse depth resolution Improving resolution Improving PSF (confocal,

multiphoton. . . ),

single molecule microscopy, Deconvolution. from Griffa et al. (2010).

2 / 1

SLIDE 3

Wide Field Fluorescence Microscopy

Uniform illumination of the whole specimen, Imaging at the emission wavelenght, Moving the focal plane produces a 3D representation

f the specimen.

Coarse depth resolution Improving resolution Improving PSF (confocal,

multiphoton. . . ),

single molecule microscopy, Deconvolution. from Griffa et al. (2010).

2 / 1

SLIDE 4

Wide Field Fluorescence Microscopy

Uniform illumination of the whole specimen, Imaging at the emission wavelenght, Moving the focal plane produces a 3D representation

f the specimen.

Coarse depth resolution Improving resolution Improving PSF (confocal,

multiphoton. . . ),

single molecule microscopy, Deconvolution. from Griffa et al. (2010).

2 / 1

SLIDE 5

Blind deconvolution

Blur modeled by a convolution : y = h ∗ x + n Deconvolution : Estimating the crisp image x of the specimen given the data y, the PSF h and the noise n statistics.

See [Agard & Sedat, 1983], [Sibarita, 2005] and [Sarder, 2005].

But : The PSF is not known — theoretical diffraction-limited PSF no flexibility — measured PSF with calibration beads complex, noisy — estimated directly from the blurred images Blind deconvolution Previous works by [Markam et al. 1999], [Hom et al. 2007] and [Kenig et al.

2010].

3 / 1

SLIDE 6

Blind deconvolution

Blur modeled by a convolution : y = h ∗ x + n Deconvolution : Estimating the crisp image x of the specimen given the data y, the PSF h and the noise n statistics.

See [Agard & Sedat, 1983], [Sibarita, 2005] and [Sarder, 2005].

But : The PSF is not known — theoretical diffraction-limited PSF no flexibility — measured PSF with calibration beads complex, noisy — estimated directly from the blurred images Blind deconvolution Previous works by [Markam et al. 1999], [Hom et al. 2007] and [Kenig et al.

2010].

3 / 1

SLIDE 7

Blind deconvolution

Blur modeled by a convolution : y = h ∗ x + n Deconvolution : Estimating the crisp image x of the specimen given the data y, the PSF h and the noise n statistics.

See [Agard & Sedat, 1983], [Sibarita, 2005] and [Sarder, 2005].

But : The PSF is not known — theoretical diffraction-limited PSF no flexibility — measured PSF with calibration beads complex, noisy — estimated directly from the blurred images Blind deconvolution Previous works by [Markam et al. 1999], [Hom et al. 2007] and [Kenig et al.

2010].

3 / 1

SLIDE 8

Blind deconvolution

Blur modeled by a convolution : y = h ∗ x + n Deconvolution : Estimating the crisp image x of the specimen given the data y, the PSF h and the noise n statistics.

See [Agard & Sedat, 1983], [Sibarita, 2005] and [Sarder, 2005].

But : The PSF is not known — theoretical diffraction-limited PSF no flexibility — measured PSF with calibration beads complex, noisy — estimated directly from the blurred images Blind deconvolution Previous works by [Markam et al. 1999], [Hom et al. 2007] and [Kenig et al.

2010].

3 / 1

SLIDE 9

Blind deconvolution

Blur modeled by a convolution : y = h ∗ x + n Deconvolution : Estimating the crisp image x of the specimen given the data y, the PSF h and the noise n statistics.

See [Agard & Sedat, 1983], [Sibarita, 2005] and [Sarder, 2005].

But : The PSF is not known — theoretical diffraction-limited PSF no flexibility — measured PSF with calibration beads complex, noisy — estimated directly from the blurred images Blind deconvolution Previous works by [Markam et al. 1999], [Hom et al. 2007] and [Kenig et al.

2010].

3 / 1

SLIDE 10

Blind deconvolution

Blur modeled by a convolution : y = h ∗ x + n Deconvolution : Estimating the crisp image x of the specimen given the data y, the PSF h and the noise n statistics.

See [Agard & Sedat, 1983], [Sibarita, 2005] and [Sarder, 2005].

But : The PSF is not known — theoretical diffraction-limited PSF no flexibility — measured PSF with calibration beads complex, noisy — estimated directly from the blurred images Blind deconvolution Previous works by [Markam et al. 1999], [Hom et al. 2007] and [Kenig et al.

2010].

3 / 1

SLIDE 11

Blind deconvolution

Blur modeled by a convolution : y = h ∗ x + n Deconvolution : Estimating the crisp image x of the specimen given the data y, the PSF h and the noise n statistics.

See [Agard & Sedat, 1983], [Sibarita, 2005] and [Sarder, 2005].

But : The PSF is not known — theoretical diffraction-limited PSF no flexibility — measured PSF with calibration beads complex, noisy — estimated directly from the blurred images Blind deconvolution Previous works by [Markam et al. 1999], [Hom et al. 2007] and [Kenig et al.

2010].

3 / 1

SLIDE 12

Blind deconvolution

Blur modeled by a convolution : y = h ∗ x + n Deconvolution : Estimating the crisp image x of the specimen given the data y, the PSF h and the noise n statistics.

See [Agard & Sedat, 1983], [Sibarita, 2005] and [Sarder, 2005].

But : The PSF is not known — theoretical diffraction-limited PSF no flexibility — measured PSF with calibration beads complex, noisy — estimated directly from the blurred images Blind deconvolution Previous works by [Markam et al. 1999], [Hom et al. 2007] and [Kenig et al.

2010].

3 / 1

SLIDE 13

Maximum a posteriori blind deconvolution

Estimating the most probable couple Object/PSF {x+, h+} according to the data and some a priori knowledge. Done by the minimisation of a cost function J(x, h): J(x, h) = Jdata(x, h)

likelihood

+µ Jprior(x)

bject priors

, PSF priors enforced by its parametrization. Object a priori globally smooth with few sharp edges : Hyperbolic approximation of 3D total variation : Jprior(x) =

k
∇xk2

2 + ǫ2 .

4 / 1

SLIDE 14

Maximum a posteriori blind deconvolution

Estimating the most probable couple Object/PSF {x+, h+} according to the data and some a priori knowledge. Done by the minimisation of a cost function J(x, h): J(x, h) = Jdata(x, h)

likelihood

+µ Jprior(x)

bject priors

, PSF priors enforced by its parametrization. Object a priori globally smooth with few sharp edges : Hyperbolic approximation of 3D total variation : Jprior(x) =

k
∇xk2

2 + ǫ2 .

4 / 1

SLIDE 15

Maximum a posteriori blind deconvolution

Estimating the most probable couple Object/PSF {x+, h+} according to the data and some a priori knowledge. Done by the minimisation of a cost function J(x, h): J(x, h) = Jdata(x, h)

likelihood

+µ Jprior(x)

bject priors

, PSF priors enforced by its parametrization. Object a priori globally smooth with few sharp edges : Hyperbolic approximation of 3D total variation : Jprior(x) =

k
∇xk2

2 + ǫ2 .

4 / 1

SLIDE 16

Likelihood

Gaussian noise : Jdata(x) = 1 2(y − H · x)T · C−1

noise · (y − H · x)

Uncorrelated non-stationnary Gaussian noise : Jdata(x) =

k=Pixels
λ

1 σk,λ (H · x)k − yk,λ 2 Missing pixels k −→ σk,λ = ∞. Poisson Noise ≈ non-stationnary Gaussian noise σk,λ = γ(H · x)k,λ + σ2

CCD≈ γ max(yk,λ, 0) + σ2 CCD

where γ is a quantization factor and σ2

CCD account for Gaussian

additive noise (e.g. readout noise).

5 / 1

SLIDE 17

Likelihood

Gaussian noise : Jdata(x) = 1 2(y − H · x)T · C−1

noise · (y − H · x)

Uncorrelated non-stationnary Gaussian noise : Jdata(x) =

k=Pixels
λ

1 σk,λ (H · x)k − yk,λ 2 Missing pixels k −→ σk,λ = ∞. Poisson Noise ≈ non-stationnary Gaussian noise σk,λ = γ(H · x)k,λ + σ2

CCD≈ γ max(yk,λ, 0) + σ2 CCD

where γ is a quantization factor and σ2

CCD account for Gaussian

additive noise (e.g. readout noise).

5 / 1

SLIDE 18

Likelihood

Gaussian noise : Jdata(x) = 1 2(y − H · x)T · C−1

noise · (y − H · x)

Uncorrelated non-stationnary Gaussian noise : Jdata(x) =

k=Pixels
λ

1 σk,λ (H · x)k − yk,λ 2 Missing pixels k −→ σk,λ = ∞. Poisson Noise ≈ non-stationnary Gaussian noise σk,λ = γ(H · x)k,λ + σ2

CCD≈ γ max(yk,λ, 0) + σ2 CCD

where γ is a quantization factor and σ2

CCD account for Gaussian

additive noise (e.g. readout noise).

5 / 1

SLIDE 19

PSF parametrization

PSF h defined as a function of the pupil function as Markham (1999) h(rj, z) =

k Fj,k ak(z)
2

, with rj lateral position of pixel j, F discrete Fourier transform and ak(z) pupil function at frequel k and depth z. ak(z) = ρk exp(i 2 π (φk + z ψk)) , ρk =

n βnZn

k ,

φk =

n αnZn

k ,

ψk =

(ni/λ)2 − κk2

where Zn

k the n-th Zernike polynomial [Hanser 2004] and ni the refractive

index of immersion medium. PSF parametrized by {ni, α, β}

6 / 1

SLIDE 20

PSF parametrization

PSF h defined as a function of the pupil function as Markham (1999) h(rj, z) =

k Fj,k ak(z)
2

, with rj lateral position of pixel j, F discrete Fourier transform and ak(z) pupil function at frequel k and depth z. ak(z) = ρk exp(i 2 π (φk + z ψk)) , ρk =

n βnZn

k ,

φk =

n αnZn

k ,

ψk =

(ni/λ)2 − κk2

where Zn

k the n-th Zernike polynomial [Hanser 2004] and ni the refractive

index of immersion medium. PSF parametrized by {ni, α, β}

6 / 1

SLIDE 21

PSF parametrization

PSF h defined as a function of the pupil function as Markham (1999) h(rj, z) =

k Fj,k ak(z)
2

, with rj lateral position of pixel j, F discrete Fourier transform and ak(z) pupil function at frequel k and depth z. ak(z) = ρk exp(i 2 π (φk + z ψk)) , ρk =

n βnZn

k ,

φk =

n αnZn

k ,

ψk =

(ni/λ)2 − κk2

where Zn

k the n-th Zernike polynomial [Hanser 2004] and ni the refractive

index of immersion medium. PSF parametrized by {ni, α, β}

6 / 1

SLIDE 22

PSF parametrization

PSF h defined as a function of the pupil function as Markham (1999) h(rj, z) =

k Fj,k ak(z)
2

, with rj lateral position of pixel j, F discrete Fourier transform and ak(z) pupil function at frequel k and depth z. ak(z) = ρk exp(i 2 π (φk + z ψk)) , ρk =

n βnZn

k ,

φk =

n αnZn

k ,

ψk =

(ni/λ)2 − κk2

where Zn

k the n-th Zernike polynomial [Hanser 2004] and ni the refractive

index of immersion medium. PSF parametrized by {ni, α, β}

6 / 1

SLIDE 23

PSF parametrization

PSF h defined as a function of the pupil function as Markham (1999) h(rj, z) =

k Fj,k ak(z)
2

, with rj lateral position of pixel j, F discrete Fourier transform and ak(z) pupil function at frequel k and depth z. ak(z) = ρk exp(i 2 π (φk + z ψk)) , ρk =

n βnZn

k ,

φk =

n αnZn

k ,

ψk =

(ni/λ)2 − κk2

where Zn

k the n-th Zernike polynomial [Hanser 2004] and ni the refractive

index of immersion medium. PSF parametrized by {ni, α, β}

6 / 1

SLIDE 24

PSF parametrization

PSF h defined as a function of the pupil function as Markham (1999) h(rj, z) =

k Fj,k ak(z)
2

, with rj lateral position of pixel j, F discrete Fourier transform and ak(z) pupil function at frequel k and depth z. ak(z) = ρk exp(i 2 π (φk + z ψk)) , ρk =

n βnZn

k ,

φk =

n αnZn

k ,

ψk =

(ni/λ)2 − κk2

where Zn

k the n-th Zernike polynomial [Hanser 2004] and ni the refractive

index of immersion medium. PSF parametrized by {ni, α, β}

6 / 1

SLIDE 25

PSF parametrization

Preventing some degeneracies : Centering PSF : removing phase tip-tilt α1 = α2 = 0. normalizing PSF

h(k)dk = 1 : constraining

n β2 n = 1.

Benefits of such parametrization :

ptically derived model,

require only the knowledge of the wavelength λ, the numerical aperture NA, few parameters (several tenth), no additional priors (regularization), ensure PSF positivity, taking only radial Zernike polynomials ensure axially symmetric PSF .

7 / 1

SLIDE 26

PSF parametrization

Preventing some degeneracies : Centering PSF : removing phase tip-tilt α1 = α2 = 0. normalizing PSF

h(k)dk = 1 : constraining

n β2 n = 1.

Benefits of such parametrization :

ptically derived model,

require only the knowledge of the wavelength λ, the numerical aperture NA, few parameters (several tenth), no additional priors (regularization), ensure PSF positivity, taking only radial Zernike polynomials ensure axially symmetric PSF .

7 / 1

SLIDE 27

Algorithm summary

Solution given by : {x+, n+

i , α+, β+} = arg min x,ni,α,β

Jdata(x, h(ni, α, β; y)) + µ Jprior(x)
Non-convex and badly conditioned problem.

Alternating minimization : Begin with aberrations free PSF h(0) (α = β = 0), set n = 1:

1

x(n) = arg min

x

Jdata(x, h(n−1); y) + µ Jprior(x)
2

n(n)

i

= arg min

ni

Jdata(x(n), h(ni, α(n−1), β(n−1)); y)

3

α(n) = arg min

α

Jdata(x(n), h(n(n)

i , α, β(n−1)); y)

4

β(n) = arg min

β

Jdata(x(n), h(n(n)

i , α(n), β; y) under constraint k β2 k = 1.

5

n = n + 1, go to step 1 until a certain convergence.

8 / 1

SLIDE 28

Algorithm summary

Solution given by : {x+, n+

i , α+, β+} = arg min x,ni,α,β

Jdata(x, h(ni, α, β; y)) + µ Jprior(x)
Non-convex and badly conditioned problem.

Alternating minimization : Begin with aberrations free PSF h(0) (α = β = 0), set n = 1:

1

x(n) = arg min

x

Jdata(x, h(n−1); y) + µ Jprior(x)
2

n(n)

i

= arg min

ni

Jdata(x(n), h(ni, α(n−1), β(n−1)); y)

3

α(n) = arg min

α

Jdata(x(n), h(n(n)

i , α, β(n−1)); y)

4

β(n) = arg min

β

Jdata(x(n), h(n(n)

i , α(n), β; y) under constraint k β2 k = 1.

5

n = n + 1, go to step 1 until a certain convergence.

8 / 1

SLIDE 29

Algorithm summary

Solution given by : {x+, n+

i , α+, β+} = arg min x,ni,α,β

Jdata(x, h(ni, α, β; y)) + µ Jprior(x)
Non-convex and badly conditioned problem.

Alternating minimization : Begin with aberrations free PSF h(0) (α = β = 0), set n = 1:

1

x(n) = arg min

x

Jdata(x, h(n−1); y) + µ Jprior(x)
2

n(n)

i

= arg min

ni

Jdata(x(n), h(ni, α(n−1), β(n−1)); y)

3

α(n) = arg min

α

Jdata(x(n), h(n(n)

i , α, β(n−1)); y)

4

β(n) = arg min

β

Jdata(x(n), h(n(n)

i , α(n), β; y) under constraint k β2 k = 1.

5

n = n + 1, go to step 1 until a certain convergence.

8 / 1

SLIDE 30

Algorithm summary

Solution given by : {x+, n+

i , α+, β+} = arg min x,ni,α,β

Jdata(x, h(ni, α, β; y)) + µ Jprior(x)
Non-convex and badly conditioned problem.

Alternating minimization : Begin with aberrations free PSF h(0) (α = β = 0), set n = 1:

1

x(n) = arg min

x

Jdata(x, h(n−1); y) + µ Jprior(x)
2

n(n)

i

= arg min

ni

Jdata(x(n), h(ni, α(n−1), β(n−1)); y)

3

α(n) = arg min

α

Jdata(x(n), h(n(n)

i , α, β(n−1)); y)

4

β(n) = arg min

β

Jdata(x(n), h(n(n)

i , α(n), β; y) under constraint k β2 k = 1.

5

n = n + 1, go to step 1 until a certain convergence.

8 / 1

SLIDE 31

Algorithm summary

Solution given by : {x+, n+

i , α+, β+} = arg min x,ni,α,β

Jdata(x, h(ni, α, β; y)) + µ Jprior(x)
Non-convex and badly conditioned problem.

Alternating minimization : Begin with aberrations free PSF h(0) (α = β = 0), set n = 1:

1

x(n) = arg min

x

Jdata(x, h(n−1); y) + µ Jprior(x)
2

n(n)

i

= arg min

ni

Jdata(x(n), h(ni, α(n−1), β(n−1)); y)

3

α(n) = arg min

α

Jdata(x(n), h(n(n)

i , α, β(n−1)); y)

4

β(n) = arg min

β

Jdata(x(n), h(n(n)

i , α(n), β; y) under constraint k β2 k = 1.

5

n = n + 1, go to step 1 until a certain convergence.

8 / 1

SLIDE 32

Algorithm summary

Solution given by : {x+, n+

i , α+, β+} = arg min x,ni,α,β

Jdata(x, h(ni, α, β; y)) + µ Jprior(x)
Non-convex and badly conditioned problem.

Alternating minimization : Begin with aberrations free PSF h(0) (α = β = 0), set n = 1:

1

x(n) = arg min

x

Jdata(x, h(n−1); y) + µ Jprior(x)
2

n(n)

i

= arg min

ni

Jdata(x(n), h(ni, α(n−1), β(n−1)); y)

3

α(n) = arg min

α

Jdata(x(n), h(n(n)

i , α, β(n−1)); y)

4

β(n) = arg min

β

Jdata(x(n), h(n(n)

i , α(n), β; y) under constraint k β2 k = 1.

5

n = n + 1, go to step 1 until a certain convergence.

8 / 1

SLIDE 33

Algorithm summary

Solution given by : {x+, n+

i , α+, β+} = arg min x,ni,α,β

Jdata(x, h(ni, α, β; y)) + µ Jprior(x)
Non-convex and badly conditioned problem.

Alternating minimization : Begin with aberrations free PSF h(0) (α = β = 0), set n = 1:

1

x(n) = arg min

x

Jdata(x, h(n−1); y) + µ Jprior(x)
2

n(n)

i

= arg min

ni

Jdata(x(n), h(ni, α(n−1), β(n−1)); y)

3

α(n) = arg min

α

Jdata(x(n), h(n(n)

i , α, β(n−1)); y)

4

β(n) = arg min

β

Jdata(x(n), h(n(n)

i , α(n), β; y) under constraint k β2 k = 1.

5

n = n + 1, go to step 1 until a certain convergence.

8 / 1

SLIDE 34

Blind deconvolution on simulations

Ground truth Data proposed method XY section XZ section

Simulation with depth aberrations from Kenig, Kam & Feuer, TPAMI, (2010)

9 / 1

SLIDE 35

Blind deconvolution on simulations

Ground truth Data proposed method XY section XZ section

Simulation with depth aberrations from Kenig, Kam & Feuer, TPAMI, (2010)

9 / 1

SLIDE 36

Blind deconvolution on simulations

Ground truth Data Kenig et al. Proposed method XY section XZ section RMSE 36.16 24.88 11.27

Simulation with depth aberrations from Kenig, Kam & Feuer, TPAMI, (2010)

10 / 1

SLIDE 37

Experimental results: Calibration bead

XY section XZ section

5 10 15 5 10 15

5 10 15 5 10 15 20

1e+03 2e+03 3e+03 4e+03

— Bead diameter: 2.5µm, NA = 1.4 — 2563 pixels 64.5 × 64.5 × 160nm3 from A. Griffa, N. Garin & D. Sage, G.I.T. Imaging & Microscopy, 2010.

11 / 1

SLIDE 38

Non blind deconvolution with theoretical PSF

XY section XZ section

5 10 15 5 10 15

5 10 15 5 10 15 20

1.64e+04 3.28e+04 4.92e+04 6.55e+04

— Bead diameter: 2.5µm, NA = 1.4, λ = 512nm — 2563 pixels 64.5 × 64.5 × 160nm3 from A. Griffa, N. Garin & D. Sage, G.I.T. Imaging & Microscopy, 2010.

12 / 1

SLIDE 39

Calibration bead : blind deconvolution

XY section XZ section

5 10 15 5 10 15

5 10 15 5 10 15 20

— Bead diameter: 2.5µm, NA = 1.4, λ = 512nm — 2563 voxels 64.5 × 64.5 × 160nm3 from A. Griffa, N. Garin & D. Sage, G.I.T. Imaging & Microscopy, 2010.

13 / 1

SLIDE 40

Calibration bead

5 10 0.0 0.5 1.0

Data Deconvolution Blind deconvolution

3D Radial profile of the bead

data Hyugens AutoDeblur Deconvol. proposed method parameters Lab non-blind blind transversal FWHM 2.87 2.71 2.71 2.66 2.74 2.78 axial FWHM (in µm) 4.76 4.00 4.64 4.16 3.05 2.98 Relative contrast 18 % 53 % 78 % 68 % 84 % 88 %

Performance of 3 deconvolution methods as reported by Griffa (2010) compared to the proposed method. Hyugens and AutoDeblur are commercial softwares and Deconvolution Lab is an imageJ plugin implementing (Vonesch, 2008).

14 / 1

SLIDE 41

Calibration bead: PSF

Theoretical PSF

−1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 −3 −2 −1 1 2 3 −4 −2 2 4

Estimated PSF

−1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 −3 −2 −1 1 2 3 −4 −2 2 4 −3 −2 −1 1 2 3 −4 −2 2 4 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 0.0007

XZ section XZ section YZ section

15 / 1

SLIDE 42

Experimental result: C. Elegans

C. Elegans embryo

×63, 1.4 NA oil

bjective,

DAPI + FITC + CY3, 672 × 712 × 104 voxels, voxels size 64.5 × 64.5 × 200 nm3 from A. Griffa, N. Garin & D. Sage, G.I.T. Imaging & Microscopy, 2010.

16 / 1

SLIDE 43

Experimental result: C. Elegans

17 / 1

SLIDE 44

Conclusion

An effective blind deconvolution method increase both lateral and axial resolution,

ptically motivated PSF model,

few needed parameters (NA and wavelength), But still one hyper-parameter to tune. Works in progress extending to confocal and two photons microscopy, using [Denis et al 2011] for depth variant blind deconvolution.

18 / 1

SLIDE 45

Conclusion

An effective blind deconvolution method increase both lateral and axial resolution,

ptically motivated PSF model,

few needed parameters (NA and wavelength), But still one hyper-parameter to tune. Works in progress extending to confocal and two photons microscopy, using [Denis et al 2011] for depth variant blind deconvolution.

18 / 1