[PPT] - Combined Image Reconstruction for Combined PET-MR Imaging Matthias PowerPoint Presentation

SLIDE 1

Combined Image Reconstruction for Combined PET-MR Imaging

Matthias J. Ehrhardt

University of Cambridge, UK with: Arridge, Atkinson, Barnes, Duncan, Hutton, Liljeroth, Markiewicz Ourselin, Pizarro, Thielemans (UCL, London) Kolehmainen (Kuopio, Finland)

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 2

Positron Emission Tomography and Magnetic Resonance Imaging

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 3

Positron Emission Tomography (PET)

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 4

Positron Emission Tomography (PET)

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 5

Positron Emission Tomography (PET)

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 6

Positron Emission Tomography (PET)

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 7

Magnetic Resonance Imaging (MRI)

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 8

Magnetic Resonance Imaging (MRI)

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 9

Magnetic Resonance Imaging (MRI)

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 10

Magnetic Resonance Imaging (MRI)

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 11

Magnetic Resonance Imaging (MRI)

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 12

Combined PET-MR Imaging

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 13

Combined PET-MR Imaging

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 14

Combined PET-MR Imaging

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 15

Part I: Utilizing Resolution of MRI Part II: Joint PET-MRI Reconstruction

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 16

Part I: Utilizing Resolution of MRI Part II: Joint PET-MRI Reconstruction

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 17

PET Reconstruction

PET data MLEM TV MRI

?

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 18

MAP reconstruction and Total Variation

MAP reconstruction

u∗ ∈ argmin

u

L(Au + r, b) + αR(u)
◮ total variation Rudin, Osher, Fatemi 1992

R(u) = TV(u) =

Ω

|∇u| R(u) = TVβ(u) =

Ω
β2 + |∇u|21/2

edge-preserved reconstruction

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 19

MAP reconstruction and Anatomical Information

MAP reconstruction with Anatomical Information

u∗ ∈ argmin

u

L(Pu + r, b) + αR(u|v)
We want

1) Convexity: R(u|v) should be convex in u 2) No Segmentation: should not need a segmentation of v 3∗) Total Variation: R(u|v = const) = TV(u) PET with MRI/CT: Leahy and Yan 1991, Baete et al 2004, Pedemonte et al

2011, Bowsher et al 2004, Kazantsev et al 2014, Nuyts 2007, Somayayula et al 2005 2011, Tang and Rahmim 2009 2015 (Mutual information/ Entropy), Jiao et al 2015

EIT with CT: Kaipio et al 1999

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 20

Parallel Level Set Prior

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 21

Parallel Level Set Prior

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 22

Parallel Level Set Prior

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 23

Parallel Level Set Prior

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 24

Parallel Level Set Prior

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 25

Parallel Level Set Prior

∇u, ∇v = cos(θ)|∇u||∇v|

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 26

Parallel Level Set Prior

∇u, ∇v = cos(θ)|∇u||∇v|

Measure Similar Structures

S(u) :=

|∇u|2 − ∇u, ∇v/|∇v|21/2

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 27

Parallel Level Set Prior

∇u, ∇v = cos(θ)|∇u||∇v|

Measure Similar Structures

S(u) :=

|∇u|2 − ∇u, ξ21/2

◮ ξ := ∇v/|∇v|η,

|∇v|η :=

|∇v|2 + η2,

η > 0

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 28

Parallel Level Set Prior

∇u, ∇v = cos(θ)|∇u||∇v|

Measure Similar Structures

S(u) :=

|∇u|2 − ∇u, ξ21/2

◮ ξ := ∇v/|∇v|η,

|∇v|η :=

|∇v|2 + η2,

η > 0

◮ 0 ≤ S(u) ≤ |∇u| ◮ S(u) = 0 ⇔ u ∼ v (∇u ∇v) m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 29

Asymmetric Parallel Level Sets

S(u) :=

|∇u|2 − ∇u, ξ21/2

Asymmetric Parallel Level Sets Prior

P(u|v) :=

Ω
β2 + |∇u|2 − ∇u, ξ21/2

, β > 0 This is convex, does not need a segmentation and reduces to total variation.

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 30

Numerical Results

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 31

Other Methods for Anatomical Information

TVJ (u|v) :=

Ω
β2 + |∇u|2 + γ|∇v|21/2

, γ > 0

Sapiro and Ringach IEEE TIP 1996; Haber and Holtzman-Gazit Surveys in Geophysics 2013; Ehrhardt et al Inv Probl 2015, Lu et al Phys Med Bio 2015

B(u|v) := 1 2

i
j∈N(i)

ωi,j(v)(ui − uj)2, k ∈ N

Bowsher et al IEEE NSS-MIC 2004

D(u|v) :=

Ω
β2 + |∇u|21/2

− ∇u, ξ

Kazantsev et al Sensing and Imaging 2014

K(u|v) := 1 2

Ω

|∇u|2 − ∇u, ξ2

Kaipio et al Inv Prob 1999 m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 32

Summary of Methods

TVJ B D K P reduces to total variation ✓ ✗ ✓ ✗ ✓ edge location dependent ✓ ✓ ✓ ✓ ✓ edge orientation dependent ✗ ✗ ✓ ✓ ✓ allows negative edge correlation

✗

✓ ✓

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 33

Data

a) b) c)

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 34

Software Phantom: Quantitative Results

relative ℓ2-error [%] whole phantom 30 40 50 SSIM [%] whole phantom 70 80 relative ℓ2-error [%] grey matter low high 20 30 40 regularization relative ℓ2-error [%] right hot lesion low high 20 30 40 TV TV

J

B D K P† MLEM

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 35

Software Phantom: Normal Recon v Anatomical Prior

PET gr. truth MRI MLEM TV P

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 36

Software Phantom: Compare Anatomical Priors

PET gr. truth MRI TVJ B D K P

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 37

Software Phantom: Close-Ups

TV

J

B D K P† MRI side info PET ground truth

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 38

Software Phantom: Bias vs SD

mean standard deviation whole phantom 0.01 0.02 0.03 0.01 0.02 0.03 grey matter 0.03 0.06 0.09 0.03 0.06 0.09 reg → 0 reg → ∞ white matter 0.03 0.06 0.09 0.03 0.06 0.09 mean absolute bias lesions 0.13 0.26 0.39 0.13 0.26 0.39 cold left hot right hot TV TV

J

B D K P† MLEM

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 39

Software Phantom: Bias vs SD

≥ 0.5 ≤ -0.5 ≥ 0.25 MLEM TV TV

J

B D K P† bias standard deviation

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 40

Software Phantom: Bias vs SD

MLEM TV bias standard deviation ≥ 0.5 ≤ -0.5 ≥ 0.25 K P† m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 41

Hardware Phantom: Compare Anatomical Priors

MRI TVJ B D K P

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 42

Hardware Phantom: Close-Ups

MLEM TV TV

J

B D K P† MRI side info

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 43

Patient Data: Normal Recon v Anatomical Prior

MRI MLEM TV P

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 44

Patient Data: Compare Anatomical Priors

MRI TVJ B D K P

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 45

Patient Data: Close-Ups

MLEM TV TV

J

B D K P† MRI side info

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al 2016 (under review)

SLIDE 46

Conclusions of Part I

◮ new prior that can incorporate anatomical structure

◮ convex, no segmentation and reduces to total variation ◮ based on directions, not only magnitude ◮ handles arbitrary intensities, no need for positive correlation ◮ better in quality measures (ℓ2-error, SSIM, bias-vs-SD

trade-off)

◮ reduces bias of total variation (similar to Bregman iterations)

MLEM P

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 47

Part I: Utilizing Resolution of MRI Part II: Joint PET-MRI Reconstruction

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 48

Data Acquisition in MRI

◮ sequential sampling of Fourier

coefficients

◮ less data

⇒ shorter acquisition time ⇒ motion, patient comfort, money

◮ higher spatial resolution m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 49

Joint Reconstruction

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 50

Joint Reconstruction

?

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al Inverse Problems 2015

SLIDE 51

Joint Reconstruction

?

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al Inverse Problems 2015

SLIDE 52

Joint Reconstruction

? ?

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al Inverse Problems 2015

SLIDE 53

Joint Reconstruction

? ?

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al Inverse Problems 2015

SLIDE 54

Problem Set Up

◮ Reconstruct jointly PET and MRI ◮ Two modalities with different characteristics

MRI:

◮ Undersampled Fourier data with Gaussian noise ◮ Forward operator is not injective but pseudo inverse is

well-conditioned

PET:

◮ Blurry Radon data with Poisson noise ◮ Forward operator compact, inverse is ill-conditioned

◮ two problems coupled by underlying anatomy m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al Inverse Problems 2015

SLIDE 55

Framework for Joint Reconstruction

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al Inverse Problems 2015

SLIDE 56

Framework for Joint Reconstruction

p(u, v|f , g) u v f g

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al Inverse Problems 2015

SLIDE 57

Framework for Joint Reconstruction

p(u, v|f , g) ∝ p(f , g|u, v)p(u, v) u v f g

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al Inverse Problems 2015

SLIDE 58

Framework for Joint Reconstruction

p(u, v|f , g) ∝ p(f , g|u, v)p(u, v) = p(f |u, v)p(g|u, v)p(u, v) u v f g

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al Inverse Problems 2015

SLIDE 59

Framework for Joint Reconstruction

p(u, v|f , g) ∝ p(f , g|u, v)p(u, v) = p(f |u, v)p(g|u, v)p(u, v) = p(f |u)p(g|v)p(u, v) u v f g

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al Inverse Problems 2015

SLIDE 60

Framework for Joint Reconstruction

p(u, v|f , g) ∝ p(f , g|u, v)p(u, v) = p(f |u, v)p(g|u, v)p(u, v) = p(f |u)p(g|v)p(u, v) u v f g minimize

u,v

− log p(f |u) − log p(g|v) − log p(u, v)
m.j.ehrhardt@damtp.cam.ac.uk

Ehrhardt et al Inverse Problems 2015

SLIDE 61

Framework for Joint Reconstruction

p(u, v|f , g) ∝ p(f , g|u, v)p(u, v) = p(f |u, v)p(g|u, v)p(u, v) = p(f |u)p(g|v)p(u, v) u v f g minimize

u,v

− log p(f |u) − log p(g|v) − log p(u, v)
∝ KL(Au + b; f ) +

1 2σ2 Bv − g2 − log p(u, v) KL(x; y) :=

j xj − yj + yj log(yj/xj) m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al Inverse Problems 2015

SLIDE 62

Parallel Level Sets

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 63

Joint Parallel Level Set Prior

∇u, ∇v = cos(θ)|∇u||∇v|

Measure Similar Structures

S(u) =

|∇u|2 − ∇u, ∇v/|∇v|21/2

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt and Arridge IEEE TIP 2014, Ehrhardt 2015

SLIDE 64

Joint Parallel Level Set Prior

∇u, ∇v = cos(θ)|∇u||∇v|

Measure Similar Structures

S(u, v) = |∇u| |∇v| − | ∇u, ∇v |

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt and Arridge IEEE TIP 2014, Ehrhardt 2015

SLIDE 65

Joint Parallel Level Set Prior

∇u, ∇v = cos(θ)|∇u||∇v|

Measure Similar Structures

S(u, v) = |∇u| |∇v| − | ∇u, ∇v |

◮ S(u, v) ≥ 0 ◮ S(u, v) = 0 ⇔ u ∼ v (∇u ∇v) m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt and Arridge IEEE TIP 2014, Ehrhardt 2015

SLIDE 66

Joint Parallel Level Set Prior

∇u, ∇v = cos(θ)|∇u||∇v|

Measure Similar Structures

S(u, v) =

Ω

|∇u| |∇v| − | ∇u, ∇v |

◮ S(u, v) ≥ 0 ◮ S(u, v) = 0 ⇔ u ∼ v (∇u ∇v almost everywhere) m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt and Arridge IEEE TIP 2014, Ehrhardt 2015

SLIDE 67

Joint Parallel Level Set Prior

Recall, S(u, v) =

|∇u| |∇v| − |

∇u, ∇v |.

Structure is Intensity Invariant

Let f ∈ C 1(R, R) (with f injective). Then, u ∼ v ⇒ (⇔) u ∼ v ◦ f Proof: At almost every location x, there is ∇(v ◦ f )(x) = f ′(v(x))∇v(x) = f ′(v(x))λ(x)∇u(x) = ˜ λ(x)∇u(x).

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt and Arridge IEEE TIP 2014, Ehrhardt 2015

SLIDE 68

Joint Parallel Level Set Prior

Recall, S(u, v) =

|∇u||∇v| − |

∇u, ∇v |.

Asymptotics

◮ For |∇v| ≈ 0, there is

S(u, v) ≈ 0.

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt and Arridge IEEE TIP 2014, Ehrhardt 2015

SLIDE 69

Joint Parallel Level Set Prior

Recall, S(u, v) =

|∇u||∇v| − |

∇u, ∇v |.

Asymptotics

◮ For |∇v| ≈ 0, there is

S(u, v) ≈ 0.

Parallel Level Sets Prior

Sβ(u, v) =

|∇u|β|∇v|β − |

∇u, ∇v |β2 with “smoothed” norm |x|β =

β2 + |x|2.

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt and Arridge IEEE TIP 2014, Ehrhardt 2015

SLIDE 70

Joint Parallel Level Set Prior

Recall, Sβ(u, v) =

|∇u|β|∇v|β − |

∇u, ∇v |β2.

Asymptotics

◮ For |∇v| ≈ 0, there is

S(u, v) ≈ 0.

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt and Arridge IEEE TIP 2014, Ehrhardt 2015

SLIDE 71

Joint Parallel Level Set Prior

Recall, Sβ(u, v) =

|∇u|β|∇v|β − |

∇u, ∇v |β2.

Asymptotics

◮ For |∇v| ≪ β, there is |∇v|β ≈ β, hence

Sβ(u, v) ≈

β|∇u|β + const = β TVβ(u) + const.

◮ For |∇u|, |∇v| ≫ β, there is |x|β ≈ |x|, hence

Sβ(u, v) ≈ S(u, v).

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt and Arridge IEEE TIP 2014, Ehrhardt 2015

SLIDE 72

Parallel Level Set Prior

Generalization

ϕ, ψ : [0, ∞) → [0, ∞), ϕ(0) = 0, both monotonically increasing Sϕ,ψ(u, v) =

ϕ
ψ
|∇u|β |∇v|β
− ψ
|

∇u, ∇v |β2

Special cases

◮ linear parallel level sets: ϕ(x), ψ(x) = x ◮ quadratic parallel level sets (Nambu functional):

ϕ(x) = √x, ψ(x) = x2

◮ cross-gradients: β = 0, ϕ(x) = x, ψ(x) = x2 m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt and Arridge IEEE TIP 2014, Ehrhardt et al Inverse Problems 2015, Gallardo and Meju Geophysical Research Letters 2003; Sochen et al IEEE TIP 1998

SLIDE 73

Evolution of Test Data

total variation quadratic PL linear PL

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt 2015

SLIDE 74

Generated Diffusion

Diffusivity of Parallel Level Sets

The derivative of Sβ with respect to u can be written as DSβ[u] = − div

K

∇u

.

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al Inverse Problems 2015

SLIDE 75

Generated Diffusion

Diffusivity of Parallel Level Sets

Let Rv be Gauge coordinates for v. Then the derivative of Sβ with respect to u can be written as DSβ[u] = − div

K

∇u

m.j.ehrhardt@damtp.cam.ac.uk

Ehrhardt et al Inverse Problems 2015

SLIDE 76

Generated Diffusion

Diffusivity of Parallel Level Sets

Let Rv be Gauge coordinates for v. Then the derivative of Sβ with respect to u can be written as DSβ[u] = − div

RvΛRv T∇u
with Λ = Diag(λ⊥, λ, . . . , λ).

◮ form of derivative independent of ϕ, ψ ◮ only λ⊥ and λ depend on ϕ and ψ m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al Inverse Problems 2015

SLIDE 77

Numerical Results

m.j.ehrhardt@damtp.cam.ac.uk

SLIDE 78

MRI sampling: full

Total Variation Quadratic PL Linear PL Total Variation Quadratic PL Linear PL

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al Inverse Problems 2015

SLIDE 79

MRI sampling: full

Total Variation Quadratic PL Linear PL Total Variation Quadratic PL Linear PL

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al Inverse Problems 2015

SLIDE 80

MRI sampling: 20 radial spokes

Total Variation Quadratic PL Linear PL Total Variation Quadratic PL Linear PL

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al Inverse Problems 2015

SLIDE 81

MRI sampling: 20 radial spokes

Total Variation Quadratic PL Linear PL Total Variation Quadratic PL Linear PL

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al Inverse Problems 2015

SLIDE 82

MRI sampling: 15 radial spokes

Total Variation Quadratic PL Linear PL Total Variation Quadratic PL Linear PL

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al Inverse Problems 2015

SLIDE 83

MRI sampling: 15 radial spokes

Total Variation Quadratic PL Linear PL Total Variation Quadratic PL Linear PL

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al Inverse Problems 2015

SLIDE 84

MRI sampling: uniform spiral

Total Variation Quadratic PL Linear PL Total Variation Quadratic PL Linear PL

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al Inverse Problems 2015

SLIDE 85

MRI sampling: uniform spiral

Total Variation Quadratic PL Linear PL Total Variation Quadratic PL Linear PL

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al Inverse Problems 2015

SLIDE 86

Quantitative Results

f u l l r a d i a l 2 r a d i a l 1 5 u n i . s p i r a l n

n
u

n i . s p i r a l 2nd l i n e 10 20 30 relative ℓ2-error [%] PET MLEM zerofill TV QPL LPL f u l l r a d i a l 2 r a d i a l 1 5 u n i . s p i r a l n

n
u

n i . s p i r a l 2nd l i n e 10 20 30 MRI

m.j.ehrhardt@damtp.cam.ac.uk Ehrhardt et al Inverse Problems 2015

SLIDE 87

Conclusions

◮ Part I: new prior incorporates anatomical structure

MLEM P

◮ Part II: Joint Reconstruction

◮ Parallel Level Set prior encodes joint structure ◮ Minimizing PLS yields structurally coupled anisotropic diffusion ◮ Combining two inverse problems can be beneficial to both