Gabriel Peyré
June 27, 2014
500

# Des mathématiques appliquées, appliquables, et une application à l’imagerie

Exposé en l'honneur de Michel Pierre

June 27, 2014

## Transcript

1. ### c J. Malick Des mathématiques appliquées, appliquables, et une application

à l’imagerie Gabriel Peyré www.numerical-tours.com

8. ### Overview • Approximation in an Ortho-Basis • Compression and Denoising

• Sparse Inverse Problem Regularization
9. ### Orthogonal basis { m }m of L2([0, 1]d) 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
11. ### 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
12. ### 1-D Wavelet Basis Wavelets: j,n (x) = 1 2j/2 x

2jn 2j Position n, scale 2j, m = (n, j).
13. ### 1-D Wavelet Basis Wavelets: j,n (x) = 1 2j/2 x

2jn 2j Position n, scale 2j, m = (n, j).
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 2-D Fourier Basis Basis { m (x)}m of L2([0, 1]) m1 m2
15. ### 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 f(x) f, m1,m2 Fourier transform x m
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. ### 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)

22. ### Overview • Approximation in an Ortho-Basis • Compression and Denoising

• Sparse Inverse Problem Regularization
23. ### Compression by Transform-coding Image f f forward a[m] = ⇥f,

m ⇤ R transform ˜ a[m]
24. ### Compression by Transform-coding Image 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] bin T q[m] Z ˜ a[m]
25. ### Compression by Transform-coding Image 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] bin T q[m] Z ˜ a[m]
26. ### Compression by Transform-coding Image 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] bin T q[m] Z ˜ a[m]
27. ### Compression by Transform-coding Image 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] bin T q[m] Z ˜ a[m]
28. ### Compression by Transform-coding Image 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 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] bin T q[m] Z Zoom on f ˜ a[m] fR , R =0.2 bit/pixel
29. ### Compression by Transform-coding Image 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 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] bin T q[m] Z ||f fM ||2 = O(M ) =⇥ ||f fR ||2 = O(log (R)R ) Theorem: Zoom on f ˜ a[m] fR , R =0.2 bit/pixel

32. ### Noisy f Denoising thresh. f = N 1 m=0 f,

m ⇥ m ˜ f = | f, m ⇥|>T f, m ⇥ m
33. ### Noisy f 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 )
34. ### Overview • Approximation in an Ortho-Basis • Compression and Denoising

• Sparse Inverse Problems Regularization
35. ### Recovering f0 2 RN from noisy observations: y = Kf0

+ w 2 RP K : RN ! RP with P ⌧ N (missing information) Inverse Problems
36. ### Examples: Inpainting, super-resolution, . . . Recovering f0 2 RN

from noisy observations: y = Kf0 + w 2 RP K : RN ! RP with P ⌧ N (missing information) Inverse Problems f0 K K

RQ
39. ### L1 Regularization observations w coe cients image K x0 RN

f0 = x0 RQ y = Kf0 + w RP
40. ### 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
41. ### 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
42. ### 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.
43. ### Inverse Problem Regularization observations y parameter Estimator: x(y) depends only

on Observations: y = x0 + w 2 RP .
44. ### 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 .
45. ### 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
46. ### 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
47. ### J0( x ) = # { m ; xm 6=

0} “Ideal” sparsity: Sparse Priors
48. ### Sparse regularization: x? 2 argmin x 1 2 || y

x ||2 + J0( x ) J0( x ) = # { m ; xm 6= 0} “Ideal” sparsity: Sparse Priors
49. ### 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
50. ### 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
51. ### 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 ?
52. ### 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 ?
53. ### 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
54. ### 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
55. ### 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
56. ### x argmin x=y m |xm | x x = y

Noiseless Sparse Regularization Noiseless measurements: y = x0
57. ### 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
58. ### 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
59. ### y = Kf0 + w (Kf)i = ⇢ 0 if

i 2 ⌦, fi if i / 2 ⌦. Original f0 . Original Measures Inpainting Problem
60. ### y = Kf0 + w (Kf)i = ⇢ 0 if

i 2 ⌦, fi if i / 2 ⌦. Original f0 . Original Measures Inpainting Problem f ? = x ? Recovered

63. ### • Mathematical challenges: – Theoretical performance guarantees. – Fast algorithms.

Conclusion • Sparse regularization: denoising ! inverse problems.