Xray-CT image reconstruction from few projections

Slides for PhD defense

Transcript

1. Introduction Blob Simulations Real data Perspectives Background TV Besov CS 1 Literature review and Compressed Sensing theory Background Regularization by TV minimization Regularization by Besov minimization Compressed Sensing and CT 2 Blob Basis for CS and CT General properties of blob Single scale Gaussian blob image Tight frame by a multiscale blob system Multiscale Mexican hat blob 3 2D Numerical Experiments - Simulated data TV reconstruction results `1 reconstruction results TV-`1 reconstruction results 4 2D Numerical Experiments - Real data 5 Conclusion and Persperctives 2 / 79
2. Introduction Blob Simulations Real data Perspectives Background TV Besov CS Outline 1 Literature review and Compressed Sensing theory Background Regularization by TV minimization Regularization by Besov minimization Compressed Sensing and CT 2 Blob Basis for CS and CT General properties of blob Single scale Gaussian blob image Tight frame by a multiscale blob system Multiscale Mexican hat blob 3 2D Numerical Experiments - Simulated data TV reconstruction results `1 reconstruction results TV-`1 reconstruction results 4 2D Numerical Experiments - Real data 5 Conclusion and Persperctives 3 / 79
3. Introduction Blob Simulations Real data Perspectives Background TV Besov CS Section 1 : Literature review and Compressed Sensing theory 1 Literature review and Compressed Sensing theory Background Regularization by TV minimization Regularization by Besov minimization Compressed Sensing and CT 2 Blob Basis for CS and CT General properties of blob Single scale Gaussian blob image Tight frame by a multiscale blob system Multiscale Mexican hat blob 3 2D Numerical Experiments - Simulated data TV reconstruction results `1 reconstruction results TV-`1 reconstruction results 4 2D Numerical Experiments - Real data 5 Conclusion and Persperctives 4 / 79
4. Introduction Blob Simulations Real data Perspectives Background TV Besov CS X-ray CT reconstruction X-ray transform P : L2(IRd) ! L2(T ) Pf(✓, y) = P✓f(y) 4 = Z I R f(y + t✓)dt Object domain : f 2 X ✓ L2(IRd) Sinogram domain : Pf 2 L2(T ), T = {(✓, y), ✓ 2 Sd 1, y 2 ✓?} Inversion of P by Filtered Back-Projection (FBP) : f = P⇤I 1Pf Discrete X-ray projector A : X ! IRM M = P ⇥ D, P projections and D detector bins. Inversion of Af = b for essentially B-band-limited function f Shannon sampling condition : one needs P B/⇡ projections, D B detector bins. 5 / 79
5. Introduction Blob Simulations Real data Perspectives Background TV Besov CS Classical algorithms : instability and artifacts for P small Minimum square error (MSE) reconstruction f⇤ = arg min f2X kAf bk2 Solved by Conjugate Gradient (CG) with 250 iterations. P = 384 projections : (a) SNR(b)=50dB, SNR=14.19dB (b) SNR(b)=25dB, SNR=9.41dB 6 / 79
6. Introduction Blob Simulations Real data Perspectives Background TV Besov CS Classical algorithms : instability and artifacts for P small Minimum square error (MSE) reconstruction f⇤ = arg min f2X kAf bk2 Solved by Conjugate Gradient (CG) with 250 iterations. P = 64 projections : (c) SNR(b)=50dB, SNR=7.18dB (d) SNR(b)=25dB, SNR=-9.86dB 7 / 79
7. Introduction Blob Simulations Real data Perspectives Background TV Besov CS Few projections problem Problem Reconstruct a d-dimensional (d = 2, 3) object f from few projections (P ⌧ B/⇡) : Study the feasibility on generic object f Provide reconstruction algorithms validated on simulated and real data Linear data model b = Af + n, with ni ⇠ N(0, m 1 i ), mi : mean photon number Regularization in a continuous space Maximum A Posteriori estimation : ˆ fMAP = arg min f2X 1 2 kAf bk2 + R(f) Equivalent form : ˆ fMAP = ( arg minf2X R(f) s.t. kAf bk2  ✏ Discretization of the solution space X X ⇢ L2(IRd) is spanned by a Riesz basis or a frame, e.g. : f(x) = X k2Z Zd fk (x xk) Finite dimension : f = [f1 . . . , fN ] A : IRN ! IRM , with Am,n = R Lm n(x)dx 8 / 79
8. Introduction Blob Simulations Real data Perspectives Background TV Besov CS Few projections problem Problem Reconstruct a d-dimensional (d = 2, 3) object f from few projections (P ⌧ B/⇡) : Study the feasibility on generic object f Provide reconstruction algorithms validated on simulated and real data Linear data model b = Af + n, with ni ⇠ N(0, m 1 i ), mi : mean photon number Regularization in a discrete space Maximum A Posteriori estimation : ˆ fMAP = arg min f2I RN 1 2 kAf bk2 + R(f) Equivalent form : ˆ fMAP = ( arg minf2I RN R(f) s.t. kAf bk2  ✏ Discretization of the solution space X X ⇢ L2(IRd) is spanned by a Riesz basis or a frame, e.g. : f(x) = X k2Z Zd fk (x xk) Finite dimension : f = [f1 . . . , fN ] A : IRN ! IRM , with Am,n = R Lm n(x)dx 9 / 79
9. Introduction Blob Simulations Real data Perspectives Background TV Besov CS Existence of solution Admissible priors R(·) 1 convex : R( f + (1 )g)  R(f) + (1 )R(g) 2 proper : {f|R(f) = 1} 6= X 3 lower semi-continuous (l.s.c.) : lim infn R(fn) R(f), for fn ! f 4 homogeneous : R( f) = | |R(f) 5 R(f) > 0 for any f 2 Ker A, f 6= 0 Proposition There exists a global minimizer to the optimization problem ˆ fMAP = arg min f2I RN 1 2 kAf bk2 + R(f) Examples of popular regularizations Tikhonov regularization : R(f) = 1kfk2 + 2 P i,j krfi,jk2 Edge Preserving Prior : R(f) = P i,j '((D1f)i,j) + '((D2f)i,j) Regularization in BV or Besov space 10 / 79
10. Introduction Blob Simulations Real data Perspectives Background TV Besov CS Section 1 : Literature review and Compressed Sensing theory 1 Literature review and Compressed Sensing theory Background Regularization by TV minimization Regularization by Besov minimization Compressed Sensing and CT 2 Blob Basis for CS and CT General properties of blob Single scale Gaussian blob image Tight frame by a multiscale blob system Multiscale Mexican hat blob 3 2D Numerical Experiments - Simulated data TV reconstruction results `1 reconstruction results TV-`1 reconstruction results 4 2D Numerical Experiments - Real data 5 Conclusion and Persperctives 11 / 79
11. Introduction Blob Simulations Real data Perspectives Background TV Besov CS BV space and regularization by TV minimization Bounded Variation space BV (⌦) for CT object f 2 BV (⌦) , kfkBV 4 = Z ⌦ |f| + TV (f) < 1 2D discrete TV norm of a pixel image kfkT V 4 = n X i=1 n X j=1 q (fi,j fi+1,j)2 + (fi,j fi,j+1)2 Regularization by TV minimization min f2I RN 1 2 kAf bk2 + kfkT V , or min f2I RN kfkT V s.t. kAf bk2  " Minimizing the number of active pixels For f 2 BV (⌦), proportion of active pixels . n 2kfkT V ' n 1kfkT V ! 0 0 1 12 / 79
12. Introduction Blob Simulations Real data Perspectives Background TV Besov CS Reconstruction by TV minimization Solved by TVAL3 algorithm. P = 64 Projections. (a) SNR=17.68dB, SNR(b)=50dB (b) SNR=15.26dB, SNR(b)=25dB 13 / 79
13. Introduction Blob Simulations Real data Perspectives Background TV Besov CS Section 1 : Literature review and Compressed Sensing theory 1 Literature review and Compressed Sensing theory Background Regularization by TV minimization Regularization by Besov minimization Compressed Sensing and CT 2 Blob Basis for CS and CT General properties of blob Single scale Gaussian blob image Tight frame by a multiscale blob system Multiscale Mexican hat blob 3 2D Numerical Experiments - Simulated data TV reconstruction results `1 reconstruction results TV-`1 reconstruction results 4 2D Numerical Experiments - Real data 5 Conclusion and Persperctives 14 / 79
14. Introduction Blob Simulations Real data Perspectives Background TV Besov CS Besov space and wavelet basis Besov space B↵ p,q = B↵ q (Lp(IRd)) 0 < ↵ < 1 : order of smoothness 1  p, q  1 : tuning parameters, p = q to simplify Let c(l) j,k = hf, (l) j,k i the wavelet coe cient. Besov norm equivalence (s = ↵ + d(1/2 1/p)) : kfkp B↵ p,p ⇠ X k2Z Zd c(0) j0,k p + X j j0 X k2Z Zd 2d 1 X l=1 2jps c(l) j,k p 4 = kckp p,↵ Natural images are compressible (p < 2), and the energy is concentrated on a few large wavelet coe cients. Discretizing f on Cartesian lattice : f = [f1, . . . , fN ]> Discrete Wavelet Transform (DWT) : f = Wc Regularization by Besov norm minimization with p < 2 min c 1 2 kAWc bk2 + kckp p,↵ , for c 2 IRN Di↵erentiable for p > 1, solved by nonlinear CG. 15 / 79
15. Introduction Blob Simulations Real data Perspectives Background TV Besov CS Reconstruction by Besov norm minimization Using Daubechies-6 wavelet, SNR(b)=50db, p = 1.1, ↵ ' 0.81. (a) P =128, SNR=17.79dB (b) P =64, SNR=14.21dB Pseudo Gibbs artifacts are amplified for small number of projection. 16 / 79
16. Introduction Blob Simulations Real data Perspectives Background TV Besov CS Comparison with TV minimization Zoomed view on ROI. P=128 projections. (a) TV, SNR=14.61dB (b) Db6, SNR=17.79dB TV : ”Texture killer”, low contrast regions are smoothed out. Besov : Preservation of low contrast regions, no stairwise artifacts of TV. 17 / 79
17. Introduction Blob Simulations Real data Perspectives Background TV Besov CS Comparison with TV minimization Zoomed view on ROI. P=128 projections. (a) TV, SNR=14.61dB (b) Db6, SNR=17.79dB TV : ”Texture killer”, low contrast regions are smoothed out. Besov : Preservation of low contrast regions, no stairwise artifacts of TV. 18 / 79
18. Introduction Blob Simulations Real data Perspectives Background TV Besov CS Section 1 : Literature review and Compressed Sensing theory 1 Literature review and Compressed Sensing theory Background Regularization by TV minimization Regularization by Besov minimization Compressed Sensing and CT 2 Blob Basis for CS and CT General properties of blob Single scale Gaussian blob image Tight frame by a multiscale blob system Multiscale Mexican hat blob 3 2D Numerical Experiments - Simulated data TV reconstruction results `1 reconstruction results TV-`1 reconstruction results 4 2D Numerical Experiments - Real data 5 Conclusion and Persperctives 19 / 79
19. Introduction Blob Simulations Real data Perspectives Background TV Besov CS Compressed Sensing : sparsity, measurements and `1 minimization Sparsity and best s-terms approximation of IRN vector s-sparse : number of non zero entries kxk 0  s Compressible : the n-th largest entry |x(n) | . n 1/p, with p < 2 Best s-terms approximation : xs 4 = arg minkyk0 s kx yk 2 Linear measurements A : IRN ! IRM b = Ax + n, with (Ax)m = ham, xi A satisfies RIP with constant 2s < 1 Minimal measurments : M & s log(N/M) Reconstruction by `1 minimization min y2I RN kyk1 s.t. kAy bk2  ✏ (P✏ 1 ) Error of `1 minimization If A has the RIP constant 2s < p 2 1, then the solution x⇤ of (P✏ 1 ) obeys : kx⇤ xk2  C0 kx xsk 1 p s + C1✏ Reconstruction by `1 minimization is stable and as good as the best s-terms approximation. 20 / 79
20. Introduction Blob Simulations Real data Perspectives Background TV Besov CS Compressed Sensing : construction of the sensing matrix A Random construction of a good sensing matrix A Take an orthogonal system U Construct A by selecting uniformly at random M rows of U With high probability, A has the RIP constant 2s < p 2 1, if : M & sµ2(U)(log N)4 with µ(U) 4 = p N maxm,n|Um,n|, the incoherence of U. Examples and applications U = discrete Fourier transform : CT/MRI image reconstructions. U = Wavelet basis : Inpainting and image restorations. 21 / 79
21. Introduction Blob Simulations Real data Perspectives Background TV Besov CS Implications in CT reconstruction TV minimization using discrete Fourier measurements A as projector min y kykT V , s.t. Ay = Ax (c) 30 random projections (d) SNR=29.10dB 22 / 79
22. Introduction Blob Simulations Real data Perspectives Background TV Besov CS X-ray projector’s encoding capacity Monte-Carlo tests on a realistic fan-beam projector A ⇠ M ⇥ N M = P ⇥ 128, N = 64 ⇥ 64. P 2 {32, . . . , 64} projections equally distributed on [0, 2⇡). (a) (P1) (b) (P1)+ Wavelet Small P su ces to capture the information of a sparse vector x. 23 / 79
23. Introduction Blob Simulations Real data Perspectives Background TV Besov CS State-of-the-arts CS framework for CT CS framework Sparse representation : x = Dc, D is a (redundant) system. Nonlinear algorithms : ( minc kck 1 s.t. kADc yk2  ✏ minxkxkT V s.t. kAx yk2  ✏ Sparser is x (wrt D), smaller is the reconstruction error. Why we use pixel as basis ? Most of the CS methods in the CT literature are based on pixel : Simple representations and calculations on computer. Natural interface with many fast transforms D : FFT, Cosinus, Wavelet, Curvelet, etc. Why we should not use pixel as basis in CT ? Poor space-frequency localization. Discretization or approximation of the projector A. Heavy computational charge of D and A. 24 / 79
24. Introduction Blob Simulations Real data Perspectives Properties Gaussian blob Tight frame Outline 1 Literature review and Compressed Sensing theory Background Regularization by TV minimization Regularization by Besov minimization Compressed Sensing and CT 2 Blob Basis for CS and CT General properties of blob Single scale Gaussian blob image Tight frame by a multiscale blob system Multiscale Mexican hat blob 3 2D Numerical Experiments - Simulated data TV reconstruction results `1 reconstruction results TV-`1 reconstruction results 4 2D Numerical Experiments - Real data 5 Conclusion and Persperctives 25 / 79
25. Introduction Blob Simulations Real data Perspectives Properties Gaussian blob Tight frame Section 2 : Blob Basis for CS and CT 1 Literature review and Compressed Sensing theory Background Regularization by TV minimization Regularization by Besov minimization Compressed Sensing and CT 2 Blob Basis for CS and CT General properties of blob Single scale Gaussian blob image Tight frame by a multiscale blob system Multiscale Mexican hat blob 3 2D Numerical Experiments - Simulated data TV reconstruction results `1 reconstruction results TV-`1 reconstruction results 4 2D Numerical Experiments - Real data 5 Conclusion and Persperctives 26 / 79
26. Introduction Blob Simulations Real data Perspectives Properties Gaussian blob Tight frame Radial blob Radial blob : IRd ! IR Radial symmetric : (x) = b (kxk), x 2 IRd 2 L1(IRd) \ L2(IRd) and C1 smooth Well localized in both frequency and space domain Fourier transform : ˆ(!) = bˆ(k!k), ! 2 IRd Gaussian family blob (x) = exp( ↵kxk2), ↵ > 0, and other examples Exponential decay and optimal spatial-frequency localization Infinite support but can be truncated in numerical implementation Simple analytical expressions of many transforms E cient numerical implementation Parameter ↵ is easily adjusted according to the bandwidth and the lattice’s sampling step 27 / 79
27. Introduction Blob Simulations Real data Perspectives Properties Gaussian blob Tight frame Image representation by blob Construct a shift invariant space X using blob(s) a blob or a family of blobs { j}j2Z Z the lattice(s) on which the blobs are translated Image model 1 : Shift invariant space by a single blob A Gaussian blob + a d-dimensional hexagonal lattice L : f(x) = X k2Z Zd fk (x xk), 8f 2 X Image model 2 : (Tight) frame by a multiscale blob system A family of blobs j + d-dimensional hexagonal lattices Lj : f(x) = X j2Z Z X k2Z Zd fj,k j(x xj k ), 8f 2 X Reconstruction of f , determination of the representation coe cients 28 / 79
28. Introduction Blob Simulations Real data Perspectives Properties Gaussian blob Tight frame Hexagonal lattice for blob image Lattice L G = [g1, . . . gd] ⇠ d ⇥ d, an invertible matrix Lattice generated by G : L 4 = {xk = Gk, k 2 Z Zd} Fundamental domain of L Dual lattice : L⇤ 4 = {xk = G >k, k 2 Z Zd} Hexagonal lattice in 2D h : sampling step g1 = [1, 0]>, g2 = [1/2, p 3/2]>, G = h[g1, g2] More e cient in sampling band-limited function than Cartesian lattice (86.6% sampling density) Higher angular frequency (six neighbors) and better image visual quality than the Cartesian lattice (c) h = 1 (d) h = 0.5 29 / 79
29. Introduction Blob Simulations Real data Perspectives Properties Gaussian blob Tight frame Blob-driven X-ray projector Abel transform of a radial function P✓ (y) = Ab (t) 4 = 2 Z +1 t b (r) r p r2 t2 dr Analytical expression of Ab for the Gaussian blob : Ab (t) = p ⇡/↵ exp( ↵t2) X-ray transform of a blob image f(x) = P k2Z Zd fk (x xk) P✓f(y) = X k2Z Zd fkAb (ky ⇧✓? xk | {z } xk hxk,✓i✓ k) X-ray projector A Amk the contribution of the k-th blob to the m-th detector bin : (Af)m = X k2Z Z Amkfk, (A⇤b)k = M X m=1 Akmbm Easy parallelization of A on GPU : Blob-driven projector 30 / 79
30. Introduction Blob Simulations Real data Perspectives Properties Gaussian blob Tight frame Control of bandwidth Singular Value Decomposition of the X-ray transform Pf = X n2I N nhf, univn Ill-posedness : n ! 0 with n un : high frequency oscillating functions for n large For f band-limited : hf, uni = h ˆ f, ˆ uni ⇠ 0 for n large The conditionning of P is improved if restricted on X of (essentially) band-limited functions. Control the bandwidth of f by basis function ˆ f(!) = ˆ(!) X k2Z Zd fk exp( 2⇡ihxk, !i) Poor spatial-frequency localization of pixel basis Large pixel is needed to control the bandwith, while the L2 approximation error kf fhk decays as O(h) for f 2 H1. 31 / 79
31. Introduction Blob Simulations Real data Perspectives Properties Gaussian blob Tight frame Pixel vs. Gaussian blob MSE reconstructions by CG using 100 iterations. (a) Pixel : 384 ⇥ 384 (b) Blob : ⇠ 3842 blobs Reconstruction using small pixel is unstable. 32 / 79
32. Introduction Blob Simulations Real data Perspectives Properties Gaussian blob Tight frame Pixel vs. Gaussian blob MSE reconstructions by CG using 100 iterations. (c) Pixel : 256 ⇥ 256 (d) Blob : ⇠ 3842 blobs Reconstruction using large pixel is stable. 33 / 79
33. Introduction Blob Simulations Real data Perspectives Properties Gaussian blob Tight frame Pixel vs. Gaussian blob Zoomed view on pixel and Gaussian blob reconstructions. (e) Pixel : 256 ⇥ 256 (f) Blob : ⇠ 3842 blobs Large pixel : poor spatial resolution and zigzag artifacts. Blob : lower bandwidth but better image quality. 34 / 79
34. Introduction Blob Simulations Real data Perspectives Properties Gaussian blob Tight frame Section 2 : Blob Basis for CS and CT 1 Literature review and Compressed Sensing theory Background Regularization by TV minimization Regularization by Besov minimization Compressed Sensing and CT 2 Blob Basis for CS and CT General properties of blob Single scale Gaussian blob image Tight frame by a multiscale blob system Multiscale Mexican hat blob 3 2D Numerical Experiments - Simulated data TV reconstruction results `1 reconstruction results TV-`1 reconstruction results 4 2D Numerical Experiments - Real data 5 Conclusion and Persperctives 35 / 79
35. Introduction Blob Simulations Real data Perspectives Properties Gaussian blob Tight frame Shift invariant space generated by a single Gaussian blob Proposition (Riesz basis by Gaussian blob) The translation of Gaussian blob { (x xk)} on an hexagonal lattice is a Riesz basis of the spanned shift invariant space X. Stability of blob representation : { (x xk)} is a Riesz basis of X There exists 0 < A  B < 1 : A X k |fk|2  kfk2  B X k |fk|2, 8f 2 X Unique representation coe cient {fk}k2Z Zd for f 2 X Small pertubation in {fk}k2Z Zd , small pertubation in f Numerically { (x xk)} can be manipulated as an orthogonal basis. 36 / 79
36. Introduction Blob Simulations Real data Perspectives Properties Gaussian blob Tight frame Density of X in L2(IRd) Higher resolution by smaller basis and denser lattice Lattice generated by matrix hG : Lh Dilated basis function : h(x) = h d/2 (x/h) Xh = Span{ h(x yk), k 2 Z Zd}, with yk = hGk Question : Is Xh dense in L2(IRd) as h ! 0 ? Admissible basis function 1 Sampling kernel : ˆ(0) 6= 0 2 { (x xk)}k is a Riesz basis of X 3 R k!k1>R |ˆ(!)|2d! = O(R p), for some p > 0 4 |ˆ(x + y) ˆ(x)|  Ckxk p0 kyk, 8y 2 IRd, for some p0 > d 37 / 79
37. Introduction Blob Simulations Real data Perspectives Properties Gaussian blob Tight frame Density of X in L2(IRd) Higher resolution by smaller basis and denser lattice Lattice generated by matrix hG : Lh Dilated basis function : h(x) = h d/2 (x/h) Xh = Span{ h(x yk), k 2 Z Zd}, with yk = hGk Proposition (Density of Xh in L2(IRd)) As h ! 0, Xh generated by an admissible basis function is dense in L2(IRd) if and only if : |ˆ(G >k)| = 0, 8k 2 Z Zd, k 6= 0 (1) That is, ˆ(!) vanishes on all non zero nodes of the dual lattice. There exists L2(IRd) functions which cannot be approximated by Gaussian blob. 38 / 79
38. Introduction Blob Simulations Real data Perspectives Properties Gaussian blob Tight frame Reconstruction by TV minimization Reconstruction by TV minimization min f2I RN 1 2 kAf bk2 + µkfkT V , or min f2I RN kfkT V , s.t. kAf bk2  "2 Positive constraint fn 0 can be incorporated µ or " is set manually according to the noise level Solved by TVAL3 algorithm (augmented lagrangien) Discrete TV norm of a blob image f Approximation of R I Rd krf(x)kdx on L0 using trapezoidal rule : kfkT V 4 = X k2Z Zd q (D1f)2 k + · · · + (Ddf)2 k Di : X ! IR|L0| the i-th directional derivative of f taken on L0 : (Dif)k = @if(x0 k ), x0 k 2 L0 Implementation of Di Two possibilities : 1 Direct parallelization on GPU platform 2 FFT technique : d @if(!) = 2⇡i! ˆ f(!) Parallelization on GPU is more e cient and avoids aliasing errors. 39 / 79
39. Introduction Blob Simulations Real data Perspectives Properties Gaussian blob Tight frame Section 2 : Blob Basis for CS and CT 1 Literature review and Compressed Sensing theory Background Regularization by TV minimization Regularization by Besov minimization Compressed Sensing and CT 2 Blob Basis for CS and CT General properties of blob Single scale Gaussian blob image Tight frame by a multiscale blob system Multiscale Mexican hat blob 3 2D Numerical Experiments - Simulated data TV reconstruction results `1 reconstruction results TV-`1 reconstruction results 4 2D Numerical Experiments - Real data 5 Conclusion and Persperctives 40 / 79
40. Introduction Blob Simulations Real data Perspectives Properties Gaussian blob Tight frame Partition of unity and sampling condition for a TF Blob family A family { j(x)}j2Z Z of blobs constituted by band-pass functions of di↵erent frequency selectivities. Partition of unity (PU) X j | ˆ j(!)|2 = X j b ˆ j (k!k)2 = c0 > 0, for a.e. ! 2 IRd Sampling condition for a tight frame (TF) j are band-limited functions satisfying PU with c0 = 1 Lj : hexagonal lattice generated by matrix Gj supp ˆ j ⇢ j a fundamental domain of the dual lattice L⇤ j j,k(x) 4 = |det Gj|1/2 j(x xj k ), with xj k = Gjk 41 / 79
41. Introduction Blob Simulations Real data Perspectives Properties Gaussian blob Tight frame Tight frame of a multiscale blob system Proposition (Tight frame) The family { j,k(x)}j2Z Z,k2Z Zd constitutes a tight frame of L2(IRd) : f(x) = X j2Z Z X k2Z Zd hf, j,ki | {z } analysis coe cient j,k, 8f 2 L2(IRd) Any f 2 L2(IRd) can be synthesized by f(x) = X j,k fj,k |{z} synthesis coe cient j,k(x), with {fj,k} 2 `2(Z Z) Redundant system : {fj,k} may be di↵erent from hf, j,ki Construct { j,k} such that hf, j,ki has fast decay Find a sparse solution {fj,k} through `1 minimization {hf, j,ki} as a theoretical guarantee on the decaying rate of {fj,k} 42 / 79
42. Introduction Blob Simulations Real data Perspectives Properties Gaussian blob Tight frame Construction of multiscale blob system Step 1 : mother blob (x) and mother hexagonal lattice Choose the mother blob (x) of n vanishing moments. For positive 0  k1 + · · · + kd < n : Z I Rd (x)xkdx = 0, with xk 4 = xk1 1 · · · xkd d Construct mother hexagonal lattice L0 by the matrix G0, such that ˆ is supported in a fundamental domain of L⇤ 0 . Step 2 : multiscale band-pass blobs and hexagonal lattices Define the blob family j(x) = jd ( jx), > 1 fulfilling PU. Construct the hexagonal lattice Lj by matrix Gj = jG0 Translatation and normalization : j,k(x) = |det Gj|1/2 j(x xj k ). { j,k(x)} is a tight frame of L2(IRd). hf, j,ki decays as O( nj) inside polynomial regions Significant coe cients are around the discontinuities of f N-th largest coe cient |hf, j,ki| decays as Nd/(2 2d) 43 / 79
43. Introduction Blob Simulations Real data Perspectives Properties Gaussian blob Tight frame Mexcian hat blob Profile b ˆ in the frequency domain Define the band-pass profile function in frequency domain : b ˆ (s) = c0s2 exp( c1s2), c0 > 0, c1 > 0 Proposition (Relation between b and b ˆ by Abel transform) A · · · A | {z } d 1 b ˆ (s) = c b (s), s 2 IR 2D Mexican hat blob We have the following Fourier transform pairs in IR2 : (x) = (1 ↵kxk2) exp( ↵kxk2) ˆ (!) = ⇡3↵ 2k!k2 exp( ⇡2↵ 1k!k2) Mexican hat blob has two vanishing moments. 44 / 79
44. Introduction Blob Simulations Real data Perspectives Properties Gaussian blob Tight frame Mexican hat blob - Visualization Blob family and PU Partition of unity holds for  p 2. Multiscale blob family : j(x) = jd ( jx) 2D mother Mexican hat blob with the parameter ↵ = 1 : (a) 2D Mexican hat blob profile (b) 2D Mexican hat blob 45 / 79
45. Introduction Blob Simulations Real data Perspectives Properties Gaussian blob Tight frame Reconstruction by `1 minimization of blob coe cients Reconstruction by `1 minimization min f2I RN 1 2 kAf bk2 + µkfk1, or min f2I RN kfk1, s.t. kAf bk2  "2 with kfk1 4 = P n |fn| A operates directly on blob coe↵cients f : no needs for sparsifying transform. µ or " is set manually according to noise level Solved by FISTA and ADM (Alternating Direction Method) Modify the multiscale blob system for numerical implementation Only a finite number of scales is used in practice : f(x) = X k2Z Zd f0,k 0(x x0 k ) | {z } low-pass blob + 1 X j= S X k2Z Zd fj,k j(x xj k ) | {z } band-pass blob For Mexican hat blob system 0 is a Gaussian blob The modified system is no more a tight frame 46 / 79
46. Introduction Blob Simulations Real data Perspectives TV `1 TV-`1 Outline 1 Literature review and Compressed Sensing theory Background Regularization by TV minimization Regularization by Besov minimization Compressed Sensing and CT 2 Blob Basis for CS and CT General properties of blob Single scale Gaussian blob image Tight frame by a multiscale blob system Multiscale Mexican hat blob 3 2D Numerical Experiments - Simulated data TV reconstruction results `1 reconstruction results TV-`1 reconstruction results 4 2D Numerical Experiments - Real data 5 Conclusion and Persperctives 47 / 79
47. Introduction Blob Simulations Real data Perspectives TV `1 TV-`1 Sinogram simulation and image quality assessment Sinogram simulation Number of projections : P, number of detector pixels : 512 Apply the parallel beam projector A on the pixel image g : y = Ag Sample the Poisson distribution Ii ⇠ Poisson(I0 exp( yi)) Noisy sinogram b : bi = log I0 log Ii. SNR of b : SNR(b) 4 = 10 log10 ✓ Ekyk2 Eky bk2 ◆ Image quality assessment The blob image f is resampled to the same dimension of g. SNR : SNR(f, g) 4 = 20 log10 ✓ kg N 1 P i gik kh gk ◆ Streak Index (SI, smaller is better) : SI(f, g) 4 = N 1kf gkT V 48 / 79
48. Introduction Blob Simulations Real data Perspectives TV `1 TV-`1 Section 3 : 2D Numerical Experiments - Simulated data 1 Literature review and Compressed Sensing theory Background Regularization by TV minimization Regularization by Besov minimization Compressed Sensing and CT 2 Blob Basis for CS and CT General properties of blob Single scale Gaussian blob image Tight frame by a multiscale blob system Multiscale Mexican hat blob 3 2D Numerical Experiments - Simulated data TV reconstruction results `1 reconstruction results TV-`1 reconstruction results 4 2D Numerical Experiments - Real data 5 Conclusion and Persperctives 49 / 79
49. Introduction Blob Simulations Real data Perspectives TV `1 TV-`1 Reconstruction quality vs. lattice sampling step h P = 64 projections, SNR(b)=50dB Choose the sampling step h such that N ⇠ 128 ⇥ 84 The same TV minimization problem is solved Blob : SNR=16.28dB, SI=0.0217 Pixel : SNR=13.74dB, SI=0.0348 (a) Blob, h=1.460, SNR=16.28dB, SI=0.0217, T =63.14s (b) Pixel, h=1.309, SNR=13.74dB, SI=0.0348, T =12.47s 50 / 79
50. Introduction Blob Simulations Real data Perspectives TV `1 TV-`1 Reconstruction quality vs. lattice sampling step h P = 64 projections, SNR(b)=50dB Choose the sampling step h such that N ⇠ 256 ⇥ 168 The same TV minimization problem is solved Blob : SNR=18.98dB, SI=0.017 Pixel : SNR=15.71dB, SI=0.0289 (c) Blob, h=0.730, SNR=18.98dB, SI=0.017, T =98.91s (d) Pixel, h=0.654, SNR=15.71dB, SI=0.0289, T =51.89s 51 / 79
51. Introduction Blob Simulations Real data Perspectives TV `1 TV-`1 Reconstruction quality vs. lattice sampling step h (a) SNR (b) SI Blob is less sensitive to the variation of h(N). Blob achieves the same visual quality with less unknowns. Computation time for blob is close to pixel. Similar results hold for other medical oriented phantoms. 52 / 79
52. Introduction Blob Simulations Real data Perspectives TV `1 TV-`1 Reconstruction quality vs. lattice sampling step h (a) SNR (b) Time Blob is less sensitive to the variation of h(N). Blob achieves the same visual quality with less unknowns. Computation time for blob is close to pixel in middle and high dimension. Similar results hold for other medical oriented phantoms. 53 / 79
53. Introduction Blob Simulations Real data Perspectives TV `1 TV-`1 Reconstruction from few projections P = 64, 32 projections (a) TV, P =64, SNR=19.78dB, SI=0.016 (b) TV, P =32, SNR=17.65dB, SI=0.019 TV minimization removes e ciently the streak artifacts. ”Cartoon” artifacts are amplified when P is small. 54 / 79
54. Introduction Blob Simulations Real data Perspectives TV `1 TV-`1 Reconstruction from few projections TV reconstruction quality vs. number of projections on di↵erent phantoms : (a) Shepp-Logan (b) Abdomen (c) Thorax 55 / 79
55. Introduction Blob Simulations Real data Perspectives TV `1 TV-`1 Reconstruction from few projections TV reconstruction quality vs number of projections on di↵erent phantoms : (a) SNR (b) SI Sharp transitions in the SNR plot Monotone evolution for the Shepp-Logan and the Abdomen phantoms Evolution stagnated when P becomes large 56 / 79
56. Introduction Blob Simulations Real data Perspectives TV `1 TV-`1 Section 3 : 2D Numerical Experiments - Simulated data 1 Literature review and Compressed Sensing theory Background Regularization by TV minimization Regularization by Besov minimization Compressed Sensing and CT 2 Blob Basis for CS and CT General properties of blob Single scale Gaussian blob image Tight frame by a multiscale blob system Multiscale Mexican hat blob 3 2D Numerical Experiments - Simulated data TV reconstruction results `1 reconstruction results TV-`1 reconstruction results 4 2D Numerical Experiments - Real data 5 Conclusion and Persperctives 57 / 79
57. Introduction Blob Simulations Real data Perspectives TV `1 TV-`1 Sparsity of the `1 reconstruction - scale decomposition Mex-4, P = 128 projections, SNR(b)=50db, Lung phantom, N = 289243. (a) Scale 0, ⌧=84.7% (b) Scale -1, ⌧=20.7% (c) Scale -2, ⌧=16.3% (d) Scale -3, ⌧=4.9% 58 / 79
58. Introduction Blob Simulations Real data Perspectives TV `1 TV-`1 Sparsity of the `1 reconstruction Mex-4, P = 128 projections, SNR(b)=50db, Lung phantom, N = 289243. (a) `1, ⌧=8.9%, SNR=20.36dB, SI=0.0159 (b) TV, ⌧=100%, SNR=20.60dB, SI=0.0149 `1 and TV minimizations give similar image qualities for large P. TV minimization doesn’t enhance the sparsity of blob coe cients. 59 / 79
59. Introduction Blob Simulations Real data Perspectives TV `1 TV-`1 Sparsity of the `1 reconstruction : zoomed view Mex-4, P = 256 projections, SNR(b)=50db, Abdomen phantom (a) Phantom (b) `1, P =256 (c) TV, P =256 Oscillating patterns and low contrast regions Always removed by TV minimization (”texture killer”), even for large P. Well preserved by `1 minimization for large P. 60 / 79
60. Introduction Blob Simulations Real data Perspectives TV `1 TV-`1 TV vs `1 minimization Abdomen and Thorax phantoms : P 2 [32, 256], SNR(b)=50dB. (a) SNR (b) SI `1 minimization compared with TV Substantially better than TV for large P No ”texture killer” behavior Monotonical improvement of the image quality Far less e cient than TV for small P 61 / 79
61. Introduction Blob Simulations Real data Perspectives TV `1 TV-`1 `1 reconstruction using few projections : scale image Mex-4, P =48 and 128, Abdomen phantom, scale -2 (a) P =48, scale -2 (b) P =128, scale -2 The streak artifacts corrupt even coarse scales when P is small. 62 / 79
62. Introduction Blob Simulations Real data Perspectives TV `1 TV-`1 Comparison with wavelet reconstruction `1 reconstruction based on Mex-4 and Daub-6 wavelet, Abdomen phantom, zoomed views on the ROI, P=96, SNR(b)=50dB. (a) Mex-4, SNR=18.53dB, SI=0.0111 (b) Daub-6, SNR=16.40dB, SI=0.0175 (c) Phantom The multiscale blob system is redundant, and works in a lower bandwidth than the wavelet reconstruction does, therefore less concerned by the high frequency pseudo Gibbs artifacts. The blob reconstruction su↵ers from the ”salt and pepper” noise (misplaced and isolated blobs due to energy transfert between scales). 63 / 79
63. Introduction Blob Simulations Real data Perspectives TV `1 TV-`1 Comparison with wavelet reconstruction `1 reconstruction based on Mex-4 and Daub-6 wavelet, Abdomen phantom, zoomed views on the ROI, P=96, SNR(b)=50dB. (a) Mex-4, SNR=18.53dB, SI=0.0111 (b) Daub-6, SNR=16.40dB, SI=0.0175 (c) Phantom The multiscale blob system is redundant, and works in a lower bandwidth than the wavelet reconstruction does, therefore less concerned by the high frequency pseudo Gibbs artifacts. The blob reconstruction su↵ers from the ”salt and pepper” noise (misplaced and isolated blobs due to energy transfert between scales). 64 / 79
64. Introduction Blob Simulations Real data Perspectives TV `1 TV-`1 Section 3 : 2D Numerical Experiments - Simulated data 1 Literature review and Compressed Sensing theory Background Regularization by TV minimization Regularization by Besov minimization Compressed Sensing and CT 2 Blob Basis for CS and CT General properties of blob Single scale Gaussian blob image Tight frame by a multiscale blob system Multiscale Mexican hat blob 3 2D Numerical Experiments - Simulated data TV reconstruction results `1 reconstruction results TV-`1 reconstruction results 4 2D Numerical Experiments - Real data 5 Conclusion and Persperctives 65 / 79
65. Introduction Blob Simulations Real data Perspectives TV `1 TV-`1 Improve the `1 reconstruction by TV-`1 minimization TV-`1 minimization `1 minimization augmented by a TV term (µ1 µ2) : min f2I RN 1 2 kAf bk2 + µ1kfk 1 + µ2kfkT V Abdomen phantom : reconstructions obtained using di↵erent methods, P = 128 projections, data SNR=50dB, µ1/µ2 30. (a) TV, SNR=18.52dB, SI=0.0086 (b) TV-`1, SNR=19.60dB, SI=0.0085 (c) `1, SNR=19.80dB, SI=0.0097 66 / 79
66. Introduction Blob Simulations Real data Perspectives TV `1 TV-`1 Improve the `1 reconstruction by TV-`1 minimization TV-`1 minimization `1 minimization augmented by a TV term (µ1 µ2) : min f2I RN 1 2 kAf bk2 + µ1kfk 1 + µ2kfkT V Abdomen phantom : reconstructions obtained using di↵erent methods, P = 128 projections, data SNR=50dB, µ1/µ2 30. (d) TV, SNR=18.52dB, SI=0.0086 (e) TV-`1, SNR=19.60dB, SI=0.0085 (f) `1, SNR=19.80dB, SI=0.0097 67 / 79
67. Introduction Blob Simulations Real data Perspectives Outline 1 Literature review and Compressed Sensing theory Background Regularization by TV minimization Regularization by Besov minimization Compressed Sensing and CT 2 Blob Basis for CS and CT General properties of blob Single scale Gaussian blob image Tight frame by a multiscale blob system Multiscale Mexican hat blob 3 2D Numerical Experiments - Simulated data TV reconstruction results `1 reconstruction results TV-`1 reconstruction results 4 2D Numerical Experiments - Real data 5 Conclusion and Persperctives 68 / 79
68. Introduction Blob Simulations Real data Perspectives Medical application Dose reduction by few projections Decalcification detection for aged person by X-ray imaging. (a) View from y axis (b) View from z axis Figure: A trunk of human spine fixed in resins. 69 / 79
69. Introduction Blob Simulations Real data Perspectives Medical application : scale decomposition TV-`1 reconstruction, Mex-4, using partial data : P = 160 projections in [0, 2⇡) and 512 detector bins Number of blobs N = 202522 Scale 0, -1 : (a) Scale 0, ⌧ =90.28% (b) Scale -1, ⌧ =35.41% 70 / 79
70. Introduction Blob Simulations Real data Perspectives Medical application : scale decomposition TV-`1 reconstruction, Mex-4, using partial data : P = 160 projections in [0, 2⇡) and 512 detector bins Number of blobs N = 202522 Scale -2, -3 : (c) Scale -2, ⌧ =28.77% (d) Scale -3, ⌧ =11.34% 71 / 79
71. Introduction Blob Simulations Real data Perspectives Medical application : reconstructions TV-`1 reconstruction, Mex-4, using P = 160 projections and 512 detector bins Tikhonov regularization, pixel, using P = 360 projections and 1024 detector bins (a) TV-`1, P =160, ⌧=16.40% (b) Tikhonov, pixel, P =360, 100 CG iterations 72 / 79
72. Introduction Blob Simulations Real data Perspectives Medical application : comparison Zoomed views on the ROI : (a) Tikhonov, Pixel, P =360 (b) TV, GS, P =160 (c) TV-`1, Mex-4, P =160 Number of projections can be reduced by half Better image quality and spatial resolution than classical methods Similar results have been observed in Micro-CT and Electron Tomography applications. 73 / 79
73. Introduction Blob Simulations Real data Perspectives Medical application : comparison Zoomed views on the ROI : (a) Tikhonov, Pixel, P =360 (b) TV, GS, P =160 (c) TV-`1, Mex-4, P =160 Number of projections can be reduced by half Better image quality and spatial resolution than classical methods Similar results have been observed in Micro-CT and Electron Tomography applications. 74 / 79
74. Introduction Blob Simulations Real data Perspectives Outline 1 Literature review and Compressed Sensing theory Background Regularization by TV minimization Regularization by Besov minimization Compressed Sensing and CT 2 Blob Basis for CS and CT General properties of blob Single scale Gaussian blob image Tight frame by a multiscale blob system Multiscale Mexican hat blob 3 2D Numerical Experiments - Simulated data TV reconstruction results `1 reconstruction results TV-`1 reconstruction results 4 2D Numerical Experiments - Real data 5 Conclusion and Persperctives 75 / 79
75. Introduction Blob Simulations Real data Perspectives Conclusions CS framework for few projections problem Sparse modeling/representation of functions Reconstruction via TV/`1 minimization Main contributions Image representation models by blob Adapdation of TV/`1 minimizations Advantages of blob over pixel More stable reconstruction by working in lower bandwidth Equivalent or better reconstruction quality on medical oriented phantoms Easy parallelization and no needs for sparsifying transforms Choice of the image model and the reconstruction method Blocky objects and small P : Gaussian blob + TV minimization Medical objects and large P : Multiscale blob system + `1 minimization P in a middle range : Multiscale blob system + TV-`1 minimization A large application field can be covered. 76 / 79
76. Introduction Blob Simulations Real data Perspectives Perspectives 3D generalization Our theories and algorithms have been developed for d-dimensional case, and applied directly for 3D case. High dimension in 3D : special cares are needed for the parallelization on GPU. Contour blob Multiscale blob system is a mimic of the Wavelet basis Slow decaying rate of N-th largest coe↵cient in 2D : ⇠ N 1/p, with p ⇠ 1.6 Radial symmetry of blob is a limitation 77 / 79
77. Introduction Blob Simulations Real data Perspectives Contour blob Contour blob by tensorial product Gaussian : g(t) = exp( ↵1t2), Mexican hat : m(t) = exp( ↵2t2)(1 2↵2t2) (x) 4 = g(x1)m(x2) Properties Geometrical sensitivity : the coe cients of blob will be cancelled out unless they are “perfectly” aligned with the contour CT-friendly as radial blob : analytical expression for the X-ray transform 78 / 79
78. Introduction Blob Simulations Real data Perspectives Thank you for your attention ! 79 / 79