2.6k

# A Sparse Tour of Imaging Sciences

Talk at Dassault System

June 16, 2014

## Transcript

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)

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

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

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:

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

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

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 ﬁdelity 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 ﬁdelity 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 ﬁdelity 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 ﬁdelity y = x0 + w Noisy measurements:

x argmin x RQ 1 2 ||y x||2 + ||x||1 Noisy Sparse Regularization
67. ### Regularization Data ﬁdelity 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 ﬁdelity 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

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(·)

sensing.