An ADMM algorithm for constrained material decomposition in spectral - - PowerPoint PPT Presentation

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

An ADMM algorithm for constrained material decomposition in spectral CT

May, 23th, 2018

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 1 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Spectral CT

Spectral CT records energy-dependant data

Log10 of data

100 200 300 400 500 600 50 100 150

Bin 1

3 4 5 6 100 200 300 400 500 600 50 100 150

Bin 2

3 4 5 6 100 200 300 400 500 600 50 100 150

Bin 3

3 4 5 6 100 200 300 400 500 600 50 100 150

Bin 4

3 3.5 4 4.5 5 5.5

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 2 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Spectral CT

Spectral CT records energy-dependant data

Log10 of data

100 200 300 400 500 600 50 100 150

Bin 1

3 4 5 6 100 200 300 400 500 600 50 100 150

Bin 2

3 4 5 6 100 200 300 400 500 600 50 100 150

Bin 3

3 4 5 6 100 200 300 400 500 600 50 100 150

Bin 4

3 3.5 4 4.5 5 5.5

Attenuation coefficients depend on energy

Energy

Attenuation Soft tissues Bones Gadolinium

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 2 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Spectral CT

Spectral CT records energy-dependant data

Log10 of data

100 200 300 400 500 600 50 100 150

Bin 1

3 4 5 6 100 200 300 400 500 600 50 100 150

Bin 2

3 4 5 6 100 200 300 400 500 600 50 100 150

Bin 3

3 4 5 6 100 200 300 400 500 600 50 100 150

Bin 4

3 3.5 4 4.5 5 5.5

Allow decomposing materials

Material maps

100 200 300 400 500 600 50 100 150

Soft tissues

10 20 30 100 200 300 400 500 600 50 100 150

Bones

2 4 6 8 10 100 200 300 400 500 600 50 100 150

Gadolinium

0.05 0.1 0.15 0.2

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 2 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Spectral CT

Spectral CT records energy-dependant data

Log10 of data

100 200 300 400 500 600 50 100 150

Bin 1

3 4 5 6 100 200 300 400 500 600 50 100 150

Bin 2

3 4 5 6 100 200 300 400 500 600 50 100 150

Bin 3

3 4 5 6 100 200 300 400 500 600 50 100 150

Bin 4

3 3.5 4 4.5 5 5.5

Solve the non-convex

− − − − − − − − − − − − − − − →

nonlinear problem s=F(a)

Allow decomposing materials

Material maps

100 200 300 400 500 600 50 100 150

Soft tissues

10 20 30 100 200 300 400 500 600 50 100 150

Bones

2 4 6 8 10 100 200 300 400 500 600 50 100 150

Gadolinium

0.05 0.1 0.15 0.2

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 2 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Introduction

Previous work A Gauss-Newton algorithm was developed [1] + Fast convergence + Better decomposition error than first order algorithms + Spatially regularized – Can have negative values

Soft tissues

50 100 150 200 250 300 20 40 60 80

Gauss-Newton

10 20 30

Bones

50 100 150 200 250 300 20 40 60 80 2 4 6

Gadolinium

50 100 150 200 250 300 20 40 60 80 0.05 0.1 0.15 0.2

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 3 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Introduction

Previous work A Gauss-Newton algorithm was developed [1] + Fast convergence + Better decomposition error than first order algorithms + Spatially regularized – Can have negative values

Soft tissues

50 100 150 200 250 300 20 40 60 80

Gauss-Newton

10 20 30

Bones

50 100 150 200 250 300 20 40 60 80 2 4 6

Gadolinium

50 100 150 200 250 300 20 40 60 80 0.05 0.1 0.15 0.2

Constrained algorithms Positivity constraints were added for materials decomposition Long and Fessler 2014 [2] Noh and Fessler 2009 [3] Sidky and Pan 2014 [4] Ding and Long 2017 [5]

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 3 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Introduction

Previous work A Gauss-Newton algorithm was developed [1] + Fast convergence + Better decomposition error than first order algorithms + Spatially regularized – Can have negative values

Soft tissues

50 100 150 200 250 300 20 40 60 80

Gauss-Newton

10 20 30

Bones

50 100 150 200 250 300 20 40 60 80 2 4 6

Gadolinium

50 100 150 200 250 300 20 40 60 80 0.05 0.1 0.15 0.2

Constrained algorithms Positivity constraints were added for materials decomposition Long and Fessler 2014 [2] Noh and Fessler 2009 [3] Sidky and Pan 2014 [4]    First oder algorithms in projection and image domain Ding and Long 2017 [5]

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 3 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Introduction

Previous work A Gauss-Newton algorithm was developed [1] + Fast convergence + Better decomposition error than first order algorithms + Spatially regularized – Can have negative values

Soft tissues

50 100 150 200 250 300 20 40 60 80

Gauss-Newton

10 20 30

Bones

50 100 150 200 250 300 20 40 60 80 2 4 6

Gadolinium

50 100 150 200 250 300 20 40 60 80 0.05 0.1 0.15 0.2

Constrained algorithms Positivity constraints were added for materials decomposition Long and Fessler 2014 [2] Noh and Fessler 2009 [3] Sidky and Pan 2014 [4]    First oder algorithms in projection and image domain Ding and Long 2017 [5] } Second order algorithm in image domain

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 3 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Proposed method

Positivity / Inequality We propose to enforce the positivity of the materials maps

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 4 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Proposed method

Positivity / Inequality We propose to enforce the positivity of the materials maps Equality Quantity of injected marker is known Should be retrieved in the decompositions

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 4 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Proposed method

Positivity / Inequality We propose to enforce the positivity of the materials maps Equality Quantity of injected marker is known Should be retrieved in the decompositions Algorithm Use of an alternating direction method of multipliers (ADMM) : update of a with a second order algorithm

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 4 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Forward problem

Denoting s ∈ RIP the data vector and a ∈ RMP the materials maps, we have : s = [s1,1, . . . , sI,1, . . . , sI,P]T a = [a1,1, ..., am,1, ..., aM,P]T

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 5 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Forward problem

Denoting s ∈ RIP the data vector and a ∈ RMP the materials maps, we have : s = [s1,1, . . . , sI,1, . . . , sI,P]T a = [a1,1, ..., am,1, ..., aM,P]T Using the standard spectral CT forward model is : si,p =

• R

n0

i (E) exp

• −

M

• m=1

am,pτm(E)

• dE

with n0

i (E) the effective spectrum, τm(E) a function representing the material

attenuation and am,p the projected mass of material m on pixel p.

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 5 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Inverse problem

We propose to solve : min

a

C(a, s) s.t

• a ≥ 0
• p am,p = cm

(1) with cm is the quantity of mth material. We chose C(a, s) = D(a, s) + αRR(a) where D(a, s) = ||s − F(a)||2

W

R(a) = ||∆asoft||2

2 + ||∇abone||1 + ||∇aGd||1 Lake Como Summer School Tom Hohweiller CREATIS Laboratory 6 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

ADMM

The Lagrangian function is minimized by an ADMM algorithm : L(a, b,α α αI, αE, s) = D(a, s) + αRR(a) + HE(a, αE) + GE(a) + HI(a, b,α α αI) + GI(a, b) + 1(b)

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 7 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

ADMM

The Lagrangian function is minimized by an ADMM algorithm : L(a, b,α α αI, αE, s) = D(a, s) + αRR(a) + HE(a, αE) + GE(a) + HI(a, b,α α αI) + GI(a, b) + 1(b) where Equality :

p am,p = cm

HE(a, αE) = αE  

P

• p=1

agd,p cgd − 1   GE(a) = βE 2  

P

• p=1

agd,p cgd − 1  

2

with αE ∈ R and βE ∈ R.

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 7 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

ADMM

The Lagrangian function is minimized by an ADMM algorithm : L(a, b,α α αI, αE, s) = D(a, s) + αRR(a) + HE(a, αE) + GE(a) + HI(a, b,α α αI) + GI(a, b) + 1(b) where Inequality : a ≥ 0 Inequalities need an auxiliary variable, b, such that the constraint becomes a = b and b ≥ 0. HI(a, b,α α αI) = α α αT

I (b − a)

GI(a, b) = βI 2 ||b − a||2

2

and 1(b) =

• 0,

if b ≥ 0 ∞,

• therwise

with αI ∈ RMP and βI ∈ R.

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 7 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Minimization scheme

Minimization of L(a, b,α α αI, αE, s) is done by finding the saddle point [6], through alternating updates : Update Algorithm used aℓ+1 ∈ argmin

a

L(a, bℓ, αℓ

E,α

α αℓ

I)

Iterative Gauss-Newton bℓ+1 ∈ argmin

b

L(aℓ+1, b, αℓ

E,α

α αℓ

I)

Proximal algorithm αℓ+1

E

∈ argmax

αE

L(aℓ+1, bℓ+1, αE,α α αℓ

I)

Ascent gradient α α αℓ+1

I

∈ argmax

α α αI

L(aℓ+1, bℓ+1, αℓ+1

E

,α α αI) Ascent gradient Moreover, at each iteration ℓ, Lagrangian parameters are increasing βℓ+1

E

= ωβℓ

E

βℓ+1

I

= ωβℓ

I

Until they reach the maximum allow value βmax.

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 8 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Parameters

Phantom Thorax : 84 × 306 pixels 3 materials : soft tissues, bones and gadolinium (portal vein) 4 energy bins Algorithm a0 = 0g.cm−2 αR = 1, α0

E = 0, α

α α0

I = 0

β0

E = 100, β0 I = 10−2, ω = 1.5

βmax = 1010 Stopping criteria

The inner loop stops if

the relative decrease of the cost function is too low (10−3) it made more than 30 iterations

The outer loop stops if

the equality is respected : | P

p aℓ

gd,p

cgd − 1

• | < 10−3

the inequality is respected : ||aℓ − bℓ|| < 10−3

Decomposition for θ = [0 ˚ , 359 ˚ ]

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 9 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Projection domain

Ground truth

Soft tissues

5 10 15 20 25 30 35 40

Bones

2 4 6 8 10 12 14 16

Gadolinium

0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4

Unconstrained algorithm (GN)

5 10 15 20 25 30 35 40

• 2

2 4 6 8

• 0.05

0.05 0.1 0.15 0.2 0.25 0.3

Constrained algorithm (ADMM)

5 10 15 20 25 30 35 40 1 2 3 4 5 6 7 8 0.05 0.1 0.15 0.2 0.25 0.3

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 10 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Projection domain

50 100 150 200 250 300 350 60 80 100 120 140 160 Evaluation of the decomposition with the ℓ2-norm : ξθ

D = M

• m

||aθ

m − atruth,θ m

||2

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 11 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Image domain

Ground truth

Soft tissues

Bones Gadolinium

Unconstrained algorithm (GN) Constrained algorithm (ADMM)

0.05 0.1 0.15 0.2 0.25 0.3

• 0.05

0.05 0.1 0.15 0.2 0.25 0.3 0.35 5 10 15 10-3

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 12 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Image domain

10 20 30 40 50 60 70 80 1 2 3 4 5 6 7 Evaluation of the reconstructed materials with the ℓ2-norm : χz

D = M

• m

||ρ ρ ρz

m − ρ

ρ ρtruth,z

m

||2

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 13 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Conclusion

Equality and inequality constraints give a better decomposition and visualization Decrease cross-talk However, it needs more iteration

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 14 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Conclusion

Equality and inequality constraints give a better decomposition and visualization Decrease cross-talk However, it needs more iteration And next : Application of this algorithm to real data

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 14 / 15

Introduction Material decomposition Constrained algorithm Numerical experiments Results Conclusion

Thank you for your attention ! Questions ? :)

References Nicolas Ducros, Juan Felipe Perez-Juste Abascal, Bruno Sixou, Simon Rit, and Françoise Peyrin, “Regularization of nonlinear decomposition of spectral x-ray projection images,” Medical physics, vol. 44, no. 9, 2017. Yong Long and Jeffrey A Fessler, “Multi-material decomposition using statistical image reconstruction for spectral ct,” IEEE transactions on medical imaging, vol. 33, no. 8, pp. 1614–1626, 2014. Joonki Noh, Jeffrey A Fessler, and Paul E Kinahan, “Statistical sinogram restoration in dual-energy ct for pet attenuation correction,” IEEE transactions on medical imaging, vol. 28, no. 11, pp. 1688–1702, 2009. Xiaochuan Pan, Buxin Chen, Zheng Zhang, Erik Pearson, Emil Sidky, and Xiao Han, “Optimization-based reconstruction exploiting spectral information in ct,” in The Third International Conference on Image Formation in X-Ray Computed Tomography, 2014, pp. 228–232. Qiaoqiao Ding, Tianye Niu, Xiaoqun Zhang, and Yong Long, “Image-domain multi-material decomposition for dual-energy ct based on correlation and sparsity of material images,” arXiv preprint arXiv :1710.07028, 2017. Jonathan Eckstein and W Yao, “Augmented lagrangian and alternating direction methods for convex optimization : A tutorial and some illustrative computational results,” RUTCOR Research Reports, vol. 32, pp. 3, 2012.

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 15 / 15

Minimization scheme

while Stopping criteria aren’t met for the

• uter loop do

while Stopping criteria aren’t met for the inner loop do solve

• Hkδak = −gk

ak+1 = ak + λkδak k = k + 1 end aℓ+1 = ak−1

Update of a : Gauss-Newton algorithm

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 1 / 1

Minimization scheme

while Stopping criteria aren’t met for the

• uter loop do

while Stopping criteria aren’t met for the inner loop do solve

• Hkδak = −gk

ak+1 = ak + λkδak k = k + 1 end aℓ+1 = ak−1 bℓ+1 = Π(aℓ+1 − α

α αℓ

I

βI )

Update of a : Gauss-Newton algorithm Update of b : Proximal algorithm

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 1 / 1

Minimization scheme

while Stopping criteria aren’t met for the

• uter loop do

while Stopping criteria aren’t met for the inner loop do solve

• Hkδak = −gk

ak+1 = ak + λkδak k = k + 1 end aℓ+1 = ak−1 bℓ+1 = Π(aℓ+1 − α

α αℓ

I

βI )

αℓ+1

E

= αℓ

E + βE HE (aℓ+1

m

) αℓ

E

Update of a : Gauss-Newton algorithm Update of b : Proximal algorithm Update of αE : Ascent gradient

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 1 / 1

Minimization scheme

while Stopping criteria aren’t met for the

• uter loop do

while Stopping criteria aren’t met for the inner loop do solve

• Hkδak = −gk

ak+1 = ak + λkδak k = k + 1 end aℓ+1 = ak−1 bℓ+1 = Π(aℓ+1 − α

α αℓ

I

βI )

αℓ+1

E

= αℓ

E + βE HE (aℓ+1

m

) αℓ

E

α α αℓ+1

I

= α α αℓ

I + βI(bℓ+1 − aℓ+1)

Update of a : Gauss-Newton algorithm Update of b : Proximal algorithm Update of αE : Ascent gradient Update of α α αI : Ascent gradient

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 1 / 1

Minimization scheme

while Stopping criteria aren’t met for the

• uter loop do

while Stopping criteria aren’t met for the inner loop do solve

• Hkδak = −gk

ak+1 = ak + λkδak k = k + 1 end aℓ+1 = ak−1 bℓ+1 = Π(aℓ+1 − α

α αℓ

I

βI )

αℓ+1

E

= αℓ

E + βE HE (aℓ+1

m

) αℓ

E

α α αℓ+1

I

= α α αℓ

I + βI(bℓ+1 − aℓ+1)

βℓ+1

E

= ωβℓ

E

βℓ+1

I

= ωβℓ

I

Update of a : Gauss-Newton algorithm Update of b : Proximal algorithm Update of αE : Ascent gradient Update of α α αI : Ascent gradient Update of the βs

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 1 / 1

Minimization scheme

while Stopping criteria aren’t met for the

• uter loop do

while Stopping criteria aren’t met for the inner loop do solve

• Hkδak = −gk

ak+1 = ak + λkδak k = k + 1 end aℓ+1 = ak−1 bℓ+1 = Π(aℓ+1 − α

α αℓ

I

βI )

αℓ+1

E

= αℓ

E + βE HE (aℓ+1

m

) αℓ

E

α α αℓ+1

I

= α α αℓ

I + βI(bℓ+1 − aℓ+1)

βℓ+1

E

= ωβℓ

E

βℓ+1

I

= ωβℓ

I

if βℓ+1

E

> βmax then βℓ+1

E

= βmax end if βℓ+1

I

> βmax then βℓ+1

I

= βmax end ℓ = ℓ + 1 end

Update of a : Gauss-Newton algorithm Update of b : Proximal algorithm Update of αE : Ascent gradient Update of α α αI : Ascent gradient Update of the βs Make sure that the values aren’t too high

Lake Como Summer School Tom Hohweiller CREATIS Laboratory 1 / 1