A Sparse Tour of Imaging Sciences

Talk at Dassault System

Transcript

1. A Sparse Tour of Imaging Sciences Gabriel Peyré www.numerical-tours.com
2. Signals, Images and More
3. Signals, Images and More
4. Signals, Images and More
5. Signals, Images and More
6. Signals, Images and More
7. Overview • Approximation in an Ortho-Basis • Compression and Denoising • Sparse Inverse Problem Regularization • Compressed Sensing • Iterative Soft Thresholding Algorithm
8. Orthogonal basis { m }m of L2([0, 1]d) Continuous signal/image f L2([0, 1]d). Orthogonal Decompositions
9. Orthogonal basis { m }m of L2([0, 1]d) f = m f, m m ||f|| = |f(x)|2dx = m | f, m ⇥|2 Continuous signal/image f L2([0, 1]d). Orthogonal Decompositions
10. Orthogonal basis { m }m of L2([0, 1]d) f = m f, m m ||f|| = |f(x)|2dx = m | f, m ⇥|2 Continuous signal/image f L2([0, 1]d). Orthogonal Decompositions m
11. 1-D Wavelet Basis Wavelets: j,n (x) = 1 2j/2 x 2jn 2j Position n, scale 2j, m = (n, j).
12. 1-D Wavelet Basis Wavelets: j,n (x) = 1 2j/2 x 2jn 2j Position n, scale 2j, m = (n, j).
13. m1,m2 Basis { m1,m2 (x1, x2 )}m1,m2 of L2([0, 1]2) m1,m2 (x1, x2 ) = m1 (x1 ) m2 (x2 ) tensor product 2-D Fourier Basis Basis { m (x)}m of L2([0, 1]) m1 m2
14. m1,m2 Basis { m1,m2 (x1, x2 )}m1,m2 of L2([0, 1]2) m1,m2 (x1, x2 ) = m1 (x1 ) m2 (x2 ) tensor product f(x) f, m1,m2 Fourier transform 2-D Fourier Basis Basis { m (x)}m of L2([0, 1]) m1 m2 x m
15. 3 elementary wavelets { H, V , D}. Orthogonal basis of L2([0, 1]2): k j,n (x) = 2 j (2 jx n) k=H,V,D j<0,2j n [0,1]2 2-D Wavelet Basis V (x) H(x) D(x)
16. 3 elementary wavelets { H, V , D}. Orthogonal basis of L2([0, 1]2): k j,n (x) = 2 j (2 jx n) k=H,V,D j<0,2j n [0,1]2 2-D Wavelet Basis V (x) H(x) D(x)
17. wavelet f, k j,n Example of Wavelet Decomposition f(x) transform x (j, n, k)
18. Discrete Computations Discrete orthogonal basis { m } of CN . f = m f, m m
19. Fast Fourier Transform (FFT), O(N log(N)) operations. Discrete Computations Discrete orthogonal basis { m } of CN . m [n] = 1 N e2i N nm f = m f, m m
20. Fast Fourier Transform (FFT), O(N log(N)) operations. Fast Wavelet Transform, O(N) operations. Discrete Wavelet basis: no closed-form expression. Discrete Computations Discrete orthogonal basis { m } of CN . m [n] = 1 N e2i N nm f = m f, m m
21. Sparse Approximation in a Basis
22. Sparse Approximation in a Basis
23. Sparse Approximation in a Basis
24. Best basis Fastest error decay ||f fM ||2 log(||f fM ||) log(M) Efficiency of Transforms Fourier DCT Local DCT Wavelets
25. Overview • Approximation in an Ortho-Basis • Compression and Denoising • Sparse Inverse Problem Regularization • Compressed Sensing • Iterative Soft Thresholding Algorithm
26. JPEG-2000 vs. JPEG, 0.2bit/pixel
27. Compression by Transform-coding Image f Zoom on f f forward a[m] = ⇥f, m ⇤ R transform
28. Compression by Transform-coding Image f Zoom on f f forward a[m] = ⇥f, m ⇤ R transform Quantization: q[m] = sign(a[m]) |a[m]| T ⇥ Z ˜ a[m] T T 2T 2T a[m] Quantized q[m] bin T q[m] Z
29. Compression by Transform-coding Image f Zoom on f f forward a[m] = ⇥f, m ⇤ R coding transform Entropic coding: use statistical redundancy (many 0’s). Quantization: q[m] = sign(a[m]) |a[m]| T ⇥ Z ˜ a[m] T T 2T 2T a[m] Quantized q[m] bin T q[m] Z
30. Compression by Transform-coding Image f Zoom on f f forward a[m] = ⇥f, m ⇤ R coding decoding q[m] Z transform Entropic coding: use statistical redundancy (many 0’s). Quantization: q[m] = sign(a[m]) |a[m]| T ⇥ Z ˜ a[m] T T 2T 2T a[m] Quantized q[m] bin T q[m] Z
31. Compression by Transform-coding Image f Zoom on f f forward Dequantization: ˜ a[m] = sign(q[m]) |q[m] + 1 2 ⇥ T a[m] = ⇥f, m ⇤ R coding decoding q[m] Z ˜ a[m] dequantization transform Entropic coding: use statistical redundancy (many 0’s). Quantization: q[m] = sign(a[m]) |a[m]| T ⇥ Z ˜ a[m] T T 2T 2T a[m] Quantized q[m] bin T q[m] Z
32. Compression by Transform-coding Image f Zoom on f f , R =0.2 bit/pixel f forward Dequantization: ˜ a[m] = sign(q[m]) |q[m] + 1 2 ⇥ T a[m] = ⇥f, m ⇤ R coding decoding q[m] Z ˜ a[m] dequantization transform backward fR = m IT ˜ a[m] m transform Entropic coding: use statistical redundancy (many 0’s). Quantization: q[m] = sign(a[m]) |a[m]| T ⇥ Z ˜ a[m] T T 2T 2T a[m] Quantized q[m] bin T q[m] Z
33. Compression by Transform-coding Image f Zoom on f f , R =0.2 bit/pixel f forward Dequantization: ˜ a[m] = sign(q[m]) |q[m] + 1 2 ⇥ T a[m] = ⇥f, m ⇤ R coding decoding q[m] Z ˜ a[m] dequantization transform backward fR = m IT ˜ a[m] m transform Entropic coding: use statistical redundancy (many 0’s). Quantization: q[m] = sign(a[m]) |a[m]| T ⇥ Z ˜ a[m] T T 2T 2T a[m] Quantized q[m] bin T q[m] Z ||f fM ||2 = O(M ) =⇥ ||f fR ||2 = O(log (R)R ) Theorem:
34. Noise in Images
35. Denoising
36. Denoising thresh. f = N 1 m=0 f, m ⇥ m ˜ f = | f, m ⇥|>T f, m ⇥ m
37. Denoising thresh. f = N 1 m=0 f, m ⇥ m ˜ f = | f, m ⇥|>T f, m ⇥ m In practice: T 3 for T = 2 log(N) Theorem: if ||f0 f0,M ||2 = O(M ), E(|| ˜ f f0 ||2) = O( 2 +1 )
38. Overview • Approximation in an Ortho-Basis • Compression and Denoising • Sparse Inverse Problems Regularization • Compressed Sensing • Iterative Soft Thresholding Algorithm
39. Recovering f0 2 RN from noisy observations: y = Kf0 + w 2 RP : RN ! RP with P ⌧ N (missing information) Inverse Problems
40. Examples: Inpainting, super-resolution, . . . Recovering f0 2 RN from noisy observations: y = Kf0 + w 2 RP : RN ! RP with P ⌧ N (missing information) f0 K Inverse Problems K
41. Kf = (p✓k )k Inverse Problems in Medical Imaging
42. Magnetic resonance imaging (MRI): Kf = (p✓k )k Kf = ( ˆ f(!))!2⌦ Inverse Problems in Medical Imaging ˆ x
43. Magnetic resonance imaging (MRI): Other examples: MEG, EEG, . . . Kf = (p✓k )k Kf = ( ˆ f(!))!2⌦ Inverse Problems in Medical Imaging ˆ x
44. L1 Regularization coe cients x0 RN
45. L1 Regularization coe cients image x0 RN f0 = x0 RQ
46. L1 Regularization observations w coe cients image K x0 RN f0 = x0 RQ y = Kf0 + w RP
47. y = K f0 + w = x0 + w 2 RP Equivalent model: L1 Regularization observations = K ⇥ ⇥ RP N w coe cients image K x0 RN f0 = x0 RQ y = Kf0 + w RP
48. y = K f0 + w = x0 + w 2 RP Equivalent model: If is invertible: 1 y = x0 + 1 w L1 Regularization observations = K ⇥ ⇥ RP N w coe cients image K x0 RN f0 = x0 RQ y = Kf0 + w RP
49. y = K f0 + w = x0 + w 2 RP Equivalent model: If is invertible: 1 y = x0 + 1 w L1 Regularization observations = K ⇥ ⇥ RP N w coe cients image K x0 RN f0 = x0 RQ y = Kf0 + w RP Problems: 1w can “explose”. can even be non-invertible.
50. Inverse Problem Regularization observations y parameter Estimator: x(y) depends only on Observations: y = x0 + w 2 RP .
51. Inverse Problem Regularization observations y parameter Example: variational methods Estimator: x(y) depends only on x ( y ) 2 argmin x 2RN 1 2 || y x ||2 + J ( x ) Data fidelity Regularity Observations: y = x0 + w 2 RP .
52. J ( x0) Regularity of x0 Inverse Problem Regularization observations y parameter Example: variational methods Estimator: x(y) depends only on x ( y ) 2 argmin x 2RN 1 2 || y x ||2 + J ( x ) Data fidelity Regularity Observations: y = x0 + w 2 RP . Choice of : tradeo ||w|| Noise level
53. J ( x0) Regularity of x0 x ( y ) 2 argmin x = y J ( x ) Inverse Problem Regularization observations y parameter Example: variational methods Estimator: x(y) depends only on x ( y ) 2 argmin x 2RN 1 2 || y x ||2 + J ( x ) Data fidelity Regularity Observations: y = x0 + w 2 RP . No noise: 0+, minimize Choice of : tradeo ||w|| Noise level
54. J0( x ) = # { m ; xm 6= 0} “Ideal” sparsity: Sparse Priors
55. Sparse regularization: x? 2 argmin x 1 2 || y x ||2 + J0( x ) J0( x ) = # { m ; xm 6= 0} “Ideal” sparsity: Sparse Priors
56. Sparse regularization: Denoising in ortho-basis: K = Id, 2 = Id, ⇤ = Id x? 2 argmin x 1 2 || y x ||2 + J0( x ) J0( x ) = # { m ; xm 6= 0} “Ideal” sparsity: Sparse Priors y
57. Sparse regularization: Denoising in ortho-basis: K = Id, 2 = Id, ⇤ = Id x? 2 argmin x 1 2 || y x ||2 + J0( x ) J0( x ) = # { m ; xm 6= 0} “Ideal” sparsity: where ˜ x = ⇤ y = {h y, m i}m, T 2 = 2 min x || ⇤ y x ||2 + T 2 J0( x ) = P m |˜ xm xm |2 + T 2 ( xm ) Sparse Priors y ?y
58. Sparse regularization: Denoising in ortho-basis: K = Id, 2 = Id, ⇤ = Id x? 2 argmin x 1 2 || y x ||2 + J0( x ) J0( x ) = # { m ; xm 6= 0} “Ideal” sparsity: where ˜ x = ⇤ y = {h y, m i}m, T 2 = 2 min x || ⇤ y x ||2 + T 2 J0( x ) = P m |˜ xm xm |2 + T 2 ( xm ) Solution: x ? m = ⇢ xm if | ˜ xm | > T, 0 otherwise. Sparse Priors y ?y ? x ? x ?
59. Sparse regularization: Denoising in ortho-basis: K = Id, 2 = Id, ⇤ = Id Non-orthogonal : NP-hard to solve. x? 2 argmin x 1 2 || y x ||2 + J0( x ) J0( x ) = # { m ; xm 6= 0} “Ideal” sparsity: where ˜ x = ⇤ y = {h y, m i}m, T 2 = 2 min x || ⇤ y x ||2 + T 2 J0( x ) = P m |˜ xm xm |2 + T 2 ( xm ) Solution: x ? m = ⇢ xm if | ˜ xm | > T, 0 otherwise. Sparse Priors y ?y ? x ? x ?
60. Image with 2 pixels: q = 0 J0 (x) = # {m \ xm = 0} J0 (x) = 0 null image. J0 (x) = 1 sparse image. J0 (x) = 2 non-sparse image. x2 Convex Relaxation: L1 Prior x1
61. Image with 2 pixels: q = 0 q = 1 q = 2 q = 3/2 q = 1/2 Jq (x) = m |xm |q J0 (x) = # {m \ xm = 0} J0 (x) = 0 null image. J0 (x) = 1 sparse image. J0 (x) = 2 non-sparse image. x2 Convex Relaxation: L1 Prior q priors: (convex for q 1) x1
62. Image with 2 pixels: q = 0 q = 1 q = 2 q = 3/2 q = 1/2 Jq (x) = m |xm |q J0 (x) = # {m \ xm = 0} J0 (x) = 0 null image. J0 (x) = 1 sparse image. J0 (x) = 2 non-sparse image. x2 J1 (x) = m |xm | Convex Relaxation: L1 Prior Sparse 1 prior: q priors: (convex for q 1) x1
63. x argmin x=y m |xm | x x = y Noiseless Sparse Regularization Noiseless measurements: y = x0
64. x argmin x=y m |xm | x x = y x argmin x=y m |xm |2 Noiseless Sparse Regularization x x = y Noiseless measurements: y = x0
65. x argmin x=y m |xm | x x = y x argmin x=y m |xm |2 Noiseless Sparse Regularization Convex linear program. Interior points, cf. [Chen, Donoho, Saunders] “basis pursuit”. Douglas-Rachford splitting, see [Combettes, Pesquet]. x x = y Noiseless measurements: y = x0
66. Regularization Data fidelity y = x0 + w Noisy measurements: x argmin x RQ 1 2 ||y x||2 + ||x||1 Noisy Sparse Regularization
67. Regularization Data fidelity Equivalence || x = y|| y = x0 + w Noisy measurements: x argmin x RQ 1 2 ||y x||2 + ||x||1 x argmin || x y|| ||x||1 Noisy Sparse Regularization x
68. Iterative soft thresholding Forward-backward splitting Algorithms: Regularization Data fidelity Equivalence || x = y|| y = x0 + w Noisy measurements: x argmin x RQ 1 2 ||y x||2 + ||x||1 x argmin || x y|| ||x||1 Noisy Sparse Regularization Nesterov multi-steps schemes. see [Daubechies et al], [Pesquet et al], etc x
69. K y = Kf0 + w Measures: (Kf)i = ⇢ 0 if i 2 ⌦, fi if i / 2 ⌦. Inpainting Problem
70. Overview • Approximation in an Ortho-Basis • Compression and Denoising • Sparse Inverse Problem Regularization • Compressed Sensing • Iterative Soft Thresholding Algorithm
71. f[n] f0 (n/N) Sampling: ˜ f L2([0, 1]d) f RN Idealization: acquisition device Discretization
72. Data aquisition: Sensors Pointwise Sampling and Smoothness ˜ f L2 f RN f[i] = ˜ f(i/N)
73. Data aquisition: Sensors ˜ f(t) = i f[i]h(Nt i) Shannon interpolation: if Supp( ˆ ˜ f) [ N , N ] h(t) = sin( t) t Pointwise Sampling and Smoothness ˜ f L2 f RN f[i] = ˜ f(i/N)
74. Data aquisition: Sensors ˜ f(t) = i f[i]h(Nt i) Natural images are not smooth. Shannon interpolation: if Supp( ˆ ˜ f) [ N , N ] h(t) = sin( t) t Pointwise Sampling and Smoothness ˜ f L2 f RN f[i] = ˜ f(i/N)
75. Data aquisition: Sensors ˜ f(t) = i f[i]h(Nt i) Natural images are not smooth. But can be compressed e ciently. Shannon interpolation: if Supp( ˆ ˜ f) [ N , N ] 0,1,0,. . . h(t) = sin( t) t Sample and compress simultaneously? Pointwise Sampling and Smoothness ˜ f L2 f RN f[i] = ˜ f(i/N) JPEG-2k
76. ˜ f Single Pixel Camera (Rice)
77. ˜ f P measures N micro-mirrors Single Pixel Camera (Rice) y[i] = f, i
78. ˜ f P/N = 0.16 P/N = 0.02 P/N = 1 P measures N micro-mirrors Single Pixel Camera (Rice) y[i] = f, i
79. Physical hardware resolution limit: target resolution f RN . ˜ f L2 f RN y RP micro mirrors array resolution CS hardware K CS Hardware Model CS is about designing hardware: input signals ˜ f L2(R2).
80. Physical hardware resolution limit: target resolution f RN . ˜ f L2 f RN y RP micro mirrors array resolution CS hardware , ... K CS Hardware Model CS is about designing hardware: input signals ˜ f L2(R2). , ,
81. Physical hardware resolution limit: target resolution f RN . ˜ f L2 f RN y RP micro mirrors array resolution CS hardware , ... f Operator K K CS Hardware Model CS is about designing hardware: input signals ˜ f L2(R2). , ,
82. Need to solve y = Kf. More unknown than equations. dim(ker(K)) = N P is huge. Inversion and Sparsity f Operator K
83. Need to solve y = Kf. More unknown than equations. dim(ker(K)) = N P is huge. Prior information: f is sparse in a basis { m }m . J (f) = Card {m \ | f, m | > } is small. Inversion and Sparsity f Operator K f, m f
84. 0 reconstruction: Minimize subject to Kf = y CS Reconstruction J0 (f) = Card {m \ f, m = 0}
85. 0 reconstruction: Minimize subject to Kf = y NP-hard to solve. CS Reconstruction J0 (f) = Card {m \ f, m = 0}
86. 0 reconstruction: Minimize subject to Kf = y 1 reconstruction: m | f, m | Polynomial-time algorithms. NP-hard to solve. CS Reconstruction J0 (f) = Card {m \ f, m = 0} Minimize subject to Kf = y y = f f, 1 f, 2
87. Theorem: [Candes, Romberg, Tao, Donoho, 2004] If f is k-sparse, i.e. J0 (f) k If P C log(N/k)k then 1-CS reconstruction is exact. Theoretical Performance Guaranties
88. Theorem: [Candes, Romberg, Tao, Donoho, 2004] If f is k-sparse, i.e. J0 (f) k If P C log(N/k)k then 1-CS reconstruction is exact. Extensions to: noisy observation y = Kf + approximate sparsity f = fk sparse + Theoretical Performance Guaranties
89. Theorem: [Candes, Romberg, Tao, Donoho, 2004] If f is k-sparse, i.e. J0 (f) k If P C log(N/k)k then 1-CS reconstruction is exact. Extensions to: noisy observation y = Kf + approximate sparsity f = fk sparse + Research problem: optimal value of C ? for N/k = 4, C log(N/k) 5. “CS is 5 less e cient than JPEG-2k” Theoretical Performance Guaranties
90. Overview • Approximation in an Ortho-Basis • Compression and Denoising • Sparse Inverse Problems Regularization • Compressed Sensing • Iterative Soft Thresholding Algorithm
91. Orthogonal-basis: ⇤ = IdN , x = ⇤ f . Regularization-based denoising: x? = argmin x 2RN 1 2 || x y ||2 + J ( x ) Denoising: y = x0 + w 2 RN , K = Id. Sparse regularization: J ( x ) = P m | xm |q (where |a|0 = (a)) Sparse Regularization Denoising
92. Orthogonal-basis: ⇤ = IdN , x = ⇤ f . Regularization-based denoising: x? = argmin x 2RN 1 2 || x y ||2 + J ( x ) Denoising: y = x0 + w 2 RN , K = Id. Sparse regularization: J ( x ) = P m | xm |q x ? m = S q T ( xm) (where |a|0 = (a)) Sparse Regularization Denoising
93. Sparse regularization: x? 2 argmin x 2RN E ( x ) = 1 2 || y x ||2 + || x ||1 E ( x, ˜ x ) = E ( x ) 1 2 || ( x ˜ x )||2 + 1 2 ⌧ || x ˜ x ||2 E (· , ˜ x ) x ˜ x S ⌧ (u) ⌧ < 1/|| ⇤ || Surrogate Functionals Surrogate functional: E(·)
94. Sparse regularization: x? 2 argmin x 2RN E ( x ) = 1 2 || y x ||2 + || x ||1 E ( x, ˜ x ) = E ( x ) 1 2 || ( x ˜ x )||2 + 1 2 ⌧ || x ˜ x ||2 Proof: E ( x, ˜ x ) / 1 2 || u x ||2 + || x ||1+ cst. argmin x E ( x, ˜ x ) = S ⌧ ( u ) Theorem: E (· , ˜ x ) x ˜ x S ⌧ (u) ⌧ < 1/|| ⇤ || where u = ˜ x ⌧ ⇤( ˜ x y ) Surrogate Functionals Surrogate functional: E(·)
95. Algorithm: x ( ` +1) = argmin x E ( x, x ( ` )) x (`+1) = S 1 ⌧ ( u (`)) Initialize x (0), set ` = 0. x (0) x (1) x (2) u (`) = x (`) ⌧ ⇤( x (`) y ) Iterative Thresholding x E(·)
96. Algorithm: x ( ` +1) = argmin x E ( x, x ( ` )) x (`+1) = S 1 ⌧ ( u (`)) Remark: x (`) 7! u (`) is a gradient descent of || x y ||2 . S 1 `⌧ is the proximal step of associated to || x ||1. Initialize x (0), set ` = 0. x (0) x (1) x (2) u (`) = x (`) ⌧ ⇤( x (`) y ) Iterative Thresholding x E(·)
97. Theorem: if ⌧ < 2 / || ⇤ ||, then x (`) ! x ?. Algorithm: x ( ` +1) = argmin x E ( x, x ( ` )) x (`+1) = S 1 ⌧ ( u (`)) Remark: x (`) 7! u (`) is a gradient descent of || x y ||2 . S 1 `⌧ is the proximal step of associated to || x ||1. Initialize x (0), set ` = 0. x (0) x (1) x (2) u (`) = x (`) ⌧ ⇤( x (`) y ) Iterative Thresholding x E(·)
98. Conclusion
99. Conclusion
100. random acquisition. optimization for reconstruction. #measures sparsity Conclusion • Compressed sensing.