Semi-discretization of Euler-Lagrange equation: du dt = ( µL ( u ) − A 2 ) u + Ab, where L ( u ) a nonlinear operator. • Explicit integration method (Euler): CFL condition imposes tiny time steps ⇒ method expensive • Semi-implicit integration method: [ I − τ ( µL ( u j − 1 ) − A 2 )] u j = u j − 1 + τAb, where u j ≈ u ( jτ ). System expensive to solve and small time steps τ required ⇒ method expensive 41
Use Perona-Malik operator D to define a nonlinear prolongation operator for the multilevel scheme. Discretization of ∂u ∂t = D ( u ) gives du dt = L ( u ) u, u ∈ W i , u (0) = P i u i − 1 . May integrate for a few, say ≤ 10, (non-tiny) time steps by explicit method. 42
Cascadic multilevel methods 43
Consider � t ∈ Ω Ω k ( s, t ) x ( s ) ds = b ( t ) , Let W 1 ⊂ W 2 ⊂ . . . ⊂ W ℓ ⊂ L 2 (Ω) nested subspaces R i : L 2 (Ω) → W i restriction operator Q ∗ prolongation operator, e.g., Q ∗ i = R ∗ → : W i L 2 (Ω) i i ˆ b i = R i ˆ A i = R i AQ ∗ b, b i = R i b, i P i : W i − 1 → W i prolongation operator 44
Algorithm 1 Multilevel Algorithm Input: A , b , δ , ℓ ≥ 1 (number of levels); Output: approximate solution ˜ x ∈ W ℓ ; Determine A i and b i for 1 ≤ i ≤ ℓ ; x 0 := 0 ; for i := 1 , 2 , . . . , ℓ do x i, 0 := P i x i − 1 ; ∆ x i,m i := IM ( A i , b i − A i x i, 0 ) ; Correction step: x i := x i, 0 + ∆ x i,m i ; endfor x := x ℓ ; ˜ 45
Noise reducing restrictions R i : Solve weighted local least-squares problems (inspired by Buades et al.): 1D problems: Define M i : W i → W i − 1 by solving � 2 ω (2 j ) � x (2 j + s ) � min − ( a 0 + a 1 s ) ( s ) ∀ j i i a 0 ,a 1 s ∈{ 0 , ± 1 } where � � 2 � ω (2 j ) � x (2 j + s ) − x (2 j ) − γ ( s ) := exp , γ > 0 i i i 46
Solution { ˆ a 1 } . a 0 , ˆ Let x ( j ) i − 1 := ˆ a 0 for all j . Define 1 ≤ i < ℓ. R i = M i +1 M i +2 . . . M ℓ , 2D problems: 47
Pixel mask for 2D problems (2j−1,2k−1) (2j−1,2k) (2j−1,2k+1) (2j,2k) (2j,2k−1) (2j,2k+1) (2j+1,2k−1) (2j+1,2k) (2j+1,2k+1) 48
Assume that � ˆ b − b � = δ. Then, under some simplifying assumptions, such as γ = 0, one can show that for 1D problems one can expect the projected errors to satisfy the bounds c i = 1 � ˆ √ b i − b i � ≤ c i δ, 1 ≤ i < ℓ, 3 c i +1 , c ℓ = 1 . For 2D problems, c i = 1 3 c i +1 . 49
Edge preserving nonlinear prolongation The operator P i consists of two parts: x i − 1 ∈ W i − 1 → L i x i − 1 ∈ W i piecewise linear interpolation L i x i − 1 is smoothed by solving an IBVP for a (discretized) Perona-Malik diffusion equation, 2 � � ∂x ∂τ = ∂ ∂ ∂ � � , ψ ′ a ≤ s ≤ b, ∂sx ∂sx x = x ( τ, s ) , � � ∂s � � � � over a short τ -interval with ψ ′ ( s ) = ρ/ ( s + ρ ), ρ > 0. This removes noise and preserves edges. 50
Theorem: Let � ˆ b i − b i � ≤ c i δ i , c i > 1 , 1 ≤ i ≤ ℓ. Assume that the P i are linear(ized) and that R ( P i ) ⊂ R ( A ∗ 2 ≤ i ≤ ℓ. i ) , Then || ˆ x i − x i,m i ( δ i ) || = 0 , 1 ≤ i ≤ ℓ, lim sup δ i ց 0 || ˆ b i − b i ||≤ c i δ i i ˆ x i = A † where ˆ b i . Thus, the multilevel is a regularization method. 51
Examples: Noise reduction by Perona-Malik integration Example 2: Signal contaminated by 1% Gaussian noise. The noise-free signal is smooth. 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0 −0.01 0 50 100 150 200 250 300 350 400 450 500 52
Example 2: Signal denoised by integrating the Perona-Malik equation 10 steps of size ∆ τ = 0 . 3. 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0 −0.01 0 50 100 150 200 250 300 350 400 450 500 53
Example 3: Signal contaminated by 1% Gaussian noise. The noise-free signal is piecewise constant. 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0 −0.01 0 100 200 300 400 500 600 54
Example 3: Signal denoised by integrating the Perona-Malik equation 10 steps of size ∆ τ = 0 . 3. 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0 −0.01 0 100 200 300 400 500 600 55
Example 4: Signal contaminated by 10% Gaussian noise. The noise-free signal is piecewise constant. 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0 −0.01 −0.02 0 100 200 300 400 500 600 56
Example 4: Signal denoised by integrating the Perona-Malik equation. 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0 −0.01 0 100 200 300 400 500 600 57
Example 5: Fredholm integral equation of the 1st kind � π 0 exp( − st ) x ( t ) dt = 2sinh( s ) 0 ≤ s ≤ π , 2; s same as in Example 1. Discretization by Galerkin method using 512 piecewise constant test and trial functions. Relative error ν . Matrix sizes for 5-level method: R 32 × 32 , A 2 ∈ R 64 × 64 , A 3 ∈ R 128 × 128 , ∈ A 1 R 256 × 256 , A = A 5 ∈ R 512 × 512 . ∈ A 4 Minimal residual Iterative method for symmetric problems. 58
P i L i ℓ ν � ˜ x − ˆ x � / � ˆ x � # iter � ˜ x − ˆ x � / � ˆ x � # iter 1 · 10 − 2 3 . 51 · 10 − 2 1 3 1 · 10 − 2 3 . 26 · 10 − 2 3 . 41 · 10 − 2 2 3 1 3 1 1 · 10 − 2 3 . 10 · 10 − 2 3 . 41 · 10 − 2 3 3 1 1 3 1 1 1 · 10 − 2 2 . 51 · 10 − 2 3 . 35 · 10 − 2 5 3 2 1 1 1 3 2 1 1 1 1 · 10 − 3 3 . 53 · 10 − 2 1 3 1 · 10 − 3 3 . 34 · 10 − 2 3 . 57 · 10 − 2 2 3 1 3 1 1 · 10 − 3 3 . 08 · 10 − 2 3 . 56 · 10 − 2 3 3 2 1 3 2 1 1 · 10 − 3 2 . 00 · 10 − 2 3 . 50 · 10 − 2 5 3 2 2 3 1 3 2 2 2 1 59
Example 6: Signal deblurring by solution of Fredhom integral equation of the 1st kind with a Gaussian convolution kernel. Discretization gives symmetric Toeplitz matrices. Relative error 10%. Matrix sizes for 3-level method: A 1 ∈ R 128 × 128 , A 2 R 256 × 256 , A = A 3 ∈ R 512 × 512 ∈ κ ( A ) = 1 · 10 17 . Iterative method: MR-II 60
Example 6: Signal contaminated by 10% Gaussian noise and blur. The noise-free signal is piecewise constant. 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0 −0.01 −0.02 0 100 200 300 400 500 600 61
Example 6: Restored signal by multilevel method with 3 levels. 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0 −0.01 0 100 200 300 400 500 600 62
P i L i � ˜ x − ˆ x � / � ˆ x � � ˜ x − ˆ x � / � ˆ x � ℓ # iter # iter 1 . 02 · 10 − 1 1 1 8 . 70 · 10 − 2 8 . 78 · 10 − 2 2 1 1 1 1 8 . 08 · 10 − 2 8 . 88 · 10 − 2 3 4 1 1 4 1 1 63
Application to image restoration 64
Example 7. Blur- and noise-free 512 × 512-pixel image. 65
Restoration of image with 5% noise, Gaussian blur, determined by 2 RRGMRES iterations. 255 Quantitative measure: PSNR= 20 log 10 u || = 27 . 98 || u ℓ − ˆ dB. 66
RRGMRES CGNR ℓ ν PSNR # iter PSNR # iter 1 · 10 − 2 1 29.77 3 32.36 9 1 · 10 − 2 2 31.14 3 2 34.09 5 7 5 · 10 − 2 1 27.98 2 28.74 3 5 · 10 − 2 2 28.80 2 1 29.90 2 2 67
Example 8. Original lizard image and image contaminated by 10% noise and motion blur. 68
Lizard restorations: CGNR CGNR-based 3-level method PSNR=25 . 11 PSNR=26 . 11 69
Example 9: Original 256 × 256 image. 70
Image contaminated by Gaussian blur and 1% noise. 71
Restoration by solving Euler-Lagrange equation, Perona-Mailk regularization operator, 500 time-steps, 1000 mvp, PSNR=28 . 84. 72
Restoration by 3-level method, Perona-Malik-based prolongation operator, 4 mvp, PSNR=35 . 85. 73
Carrying out too many iterations may give rise to instability and produce unexpected results ... 74
such as ... 75
or ... 76
Application to image segmentation 77
Algorithm 2 Enhanced Multilevel Algorithm Input: A , b , δ , ℓ ≥ 1 (number of levels); x := x ℓ ∈ W ℓ ; Output: approximate solution ˜ Determine A i and b i for 1 ≤ i ≤ ℓ ; x 0 := 0 ; φ 0 := initial contour; for i := 1 , 2 , . . . , ℓ do x i, 0 := P i x i − 1 ; φ i, 0 := S i φ i − 1 ; ∆ x i,m i := IM ( A i , b i − A i x i, 0 ) ; Correction step: x i := x i, 0 + ∆ x i,m i ; Segmentation step: φ i := GAC ( φ i, 0 , x i ) ; endfor 78
Example 10. Segmentation of noise- and blur-free image. 79
Available image with 10% noise and Gaussian blur restored by CGNR and then segmented. 80
Available contaminated images restored and segmented by 3-level CGNR-based method. Segmentation on level 2. 81
Segmentation on level 3. 82
Wavelet-based multilevel methods 83
Let φ and and ψ be scaling and wavelet functions, respectively. Define Ψ j,k := 2 j/ 2 Ψ(2 j · − k ) , φ j,k := 2 j/ 2 φ (2 j · − k ) , where j and k the scale and space parameters, respectively. Then � � � x ∈ L 2 (R) . x = ( x, φ 0 ,k ) φ 0 ,k + ( x, Ψ j,k )Ψ j,k , k ∈ Z j ≥ 0 k ∈ Z Let Ψ − 1 ,k := φ 0 ,k and define the subspaces W j = span { Ψ j,k } | k |≤ k ( j ) . Here k ( j ) finite since we consider a finite interval [ a, b ]. 84
Define Subspace: i � U i := W l ⊂ L 2 ([ a, b ]) l = − 1 Orthogonal projector: Q i : L 2 ([ a, b ]) → U i . Restricted operator A i : U i → U i : i = − 1 , 0 , 1 , . . . . A i := Q i AQ i , 85
Prolongation operator P i : U i − 1 → U i : i � � P i x := x j,k Ψ j,k ˜ j = − 1 | k |≤ k ( j ) with j ≤ i − 1 , | k | ≤ k ( i ) , x j,k , for x j,k = ˜ | k | ≤ k ( i ) . 0 , for j = i, Then i ) = N ( A i ) ⊥ . i − 1 ) ) ⊂ R ( A ∗ R ( P i | R ( A ∗ 86
Theorem: Apply multilevel method based on the CGNR or MR-II iterative methods. Then the computed approximate solution on level i lives in N ( A i ) ⊥ . Example 11. Consider the integral equation � π/ 2 κ ( s, t ) x ( s ) ds = sin 2 (3 t ) − 0 . 9 t 3 + t 2 , 0 ≤ t ≤ π/ 2 , 0 with κ ( s, t ) := exp( s cos( t )). 87
Order n i of the matrices A i . Submatrix A 1 A 2 A 3 A 4 A 5 A 6 A 7 A 8 n i 28 50 88 158 292 554 1071 2098 88
Right-hand sides without noise and with 10% noise 89
Exact and computed approximate solutions for 10% noise 90
Exact and computed approximate solutions for 1% noise 91
Exact and computed approximate solutions for 0 . 001% noise 92
CGNR ML-CGNR � x 8 ,k 8 − ˆ x � � x k − ˆ x � noise k k i Accel. � ˆ x � � ˆ x � 1 · 10 − 1 5 . 78 · 10 − 1 5 . 66 · 10 − 1 2 2 , 0 , 0 , 0 , 0 , 0 , 0 , 0 17.6 1 · 10 − 2 5 . 35 · 10 − 1 5 . 28 · 10 − 1 3 3 , 0 , 0 , 0 , 0 , 0 , 0 , 0 23.5 1 · 10 − 3 3 . 41 · 10 − 1 3 . 24 · 10 − 1 5 5 , 0 , 0 , 0 , 0 , 0 , 0 , 0 35.2 1 · 10 − 4 3 . 32 · 10 − 1 3 . 14 · 10 − 1 5 5 , 0 , 0 , 0 , 0 , 0 , 0 , 0 35.2 1 · 10 − 5 3 . 49 · 10 − 2 3 . 46 · 10 − 2 8 8 , 0 , 0 , 0 , 0 , 0 , 0 , 0 52.6 93
Alternating iterative methods 94
Continuous image degradation model � x ∈ Ω f ( x ) = Ω h ( x, y )ˆ u ( y ) dy + η ( x ) , with h the point spread function and f the available blur- and noise- contaminated image. Determine the blur- and noise-free image ˆ u by solving min u,w J ( u, w ) , where 95
��� � 2 � Ω h ( x, y ) u ( y ) dy − f ( x ) J ( u, w ) := Ω + α ( L ( u − w )( x )) 2 + β |∇ w ( x ) | � dx, using an alternating iterative method. Here L is a differential operator, e.g., the Perona-Malik operator and α > 0, β > 0 are regularization parameters. 96
Discrete setting: u ( i ) S h ( w ( i − 1) ) := argmin u ∈K ℓ {� Hu − f δ � 2 = + α � L ( i − 1) ( u − w ( i − 1) ) � 2 } , S tv ( u ( i ) ) := argmin w ∈K ℓ {� L ( i − 1) ( w − u ( i ) ) � 2 + β � w � tv } , w ( i ) = for i = 1 , 2 , 3 , . . . , where � · � tv is a discrete TV-norm, L ( i − 1) is a discretization of the operator L , and K ℓ a Krylov subspace. 97
Available blur- and noise-contaminated image (Gaussian blur, 50% noise. 98
Restored image. 99
Available blur- and noise-contaminated image (Gaussian blur, 10% noise. 100
Recommend
More recommend
