ANOMALOUS DIFFUSION, DILATION, AND EROSION IN IMAGE PROCESSING joint work with Sophia Vorderw¨ ulbecke & Bernhard Burgeth SIAM CSE 2019 (MS 62) | February 25, 2019 Andreas Kleefeld J¨ ulich Supercomputing Centre, Germany Member of the Helmholtz Association
TABLE OF CONTENTS Part 1: Introduction & motivation Part 2: Anomalous diffusion Part 3: Modified dilation & erosion Part 4: Numerical results Part 5: Summary & outlook SIAM CSE 2019 (MS 62) | February 25, 2019 Member of the Helmholtz Association Andreas Kleefeld
Part I: Introduction & motivation SIAM CSE 2019 (MS 62) | February 25, 2019 Member of the Helmholtz Association Andreas Kleefeld
INTRODUCTION & MOTIVATION General idea Time-dependent partial differential equations (PDEs) arise naturally in image processing. For example: convolution of image with Gaussian kernel which is equivalent to solving a linear diffusion equation. Other PDEs: dilation/erosion (evolution equations). Can serve as building blocks for higher morphological operations (opening, closing, gradients) or deblurring filters. SIAM CSE 2019 (MS 62) | February 25, 2019 Member of the Helmholtz Association Andreas Kleefeld
INTRODUCTION & MOTIVATION What is new? Different type of generalization of an evolution equation. ∂ α Temporal derivative of fractional order α : ∂ t α with α ∈ ( 0 , 2 ) . Definition of the fractional derivative as an extension of integration concatenated with regular differentiation (Caputo). Global information are considered. Also interesting for other applications. Up to now this approach was only considered for specific fractional orders as α = 1 / 2 and not for morphological operations. SIAM CSE 2019 (MS 62) | February 25, 2019 Member of the Helmholtz Association Andreas Kleefeld
Part II: Anomalous diffusion SIAM CSE 2019 (MS 62) | February 25, 2019 Member of the Helmholtz Association Andreas Kleefeld
ANOMALOUS DIFFUSION Mathematical model Diffusion equation: c ∂ α ∂ t α u = div ( κ grad u ) , where κ is a constant. Caputo fractional derivative: � t c ∂ α u ( m + 1 ) ( τ ) 1 ∂ t α u = ( t − τ ) α − m d τ , Γ( m + 1 − α ) 0 where m = ⌊ α ⌋ and 0 < α < 1 or 1 < α < 2 . Initial condition(s): given gray-value image and in case of super-diffusion we need a second initial condition. Boundary condition: homogeneous Neumann. SIAM CSE 2019 (MS 62) | February 25, 2019 Member of the Helmholtz Association Andreas Kleefeld
ANOMALOUS DIFFUSION Space discretization 2D-grid with h = 1 and M × N grid points. Approximation of Laplace operator with centered differences for interior nodes: κ ( u i + 1 , j + u i − 1 , j + u i , j + 1 + u i , j − 1 − 4 u i , j ) . Homogeneous Neumann boundaries for exterior nodes. SIAM CSE 2019 (MS 62) | February 25, 2019 Member of the Helmholtz Association Andreas Kleefeld
ANOMALOUS DIFFUSION Time discretization Grid of the form t k = k ∆ t , k = 0 , . . . , P with grid size ∆ t = T / P . Approximation of Caputo derivative by Gr¨ unwald-Letnikov formula: k + 1 m t = t k + 1 C ∂ α u ( t k + 1 ) n − α � � c ( α ) � � u k + 1 − ℓ Γ( n − α + 1 ) u ( m ) ( x i , y j ) , ≈ − � ℓ i , j ∂ t α � x =( x i , y j ) ℓ = 0 n = 0 where � 1 − 1 + α � = (∆ t ) − α , c ( α ) c ( α ) c ( α ) = k − 1 . 0 k k SIAM CSE 2019 (MS 62) | February 25, 2019 Member of the Helmholtz Association Andreas Kleefeld
ANOMALOUS DIFFUSION Numerical schemes Explicit: k + 1 m ( t k + 1 ) n − α � c ( α ) � u k + 1 − ℓ Γ( n − α + 1 ) u ( m ) ( x i , y j ) − ℓ i , j ℓ = 0 n = 0 u k i + 1 , j + u k i − 1 , j + u k i , j + 1 + u k i , j − 1 − 4 u k � � = κ . i , j Implicit: k + 1 m ( t k + 1 ) n − α � c ( α ) � u k + 1 − ℓ Γ( n − α + 1 ) u ( m ) ( x i , y j ) − ℓ i , j ℓ = 0 n = 0 � � u k + 1 i + 1 , j + u k + 1 i − 1 , j + u k + 1 i , j + 1 + u k + 1 i , j − 1 − 4 u k + 1 = κ . i , j SIAM CSE 2019 (MS 62) | February 25, 2019 Member of the Helmholtz Association Andreas Kleefeld
ANOMALOUS DIFFUSION Numerical schemes Explicit: u k + 1 = A u k − b ex with A = α I MN + (∆ t ) α κ · D 2 . Implicit: B u k + 1 = b im with B = − (∆ t ) − α I MN + κ · D 2 . D 2 is the 2D-Laplacian and b ex and b im are given by � k + 1 m � ( t k + 1 ) n − α u k + 1 − ℓ − c ( α ) � � Γ( n − α + 1 ) u ( m ) ( x i , y j ) b ex = (∆ t ) α , ℓ ℓ = 2 n = 0 � k + 1 m � ( t k + 1 ) n − α b im = − α (∆ t ) − α u k + u k + 1 − ℓ − c ( α ) � � Γ( n − α + 1 ) u ( m ) ( x i , y j ) . ℓ ℓ = 2 n = 0 SIAM CSE 2019 (MS 62) | February 25, 2019 Member of the Helmholtz Association Andreas Kleefeld
Part III: Modified dilation & erosion SIAM CSE 2019 (MS 62) | February 25, 2019 Member of the Helmholtz Association Andreas Kleefeld
MODIFIED DILATION & EROSION Mathematical model & discretization Dilation & erosion equation: �� ∂ u � 2 � 2 c ∂ α � ∂ u ∂ t α u = ± + . ∂ x ∂ y Approximation of Caputo fractional derivative as before. Approximation in space by first-order finite difference scheme of Rouy-Tourin: max ( − u i , j + u i − 1 , j , u i + 1 , j − u i , j , 0 ) 2 + max ( − u i , j + u i , j − 1 , u i , j + 1 − u i , j , 0 ) 2 � 1 / 2 . � SIAM CSE 2019 (MS 62) | February 25, 2019 Member of the Helmholtz Association Andreas Kleefeld
MODIFIED DILATION & EROSION Numerical schemes As before, we obtain an iterative scheme of the form u k + 1 = α u k + (∆ t ) α b dt ± (∆ t ) α � b 2 dx + b 2 dy , where k + 1 t − α u k + 1 − l + c ( α ) � k + 1 Γ( 1 − α ) u 0 b dt = − l l = 2 and the i -th, j -th entry of b dx is given by � − u k i , j + u k i − 1 , j , u k i + 1 , j − u k � max i , j , 0 and b dy analogously. SIAM CSE 2019 (MS 62) | February 25, 2019 Member of the Helmholtz Association Andreas Kleefeld
Part IV: Numerical results SIAM CSE 2019 (MS 62) | February 25, 2019 Member of the Helmholtz Association Andreas Kleefeld
NUMERICAL RESULTS Stability Linear test problem: c ∂ α u ( t ) = λ u ( t ) , λ ∈ C , ∂ t α u ( 0 ) = u 0 for 0 < α ≤ 1 , and additionally u ′ ( 0 ) = u 1 for 1 < α < 2 . Explicit method: C \{ ( 1 − z ) α / z : | z | ≤ 1 } . Implicit method: C \{ ( 1 − z ) α : | z | ≤ 1 } . SIAM CSE 2019 (MS 62) | February 25, 2019 Member of the Helmholtz Association Andreas Kleefeld
NUMERICAL RESULTS Stability 1.5 1.5 1.5 1 1 1 0.5 0.5 0.5 Imag(z) Imag(z) Imag(z) 0 0 0 -0.5 -0.5 -0.5 -1 -1 -1 -1.5 -1.5 -1.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 Real(z) Real(z) Real(z) 1.5 1.5 1.5 1 1 1 0.5 0.5 0.5 Imag(z) Imag(z) Imag(z) 0 0 0 -0.5 -0.5 -0.5 -1 -1 -1 -1.5 -1.5 -1.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 Real(z) Real(z) Real(z) Figure: Stability regions for explicit Euler method using parameters α = 0 . 4, α = 0 . 6, and α = 0 . 8 (first row) and α = 1 . 0, α = 1 . 2, and α = 1 . 4 (second row). SIAM CSE 2019 (MS 62) | February 25, 2019 Member of the Helmholtz Association Andreas Kleefeld
NUMERICAL RESULTS Stability 2 2 2 1.5 1.5 1.5 1 1 1 0.5 0.5 0.5 Imag(z) Imag(z) Imag(z) 0 0 0 -0.5 -0.5 -0.5 -1 -1 -1 -1.5 -1.5 -1.5 -2 -2 -2 0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3 Real(z) Real(z) Real(z) 2 2 2 1.5 1.5 1.5 1 1 1 0.5 0.5 0.5 Imag(z) Imag(z) Imag(z) 0 0 0 -0.5 -0.5 -0.5 -1 -1 -1 -1.5 -1.5 -1.5 -2 -2 -2 0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3 Real(z) Real(z) Real(z) Figure: Stability regions for implicit Euler method using parameters α = 0 . 4, α = 0 . 6, and α = 0 . 8 (first row) and α = 1 . 0, α = 1 . 2, and α = 1 . 4 (second row). SIAM CSE 2019 (MS 62) | February 25, 2019 Member of the Helmholtz Association Andreas Kleefeld
NUMERICAL RESULTS Stability Interval of stability is ( − 2 α , 0 ) . Implicit Euler method is A -stable for 0 < α ≤ 1 whereas we loose this property for 1 < α < 2 . Could investigate A ( θ ) stability, where θ ≤ π/ 2 will depend on α . We obtain the θ angles (in degrees ◦ ) 90, 81, 72, 63, 54, 45, 36, 27, 18, and 9 for the parameters α = 1 . 0, 1 . 1, 1 . 2, 1 . 3, 1 . 4, 1 . 5, 1 . 6, 1 . 7, 1 . 8, and 1 . 9, respectively. Hence, it appears to be that θ is given by ( 2 − α ) · 90 ◦ for 1 ≤ α < 2 (the proof remains open). SIAM CSE 2019 (MS 62) | February 25, 2019 Member of the Helmholtz Association Andreas Kleefeld
NUMERICAL RESULTS Convergence Homogeneous initial conditions: convergence order 1 Non-homogenous initial conditions: convergence order depends on α Calculation of error: c ∂ α u ( t ) = t 2 , u ( 0 ) = 0 , 0 ≤ t ≤ 1 , 1 < α ≤ 2 ∂ t α with exact solution u ( t ) = Γ( 3 + α ) t 2 + α . Γ( 3 ) Estimated convergence order (EOC): EOC = log ( E ∆ t / E ∆ t / 2 ) , where E ∆ t = | u ( 1 ) − ˜ u ∆ t ( 1 ) | . log ( 2 ) SIAM CSE 2019 (MS 62) | February 25, 2019 Member of the Helmholtz Association Andreas Kleefeld
Recommend
More recommend
