Anisotropic Diffusion in Image Processing: Implementation, Parameters, and Segmentation

Linear and anisotropic diffusion in image processing: A study on implementation, parameters and segmentation Patricia van Marlen April 12, 2018 1

Outline Introduction Part 1: - diffusion filtering methods - numerical schemes - results: methods and schemes Part 2: - time step t: theory and results - stopping time S: theory and results Part 3: - results: segmentation Conclusion and further research 2

Introduction LUMC/TUD MRI project: MRI scanner for developing countries 3

Introduction LUMC/TUD MRI project: MRI scanner for developing countries Advantages: - affordable 4

Introduction LUMC/TUD MRI project: MRI scanner for developing countries Advantages: - affordable - transportable 5

Introduction LUMC/TUD MRI project: MRI scanner for developing countries Advantages: - affordable - transportable - always ready for use - no cooling necessary 6

Introduction LUMC/TUD MRI project: MRI scanner for developing countries Advantages: - affordable - transportable - always ready for use - no cooling necessary Disadvantage: - low SNR: noise in the images 7

Introduction LUMC/TUD MRI project: MRI scanner for developing countries Advantages: - affordable - transportable - always ready for use - no cooling necessary Disadvantage: - low SNR: noise in the images Noise removal: diffusion filtering methods 8

Diffusion filtering methods: linear diffusion Diffusion equation = ( ) tu div D u D is the diffusion coefficient u is gray value Neumann boundary conditions: = / 0 u n 9

Diffusion filtering methods: linear diffusion Diffusion equation = ( ) tu div D u D is the diffusion coefficient u is gray value Neumann boundary conditions: = / 0 u n D=1, leads to the heat equation = tu u Disadvantage: space-invariant blurring 10

Diffusion filtering methods: anisotropic diffusion Diffusion equation = ( ) tu div D u D is the diffusion coefficient u is gray value Neumann boundary conditions: = / 0 u n 11

Diffusion filtering methods: anisotropic diffusion Diffusion equation = ( ) tu div D u D is the diffusion coefficient u is gray value Neumann boundary conditions: = / 0 u n D=c(| u|), leads to = ( (| |) ) tu div c u u 12

Diffusion filtering methods: Anisotropic diffusion Intraregional blurring: Within a region, u is small So: if u is small, we need c(| u|) to be high Edge preservation: Close to a boundary, u is large So: if u is large, we need c (| u|) to be close to zero 13

Diffusion filtering methods: Anisotropic diffusion 1 = (| |) c u Well-posed method: | | u 14

Diffusion filtering methods: Anisotropic diffusion 1 = (| |) c u Well-posed method: | | u 1 1 2 = (| |) | | c u u Ill-posed method: N N 15

Diffusion filtering methods: Anisotropic diffusion 1 = (| |) c u Well-posed method: | | u 1 1 2 = (| |) | | c u u Ill-posed method: N N Fulfill the threshold property: 1 = (| |) c u Perona-Malik 1: + 2 1 (| | / ) K u 2 = (| |/ ) u K (| |) c u e Perona-Malik e: 16

Numerical schemes: Forward Time, Central Space (FTCS) Time-derivative: Euler forward Space-derivatives: central difference Accuracy: O( t)+O( x2) Set x= y=1 17

Numerical schemes: Forward Time, Central Space (FTCS) Linear diffusion filtering = + + + + + 1 k k k k k k k ( 4 ) u u t u u u u u + + , , 1, 1, , , 1 , 1 i j i j i j i j i j i j i j 18

Numerical schemes: Forward Time, Central Space (FTCS) Linear diffusion filtering = + + + + + 1 k k k k k k k ( 4 ) u u t u u u u u + + , , 1, 1, , , 1 , 1 i j i j i j i j i j i j i j Anisotropic diffusion filtering = + 1 k k u u , , i j i j + 2 2 k k k k k k ( (| t c | )( ) (| | )( )) u u u c u u u + + 1/2, 1, , 1/2, 1, , i j i j i j i j i j i j ( ) + 2 2 k k k k k k (| | )( ) (| | )( ) t c u u u c u u u + + , 1/2 , 1 , , 1/2 , 1 , i j i j i j i j i j i j 19

Numerical schemes: Forward Time, Central Space (FTCS) Linear diffusion filtering = + + + + + 1 k k k k k k k ( 4 ) u u t u u u u u + + , , 1, 1, , , 1 , 1 i j i j i j i j i j i j i j Anisotropic diffusion filtering = + 1 k k u u , , i j i j + 2 2 k k k k k k ( (| t c | )( ) (| | )( )) u u u c u u u + + 1/2, 1, , 1/2, 1, , i j i j i j i j i j i j ( ) + 2 2 k k k k k k (| | )( ) (| | )( ) t c u u u c u u u + + , 1/2 , 1 , , 1/2 , 1 , i j i j i j i j i j i j = ( ) u G s u Addition of Gaussian kernel: 20

Numerical schemes: Forward Time, Central Space (FTCS) Linear diffusion filtering = + + + + + 1 k k k k k k k ( 4 ) u u t u u u u u + + , , 1, 1, , , 1 , 1 i j i j i j i j i j i j i j Anisotropic diffusion filtering with Gaussian kernel = + 1 k k u u , , i j i j + 2 2 k k k k k k ( (| t c | )( ) (| | )( )) u u u c u u u + + 1/2, 1, , 1/2, 1, , i j i j i j i j i j i j ( ) + 2 2 k k k k k k (| | )( ) (| | )( ) t c u u u c u u u + + , 1/2 , 1 , , 1/2 , 1 , i j i j i j i j i j i j = ( ) u G s u Addition of Gaussian kernel: 21

Numerical schemes: Additive Operator Splitting (AOS) The AOS scheme in 2D 1 ( ) 2 1 + = 1 k k k 2 ( ) U I tA U U 2 = 1 22

Numerical schemes: Additive Operator Splitting (AOS) The AOS scheme in 2D 1 ( ) 2 1 + = 1 k k k 2 ( ) U I tA U U 2 = 1 We have two dimensions, so Axand Ay 23

Numerical schemes: Additive Operator Splitting (AOS) The AOS scheme in 2D 1 ( ) 2 1 + = 1 k k k 2 ( ) U I tA U U 2 = 1 We have two dimensions, so Axand Ay Figure 2: a) original numbering b) transposed numbering 24

Results: methods and schemes Test problem: Gaussian noise with mean 0 and SD 0.06 Figure 3: test problem with and without noise 25

Results: methods and schemes Test problem: Gaussian noise with mean 0 and SD 0.06 Figure 3: test problem with and without noise Difference measure is the SD of the noise in the image: 1 || original u n = || u 2 filtered 26

Results: methods and schemes FTCS scheme method standard deviation Perona-Malik 1 0.0592 Perona-Malik e 0.0592 Perona-Malik 1 with GK 0.0594 Perona-Malik e with GK 0.0594 Well-posed method 0.0595 Ill-posed method 0.0595 Heat equation 0.0597 Table 1: results of the diffusion filtering methods for FTCS scheme 27

Results: methods and schemes FTCS scheme method standard deviation Perona-Malik 1 0.0592 Perona-Malik e 0.0592 Perona-Malik 1 with GK 0.0594 Perona-Malik e with GK 0.0594 Well-posed method 0.0595 Ill-posed method 0.0595 Heat equation 0.0597 Table 1: results of the diffusion filtering methods for FTCS scheme Figure 4: best result of Perona-Malik 1 28

Results: methods and schemes FTCS scheme method standard deviation Perona-Malik 1 0.0592 Perona-Malik e 0.0592 Perona-Malik 1 with GK 0.0594 Perona-Malik e with GK 0.0594 Well-posed method 0.0595 Ill-posed method 0.0595 Heat equation 0.0597 Table 1: results of the diffusion filtering methods for FTCS scheme Figure 5: unstable result of well-posed method 29

Results: methods and schemes AOS scheme method standard deviation Perona-Malik 1 with GK 0.0590 Perona-Malik e 0.0590 Perona-Malik 1 0.0592 Well-posed method 0.0593 Perona-Malik e with GK 0.0596 Ill-posed method 0.0596 Table 2: results of the diffusion filtering methods for AOS scheme 30

Results: methods and schemes AOS scheme method standard deviation Perona-Malik 1 with GK 0.0590 Perona-Malik e 0.0590 Perona-Malik 1 0.0592 Well-posed method 0.0593 Perona-Malik e with GK 0.0596 Ill-posed method 0.0596 Table 2: results of the diffusion filtering methods for AOS scheme Figure 6: best result of Perona-Malik 1 with GK 31

Results: methods and schemes conclusions - the AOS scheme gave slightly better results, but due to calculation time FTCS is used for the following parts 32

Results: methods and schemes conclusions - the AOS scheme gave slightly better results, but due to calculation time FTCS is used for the following parts - Perona-Malik 1 seems the best method and is also chosen for following parts, it is given by: 1 = (| |) c u + 2 1 (| | / ) K u 33

Results: methods and schemes conclusions - the AOS scheme gave slightly better results, but due to calculation time FTCS is used for the following parts - Perona-Malik 1 seems the best method and is also chosen for following parts, it is given by: 1 = (| |) c u + 2 1 (| | / ) K u - the well-posed and ill-posed method showed instability for FTCS scheme, but were stable for AOS scheme 34

Results: methods and schemes conclusions - the AOS scheme gave slightly better results, but due to calculation time FTCS is used for the following parts - Perona-Malik 1 seems the best method and is also chosen for following parts, it is given by: 1 = (| |) c u + 2 1 (| | / ) K u - the well-posed and ill-posed method showed instability for FTCS scheme, but were stable for AOS scheme - parameter choice has a big influence on the outcomes 35

Parameter study Investigation of - parameter K, used in the Perona-Malik functions - time step size t - number of time steps T 36

Parameter study Investigation of - parameter K, used in the Perona-Malik functions - time step size t - number of time steps T 37

Parameter study time step size t So far t was constant Adaptive time step for Perona-Malik methods ( ) k 1 | Re( / ) | u t = + ( ) k exp max t a b ( ) k Re( ) u ( ) k u change with =4, a=0.25 and b=0.75 38

Parameter study time step size t So far t was constant Adaptive time step for Perona-Malik methods ( ) k 1 | Re( / ) | u t = + ( ) k exp max t a b ( ) k Re( ) u ( ) k u change with =4, a=0.25 and b=0.75 Initial iterations (large uchange): t ~ 0.0625 Later iterations (small uchange): t ~ 0.25 So: an increasing function 39

Parameter study time step size t Figure 7: Optimal time step for Perona-Malik 1 Initial iterations: t ~ 0.25 Later iterations: t ~ 0 So: a decreasing function 40

Parameter study time step size t Adaptive time step for Perona-Malik methods ( ) k 1 | Re( / ) | u t = + ( ) k exp max t a b ( ) k Re( ) u with =4, a=0.25 and b=0.75 Initial iterations: t ~ 0.0625 Later iterations: t ~ 0.25 So: an increasing function 41

Parameter study time step size t Adaptive time step for Perona-Malik methods ( ) k | Re( / ) | u t 1 ( ) m k = exp max t a b od ( ) k Re( ) u with =4, a=1 and b=1 Initial iterations: t ~ 0.25 Later iterations: t ~ 0 So: a decreasing function 42

Results: parameters time step size t Behaviour of Blue: optimal t Red: original method Green: modified method Figure 8: behaviour of various time steps 43

Results: parameters time step size t Figure 9: results for K=0.05 Figure 8: behaviour of various time steps 44

Results: parameters time step size t Figure 9: results for K=0.05 Figure 8: behaviour of various time steps Figure 10: results for K=0.5 45

Parameter study time step size t Adaptive time step for well-posed and ill-posed method 1 k k 1.0001 ( 0) u u k = ( ) k ( 1), k t q t if , , i j i j , , i j i j with 0<q<1 46

Parameter study time step size t Adaptive time step for well-posed and ill-posed method 1 k k 1.0001 ( 0) u u k = ( ) k ( 1), k t q t if , , i j i j , , i j i j with 0<q<1 Method is based on constant total pixel value (TPV), if TPV starts to increase the time step is reduced 47

Results: parameters time step size t Adaptive time step for well-posed and ill-posed method 1 k k 1.0001 ( 0) u u k = ( ) k ( 1), k t q t if , , i j i j , , i j i j with 0<q<1 Figure 11: results for adaptive time step for various q 48

Parameter study stopping time S Optimal stopping time: number of time steps at which the image is at its best and the diffusion should stop 49

Parameter study stopping time S Optimal stopping time: number of time steps at which the image is at its best and the diffusion should stop Correlation stopping time S ( ) = arg min (0) ( ), ( ) u t u t corr u corr t 50