A Sparse Tour of Imaging Sciences

A Sparse Tour of Imaging Sciences

Talk at Dassault System

E34ded36efe4b7abb12510d4e525fee8?s=128

Gabriel Peyré

June 16, 2014
Tweet

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.