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
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
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
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
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
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
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
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
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
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
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
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 ampliﬁed for small number of projection. 16 / 79
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
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
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
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 satisﬁes 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
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
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
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
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
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
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
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 Inﬁnite 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
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
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
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
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
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
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
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
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
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
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
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
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
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 Deﬁne the blob family j(x) = jd ( jx), > 1 fulﬁlling 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 Signiﬁcant coe cients are around the discontinuities of f N-th largest coe cient |hf, j,ki| decays as Nd/(2 2d) 43 / 79
frame Mexcian hat blob Proﬁle b ˆ in the frequency domain Deﬁne the band-pass proﬁle 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
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 proﬁle (b) 2D Mexican hat blob 45 / 79
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 ﬁnite 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 modiﬁed system is no more a tight frame 46 / 79
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
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
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
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
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
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
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
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 ampliﬁed when P is small. 54 / 79
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
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
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
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
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
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
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
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
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
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
by few projections Decalciﬁcation detection for aged person by X-ray imaging. (a) View from y axis (b) View from z axis Figure: A trunk of human spine ﬁxed in resins. 69 / 79
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
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
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
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
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 ﬁeld can be covered. 76 / 79
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
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