TubeTK/Anisotropic Hybrid Diffusion with Continuous Switch

This module is an implementation of diffusion based smoothing technique developed by Mendrik et al.

Mendrik AM, Vonken EJ, Rutten A, Viergever MA, van Ginneken B. Noise reduction in computed tomography scans using 3-d anisotropic hybrid diffusion with continuous switch. IEEE Trans Med Imaging. 2009 Oct;28(10):1585-94.

Algorithm synopsis

This algorithm is based on anisotropic non-linear diffusion. The technique combines edge-preserving noise reduction while enhancing local structures. This algorithm uses a hybrid approach that combines the advantages of EED ( Edge enhancing diffusion ) and CED ( Cohnerence enhancing diffusion ).

EED focuses on edge preservation and enhancement. In EED, strong smoothing is applied along the direction of the edge while the strength of the smoothing along the other perpendicular directions depends on the gradient. The higher the gradient the lower the smoothing strength would be. Applying EED to a medical image would enhance boundaries of larger organs but would blur vessels and smaller structures.

On the other hand, CED is designed to to connect lines and improve flow-like structures. Running CED on medical images would preserve smaller structures and filter vessels but would not filter noise and plate-like structures. Therefore a hybrid technique was proposed to combine intelligently the benefits of the two techniques.

The main underlying equation in this algorithm is the anisotropic diffusion equation

${\frac {\delta u}{\delta t}}=\nabla .(D.\nabla u)$

Where

1. $\nabla$ is the divergence operator
2. $\nabla u$ is the gradient of the image $u$
3. $D$ is the diffusion tensor

The diffusion tensor $D$ allows to tune the smoothing( both the strength and direction ) across the image. $D$ is defined as a function of the structure tensor.

$J(\nabla u_{\sigma })=K_{\rho }*(\nabla u_{\sigma }\nabla u_{\sigma }^{T})$

Where $K$ is the Gaussian Kernel with standard deviation $\rho$ scale and $\nabla u_{\sigma }$ is the gradient of the image $u$ at scale $\sigma$

$D=\left[V_{1}V_{2}V_{3}\right].{\begin{bmatrix}\lambda _{1}&0&0\\0&\lambda _{2}&0\\0&0&\lambda _{3}\end{bmatrix}}.\left[V_{1}V_{2}V_{3}\right]^{T}$

$V_{1},V_{2},V_{3}$ denote the eigen vectors of the structure tensor. The eigenvalues $\lambda _{i}$ define the strength of the smoothing along the direction of the corresponding eigen vector $V_{i}$

EED, CED and Hybrid switch differ in the way they define $\lambda _{i}$.

Implementation

Using ITK framework

The implementation of this algorithm in tubetk follows ITK's finite difference solver framework. The framework uses solver filter and function. The following four main classes are implemented

1. Structure tensor generation filter ( Second moment matrix generator that is currently not available in ITK )
2. Edge-enhancing diffusion filter
3. Coherence-enhancing diffusion filter
4. Hybrid continuous diffusion filter

Current ITK diffusion filters are based on scalar or vector valued diffusion functions. But in this algorithm, we are using tensor-based diffusion function. Hence, this required doing some partial differential derivation Part 1 Part 2. This derivation is implemented in the solver function.

Class Diagram

Class hierarchy of the finite difference solver image filters

Class hierarchy of the finite difference solver function

Structure Tensor Generation Filter

Overview

The principle directions of diffusion/smoothing are based on the local structure. For this purpose, local structure tensor generator is needed. We implemented such type of filter using ITK's recursive Gaussian filter. The structure tensor $J$ is defined as follows

$J(\nabla u_{\sigma })=K_{\rho }*(\nabla u_{\sigma }\nabla u_{\sigma }^{T})$

$J(I)=K*{\begin{bmatrix}I_{x}^{2}&I_{x}I_{y}&I_{x}I_{z}\\I_{x}I_{y}&I_{y}^{2}&0\\I_{x}I_{z}&I_{y}I_{z}&I_{z}^{2}\end{bmatrix}}$

Relevant files

1. itkStructureTensorRecursiveGaussianImageFilter header file
2. itkStructureTensorRecursiveGaussianImageFilter Implementation file
3. Test

The itkStructureTensorRecursiveGaussianImageFilterTest generates primary eigen vector and eigen value images for visual inspection

tubeBasePreprocessingFiltersTests itkStructureTensorRecursiveGaussianImageFilterTest CylinderSynthetic.mha CylinderPrimaryEigenVectorImage.mha CylinderPrimaryEigenValueImage.mha

Structure tensor primary eigen vectors overlaid on the input image.

The results can be viewed best using Paraview as follows

1. Load the input image ( CylinderSynthetic.mha )
2. Apply contour filter
3. Load the primary eigen vector image ( CylinderPrimaryEigenVectorImage )
4. Apply glyph filter

You will be able to visualize the eigen vectors overlayed on the input image.

Edge Enhancing Diffusion Filter

Overview

EED is a plate enhancing diffusion filter in 3D.

Given $\lambda 1>\lambda 2>\lambda 3$, the eigen values of the structure tensor in order of decreasing magnitude, strong diffusion is performed in the direction of $V_{2}$ and $V_{3}$. The diffusion in the direction of $V_{1}$ depends on the gradient magnitude.

The diffusion tensor diagonal matrix will have the following elements

$\lambda _{e}1=1$ or $1-e^{K}$ depending on the gradient magnitude
$\lambda _{e}2=1$
$\lambda _{e}3=1$

Parameter $K$ is computed using threshold parameter, contrast parameter and gradient magnitude. If the gradient magnitude is much smaller than the contrast parameter, isotropic diffusion is performed

Relevant files

1. itkAnisotropicEdgeEnhancementDiffusionImageFilter Header file
2. itkAnisotrpicEdgeEnhancementDiffusionImageFilter Implementation file
3. Test

Input parameters

1. Contrast parameter $\lambda _{e}$ ( Default: 30 )
2. Scale parameter $\sigma$ ( Default: 1.0 )
3. Time step ( Default: 0.11 )
4. Number of iterations ( Default: 6 )

Coherence Enhancing Diffusion Filter

Overview

Given $\lambda 1>\lambda 2>\lambda 3$, the eigen values of the structure tensor in order of decreasing magnitude, CED performs diffusion in the direction of $V_{3}$ or diffusion is minimized in all three directions

The diffusion tensor diagonal matrix will have the following elements

$\lambda _{c}1=\alpha$
$\lambda _{c}2=\alpha$
$\lambda _{c}3=1$ or $\alpha +(1-\alpha ).e^{K}$ her value depending on the values of $\lambda _{2}$ or $\lambda _{3}$

Parameter $K$ is computed using the contrast parameter specified by the user.

Relevant files

1. itkAnisotropicCoherenceEnhancingDiffusionImageFilter Header file
2. itkAnisotropicCoherenceEnhancingDiffusionImageFilter Implementation file
3. Test

Input parameters

1. Contrast parameter $\lambda _{c}$ ( Default: 15.0 )
2. $\alpha$ ( Default: 0.001 )
3. Scale parameter $\sigma$ ( Default: 1.0 )
4. Time step ( Default: 0.11 )
5. Number of iterations ( Default: 6 )

Hybrid Continuous Filter

Overview

$\lambda _{h_{i}}=(1-\epsilon )*\lambda _{c_{i}}+\epsilon *\lambda _{e_{i}}$

$\epsilon$ is the EED fraction that is computed using the eigen values of the structure tensor. For detail, refer the paper.

Relevant files

1. itkAnisotropicHybridDiffusionImageFilter Header File
2. itkAnisotropicHybridDiffusionImageFilter Implementation file
3. Test

Input parameters

1. CED Contrast parameter $\lambda _{c}$ ( Default: 15.0 )
2. EED Contrast parameter $\lambda _{e}$ ( Default: 30 )
3. Hybrid Contrast parameter $\lambda _{h}$ ( Default 30.0 )
4. CED $\alpha$ ( Default: 0.001 )
5. Scale parameter $\sigma$ ( Default: 1.0 )
6. Time step ( Default: 0.11 )
7. Number of iterations ( Default: 6 )