Empirical Phase Transitions in Sparsity-Regularized Computed - - PowerPoint PPT Presentation

Empirical Phase Transitions in Sparsity-Regularized Computed Tomography

Jakob Sauer Jørgensen

Postdoc, DTU

Sparse Tomo Days, DTU, March 28, 2014

Joint work with Per Christian Hansen, DTU Emil Sidky and Xiaochuan Pan, U. Chicago Christian Kruschel and Dirk Lorenz, TU Braunschweig

Exploiting prior knowledge in CT

Discrete imaging model: Ax = b Typical CT images:

◮ Regions of homogeneous tissue. ◮ Separated by sharp boundaries.

Reconstruction by regularization: x⋆ = argmin

x

D(Ax, b) + λ · R(x) Sparsity-promoting choices:

◮ R(x) = x1

(ℓ1/basis pursuit)

◮ R(x) = xTV

(total variation)

◮ R(x) = DT x1

(analysis-ℓ1)

2 / 24

TV example: Physical head phantom, CB-CT

(Bian 2010). Courtesy of X. Pan, U. Chicago.

3 / 24

TV example: Human coronary artery, CB-CT

Courtesy of X. Pan, U. Chicago. Data collected with a bench-top CB-CT of Dr. E. Ritman at Mayo

4 / 24

Less successful CT cases for TV:

(Herman & Davidi, 2008) Original TV (J. et al, 2011) (Courtesy of S. Soltani, DTU)

5 / 24

Lack of quantitative understanding

Some fundamental questions remain unanswered:

◮ Under what conditions will reconstruction work? ◮ Robustness to noise? ◮ Which types of images? ◮ What is sufficient sampling?

Application-specific vs. general

◮ Focus on the imaging model.

6 / 24

Classical CT sampling results

Continuous image and data:

◮ Based on invertibility and stability of Radon transform etc. ◮ Fan-beam: 180◦ plus fan-angle ◮ Cone-beam: Tuy’s condition

Discrete data:

◮ Nyquist sampling ◮ Assumption of bandlimited signal ◮ (Huesmann 1977, Natterer)

Reconstruction with sparse/compressible signal assumption?

7 / 24

Compressed Sensing

Guarantees of accurate reconstruction

◮ Under suitable assumptions, a sufficiently sparse signal can be

recovered from few measurements by ℓ1-minimization.

◮ RIP, incoherence, spark, ...

For tomography?

◮ So far no practically useful guarantees. ◮ Results for certain discrete tomography cases (Petra et al.)

This study:

◮ Empirical study of sampling conditions for tomographic

reconstruction of sparse signals

◮ Recoverability of single images ◮ Worst-case vs. average case

8 / 24

Reconstruction problems

Inequality-constrained regularization: x⋆ = argmin

x

R(x) s.t. Ax − b2 ≤ ǫ Simplified reconstruction problems:

finite-difference approximation of gradient

❆ ❯

BP xBP = argmin

x

x1 s.t. Ax = b ATV xATV = argmin

x

DT x1 s.t. Ax = b Algorithms:

◮ Our interest: Reliably obtaining accurate solution, not speed. ◮ Recast as linear programs (LPs) and solve by MOSEK.

9 / 24

Non-uniqueness of solutions

Both BP and ATV can have multiple solutions for same data:

◮ 1-norm convex, but not strictly convex. ◮ Even if xorig is a minimizer, others may exist.

Consequences:

◮ Different algorithms may produce different solutions. ◮ Decision of recoverability of xorig is algorithm-dependent.

Alternative idea:

◮ Can we test for uniqueness of solution?

10 / 24

Uniqueness test for BP

Given:

◮ b = Axorig ◮ I = support(xorig) ◮ AI is A with columns I

Characterization of solution uniqueness:

◮ xorig

uniquely minimizes minx x1 s.t. Ax = b if and only if

◮ AI is injective, and ◮ ∃w :

AT

I w = sign(xorig)I

and AT

Icw∞ < 1 (Plumbley 2007, Fuchs 2004, Grasmair et al. 2011)

11 / 24

Uniqueness test for ATV

Given:

◮ b = Axorig ◮ I = support(DT xorig) ◮ DI is D with columns I

Characterization of solution uniqueness:

◮ xorig

uniquely minimizes minx DT x1 s.t. Ax = b if and only if

◮ A DT

Ic

• is injective, and

◮ ∃w, v :

Dv = AT w, vI = sign(DT

I xorig),

vIc∞ < 1

Application of (Haltmeier 2013)

12 / 24

Uniqueness testing procedure using LP

• 1. Check

injectivity:

• 2. Solve LP:
• 3. Unique iff:

BP AI t⋆ = argmin t −te ≤ AT

Icw ≤ te

AT

I w = sign(xorig)I

t⋆ < 1 ATV A DT

Ic

• t⋆ = argmin t

−te ≤ vIc ≤ te Dv = AT w vI = sign(DT

I xorig)

t⋆ < 1

13 / 24

The geometry and system matrix

◮ Disk-shaped image inscribed in

Nside × Nside square.

◮ Number of pixels:

N ≈ π 4 N 2

side

◮ Fan-beam, equi-angular views

(Nviews = 3 shown)

◮ Number of rays per view: 2Nside ◮ System matrix A size:

M = Nviews · 2Nside A N Elements Aij computed by the line intersection method (implementation: www.imm.dtu.dk/~pch/AIRtools/)

14 / 24

BP image class examples images

spikes 1-power 2-power Relative sparsity 0.1 0.3 0.5 0.7 0.9

15 / 24

Reconstruction error vs. sampling and sparsity

10 20 30 40 10

−10

10

−5

10 Number of views Nv Reconstruction error ε = 10−4 Full rank 20% nonzeros 40% nonzeros 60% nonzeros 80% nonzeros

16 / 24

Phase diagrams: spikes with BP

Reconstruction Uniqueness test

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 Relative sampling: µ = Nv / Nv

suf

Relative sparsity: κ = k / N 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 Relative sampling: µ = Nv / Nv

suf

Relative sparsity: κ = k / N 0.2 0.4 0.6 0.8 1

◮ Fraction recovered/unique of 100 instances at each point

(κ, µ) of relative sparsity and sampling.

◮ Excellent agreement of reconstruction and uniqueness test. ◮ Well-separated “no-recovery” and “full-recovery” regions. ◮ Phase transition as in compressed sensing (Donoho-Tanner).

17 / 24

Comparing image classes, BP

spikes 1-power 2-power

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 Relative sampling: µ = Nv / Nv

suf

Relative sparsity: κ = k / N 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 Relative sampling: µ = Nv / Nv

suf

Relative sparsity: κ = k / N 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 Relative sampling: µ = Nv / Nv

suf

Relative sparsity: κ = k / N 0.2 0.4 0.6 0.8 1

Average sufficient sampling

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 Relative sparsity: κ = k / N Relative sampling: µ = Nv / Nv

suf

spikes 1−power 2−power

18 / 24

Example images: altproj, trununif

altproj trununif Relative sparsity 0.1 0.3 0.5 0.7 0.9

19 / 24

Phase diagrams: altproj with ATV

Reconstruction Uniqueness test

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0.5 1

v

Relative sparsity: κ = k / N 0.5 1 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0.5 1

v

Relative sparsity: κ = k / N 0.5 1

◮ Fraction recovered/unique of 100 instances at each point

(κ, µ) of relative sparsity and sampling.

◮ Excellent agreement of reconstruction and uniqueness test. ◮ Well-separated “no-recovery” and “full-recovery” regions. ◮ Phase transition as in compressed sensing (Donoho-Tanner).

20 / 24

Comparing image classes, ATV

altproj

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0.5 1 Relative sampling: µ = Nv / N Relative sparsity: κ = k / N 0.5 1

trununif

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0.5 1 Relative sampling: µ = Nv / N Relative sparsity: κ = k / N 0.5 1

Average sufficient sampling

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0.2 0.4 0.6 0.8 1 Relative sparsity: κ = k / N Relative sampling: µ = Nv / Nv

suf

altproj trununif40

21 / 24

Time: Reconstruction vs. uniqueness test

0.2 0.4 0.6 0.8 1 10

−1

10 10

1

10

2

spfrac time in secs R 5 R 13 R 21 UT 5 UT 13 UT 21 0.5 1 1.5 2 10

−1

10 10

1

10

2

spfrac time in secs R 5 R 13 R 21 UT 5 UT 13 UT 21

ATV BP

◮ 10 repetitions at each relative sparsity and 5, 13, 21 views. ◮ Comparable time of reconstruction (R) and uniqueness test (UT).

22 / 24

A more well-known image: Shepp-Logan on disk

BP ATV xorig0 = 1988 NBP

v

= 17 Dxorig0 = 610 NATV

v

= 8

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 Relative sparsity: κ = k / N Relative sampling: µ = Nv / Nv

suf

spikes 1−power 2−power

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0.2 0.4 0.6 0.8 1 Relative sparsity: κ = k / N

v v

altproj trununif40

23 / 24

Conclusions and future work

Conclusions

◮ Empirical evidence of relation between sparsity and sampling ◮ Reconstruction and uniqueness test ◮ Phase transition from no to full recovery ◮ Small dependence on image class, mostly sparsity ◮ Additional results (not shown): limited angle, robustness to

noise, scaling with image size. Future work and open questions

◮ Extensions: Isotropic TV, more realistic image classes, ... ◮ Theoretical/compressed sensing explanation? ◮ Connection to classical CT sampling results?

24 / 24