Non-local Operational Anisotropic Diffusion Filter

TL;DR

This paper introduces a non-local weighted anisotropic diffusion filter that improves medical image quality by effectively reducing high-frequency noise in MRI and radiotherapy images, with detailed parameter analysis.

Contribution

The paper presents a novel non-local anisotropic diffusion filter that outperforms local methods and is suitable for medical imaging applications.

Findings

Enhanced image quality in MRI and radiotherapy images

Effective noise reduction demonstrated

Potential for improved radiotherapy treatment accuracy

Abstract

High-frequency noise is present in several modalities of medical images. It originates from the acquisition process and may be related to the scanner configurations, the scanned body, or to other external factors. This way, prospective filters are an important tool to improve the image quality. In this paper, we propose a non-local weighted operational anisotropic diffusion filter and evaluate its effect on magnetic resonance images and on kV/CBCT radiotherapy images. We also provide a detailed analysis of non-local parameter settings. Results show that the new filter enhances previous local implementations and has potential application in radiotherapy treatments.

Figures

Tables

Table 1. TABLE I: Maximum ADF values for constant λ 𝜆 \lambda and expression λ / | η s | 𝜆 subscript 𝜂 𝑠 \lambda/|\eta_{s}| , given the input image dimensions D 𝐷 D and the adjacency size | η s | subscript 𝜂 𝑠 |\eta_{s}| .

$𝑫$	$\left|{𝜼}_{𝒔}\right|$	$𝝀$	$𝝀/\left|{𝜼}_{𝒔}\right|$
2	4	4/5	1/5
2	8	8/7	1/7
3	6	6/7	1/7
3	18	18/13	1/13
3	26	78/47	3/47

Table 2. TABLE II: Quantitative results of application of WADF over BrainWeb Phantom MRI dataset.

$\gamma R$	IQI	MSE	PSNR
16${\gamma }^{t-1}$	0.9898	3,742	27.59
32${\gamma }^{t-1}$	0.9945	2,986	28.77
48${\gamma }^{t-1}$	0.9964	2,836	29.09
64${\gamma }^{t-1}$	0.9970	2,788	29.21
80${\gamma }^{t-1}$	0.9949	2,777	29.26
96${\gamma }^{t-1}$	0.9942	2,802	29.22

Table 3. TABLE III: Quantitative results of application of NL-WADF over BrainWeb Phantom MRI dataset.

SR; PR; PD; #R; $\gamma R$	IQI	MSE	PSNR
2.0; 1.1; 0.5; 2; 0.32	0.9968	2,719	29.46
3.0; 1.1; 0.5; 2; 0.32	0.9967	2,712	29.48
4.0; 1.1; 0.5; 2; 0.32	0.9967	2,710	29.49
4.0; 1.9; 1.0; 2; 0.32	0.9941	3,022	28.71
4.0; 1.1; 0.5; 2; 0.48	0.9956	2,749	29.45
2.0; 1.9; 2.0; 1; 0.64	0.9972	2,793	29.20
3.0; 1.9; 0.5; 1; 0.64	0.9972	2,776	29.23
4.0; 1.9; 1.0; 2; 0.64	0.9972	2,788	29.20

Table 4. TABLE IV: Qualitative evaluation of planar kV and 3D CBCT images by physicist specialist. Grades from 1 (poor) to 5 (excellent) were given to each image. Mean values are presented here.

${\gamma }^0$	Score
0.2${\sigma }_{G_s}$	4.9
0.4${\sigma }_{G_s}$	5.0
0.6${\sigma }_{G_s}$	4.9
0.8${\sigma }_{G_s}$	4.9
${\sigma }_{G_s}$	4.8

Equations

$$I_s^{t+1}\approx I_s^t+\frac{\lambda }{|{\eta }_s|}\sum _{p\in {\eta }_s}g(\nabla I_{s,p}^t,{\gamma }^t)\nabla I_{s,p}^t$$

$$g(\nabla I_{s,p}^{t},\gamma^{t})=\left\{\begin{array}[]{l}\left[1-\left(\frac{\nabla I_{s,p}^{t}}{\sqrt{5}\gamma^{t}}\right)^{2}\right]^{2},\nabla I_{s,p}^{t}\leq\gamma^{t}\sqrt{5}\\ 0,\text{otherwise.}\end{array}\color[rgb]{1,1,1}\right\}$$

$$G_s=\sum _{p\in {\eta }_s}{I}_s-I_p$$

$$I_s^{t+1}\approx I_s^t+\frac{1}{|{\eta }_s|}\sum _{p\in {\eta }_s}\frac{g(\nabla I_{s,p}^t,{\gamma }^t)}{d_{s,t}}\nabla I_{s,p}^t$$

$$I_s^{t+1}\approx I_s^t+\frac{1}{|H_s|}\sum _{p\in {\eta }_s}\frac{g(\nabla I_{s,p}^t,{\gamma }^t)}{d_{s,t}}\nabla I_{s,p}^t$$

$$SSIM(I,J)=\frac{(2{\mu }_I{\mu }_J+c_1)(2{\sigma }_{IJ}+c_2)}{({\mu }_I^2+{\mu }_J^2+c_1)({\sigma }_I^2+{\sigma }_J^2+c_2)}$$

$$c_1=0.0001|I||J|$$

$$c_2=0.0009|I||J|$$

$$MSE(I,J)=\frac{1}{|I|}\sum _s(I_s-J_s{)}^2$$

$$PSNR(I,J)=20\ast \log (\max (I^M,J^M))-10\log (MSE(I,J))$$

$$IQI(I,J)=\frac{(4{\mu }_I{\mu }_J{\sigma }_{IJ})}{({\mu }_I^2+{\mu }_J^2)({\sigma }_I^2+{\sigma }_J^2)}$$

Full text

I Introduction

Medical imaging is subject to several kinds of noise and artifacts coming from the acquisition and storage procedure [1]. The nature and appearance of the noise depend on distinct factors such as the imaging modality, the scanned body part, scanner parameters and positioning, among others [2, 3].

Despite recent improvements to the acquisition process, high-frequency filters are essential to improve medical image quality. Nevertheless, their application is not a trivial task, since some of the original signals may be removed and new artifacts might be inserted to the image while removing the noise [4, 5].

There are three generations of high-frequency noise filters: the isotropic which were proposed since the beginning of digital medical imaging [6, 7]; the anisotropic – which became popular in the 90s [8] – are still in use for its robustness to generating artifacts, and is present in some popular medical imaging tools such as FLS [9] and 3D Slicer [10]; and the non-local anisotropic filters which started to be developed around 2005 with the Non-Local Means (NLM) algorithm [11]. Even though the third generation produces better results in terms of several metrics, in the context of medical imaging, one should be even more careful since they tend to insert more artifacts to the filtered image [4].

In this paper we focus on the implementation and testing of an operational third generation filter which is an extension of the anisotropic diffusion filter (ADF) [12] for reducing addictive noise present in MRI and CT images. The proposed implementation is operational with respect to local and non-local conservativeness parameter. We also conduct an analysis of the non-local parameters in the models proposed in [13]. Our experiments include synthetic MRI images to estimate parameter values and also a qualitative evaluation of kV/CBCT radiotherapy images by specialists.

The reminder of this paper is organized as follows: Section II presents a review of related works in high-frequency image segmentation; Section III contains an explanation of ADF and its extensions; Section IV proposes our new filtering method; Section V shows our quantitative and qualitative experiments; and Section VI states our conclusions.

II Related Works

The first generation of high-frequency noise filters are isotropic, that is, the same operation is applied to each image pixel despite its contents or position. Some of the traditional isotropic filters are the arithmetic and geometrical means, median, and alpha-trimmed mean filter [7].

The second generation of high-frequency noise filters is applied over other transform domains and/or with local anisotropic functions. These filters apply distinct operations according to pixel features such as intensity frequency or local patterns. Wavelet transform and cosine transform based [14, 15], bilateral and trilateral [16, 17], and the diffusion [18, 19] are among the most popular filters. While first generation filters were simpleton implementations, the second generation filters bear several non-trivial parameters which yielded lots of tuning and evaluation papers [2].

Third generation filters are in most part extensions of the existing second-generation implemented to take advantage of non-local data. The basic idea is to use similar patches in order to estimate and remove high-frequency noise. Examples of bilateral filter extensions are the NLM [11], BM3D [20]. Third generation filters have also been applied to medical imaging in the context of MRI and CT [3].

In the literature, there is not much effort in developing a non-local ADF filter for MRI and CT images [2, 3]. Only two third generation of ADF were proposed in [21, 22]. These extensions are just a trivial addition of non-local patches to the algorithm with no mathematical or empirical selection of optimum parameters.

In [21], the authors use the original edge-stopping functions. In [22], the Tukey’s biweight robust estimator was used. Either of them estimated similar patches based on pre-filtered images. The second implementation computed adjacent patches at each iteration of the filter which is unnecessary and computationally expensive.

Therefore, there is a lack of studies of a non-local ADF which is the main subject of this paper.

III Anisotropic Diffusion Filter

The ADF was proposed by Perona and Malik [12]. It consists in applying concepts of fluid thermodynamics to filter images. Just as molecules with different temperature from its adjacency tends to exchange heat until achieving homogeneous temperature, there would have a propensity to match high-frequency noisy pixel intensities with their local neighbours. In this scenario, image edges are similar to tow large adjacent volumes of fluid with higher capacitance, demanding more time to homogenize the temperature in comparison to the noisy smaller volumes. Therefore, ADF is capable of eliminating the high-frequency noise while preserving image strong edges.

The discrete implementation of ADF is given by Equation 1:

[TABLE]

$I^{t}_{s}$ stands for the intensity of pixel s in image $I$ at the instant $t$, $\lambda$ is a diffusion rate scalar, $\gamma^{t}$ is a positive variable related to the diffusion or smoothing strength which decreases at each iteration, $\eta_{s}$ represents the set of pixels adjacent to $s$, $g(\cdot)$ is an edge stopping function (ESF) which controls the diffusion process, and $\nabla I^{t}_{s,p}$ is the image directional gradient from $s$ to $p$ at instant $t$. The directional gradient $\nabla I^{t}_{s,p}$ can be approximated by $I^{t}_{p}-I^{t}_{s}$.

The diffusion process follows the ESF which has been widely studied and tested [19, 23]. Tukey’s biweight robust estimator ESF is given in Equation 2 and is an excellent option since it does not affect edges with intensity variation over a certain threshold [13].

[TABLE]

With respect to parameter $\lambda$, the authors in [18] proposed the maximum values in order to keep a monotonic variation of the intensities at each iteration, so that no artifacts are inserted to the image. In [13], the authors corrected one of the constants for 26-adjacency. Table I shows the expected maximum values of $\lambda$ and of the fraction $\lambda/|\eta_{s}|$ according to the number of dimensions of the image ($D$) and to the adjacency size ($|\eta_{s}|$).

There were several different solutions proposed to estimate the initial $\gamma^{0}$ for the first iteration and the number of iterations $T$. Some approaches used global and/or local gradient to compute $\gamma$ [12, 19, 24] while others use a planar region (i.e. without strong edges) [25].

The approach in [13] used an adaptive method based on both the global gradient and a planar region in order to estimate the best initial $\gamma^{0}$. While the gradient indicates an upper bound $\gamma{E}$ to avoid filtering edges, the planar region provides a lower bound $\gamma{F}$ so that practically all noise is removed. In case of low signal to noise ratio (SNR)111In the experiments, Gaussian or Rician noise of around 3% of amplitude is the maximum SNR which is perfectly filtered without blurring significant edges., it may happen that $\gamma{F}>\gamma{E}$. In these cases, a value between $\gamma{E}$ and $\gamma{F}$ should be chosen. Lower $\gamma^{0}$ translates into a conservative ADF which preserves the stronger edges, while higher $\gamma^{0}$ constitutes into an aggressive filtering, removing all noise.

Also, since its proposal [12] and as corroborated by [24], it is very important to update $\gamma^{t}$ at each iteration. As showed in [13], it is fundamental to reduce $\gamma^{t}$ at each iteration in order to preserve the borders and it may be very expensive to recompute $\gamma^{t}$ at each iteration based on the noise in the planar region and on the strong edges. A conservative proposed solution was to set $\gamma^{t}=0.25\gamma^{t-1}$. That is because it is expected that the strongest noise has at most one-fourth of its adjacent pixels with a similar intensity, while for the edges, it is reasonable to expect that at least one-fourth of the pixels have distinct intensity. Therefore, reducing $\gamma$ by $0.16\gamma$ was a conservative adopted criteria222For $\gamma$ reduction values of other ESFs, please refer to [13]..

Finally, with respect to the number of iterations $T$, in [25] the authors verify the intensity of the strongest edges. That solution is not suitable for decreasing $\gamma^{t}$ with Tukey’s biweight robust estimator given in Equation 2. In [13], the authors propose to stop as $\gamma^{t}<\gamma_{\cal E}/7$ which does not depend on the noise intensity itself.

IV Non-Local Optimal Anisotropic Diffusion Filter

The non-local ADF (NL-ADF) proposed in [21, 22] do not optimize any filter parameter. We propose a novel non-local ADF with enhanced parameter optimization. It has improvements over the previous optimization procedure proposed in [13]. We also study the non-local patches influence in the filtering process.

IV-A Local Parameter Optimization

Instead of estimating the initial $\gamma^{0}$ based on the standard deviation of edgy or noisy regions, in this paper we propose a more robust and straightforward strategy. Noise pixels in flat regions have a distinct intensity value than all their adjacents. An edge pixel, on the other hand, has a similar intensity to pixels in the same side of the edge and distinct intensity from pixels in the opposite side of the edge. Therefore, instead of computing the standard deviation of pixels to estimate noise, we propose to use the sum of the directional gradients $G_{s}$ as in Equation 3:

[TABLE]

Note that the direction of $G_{s}$ is always computed toward the source pixel $s$. This way, edge pixels will have opposite edge pixel intensities neutralized, which noise pixels will display the highest values of $G_{s}$. We then select standard deviation ($\sigma_{G_{s}}$) of the 5$\%$ highest $G_{s}$ pixels as our $\gamma^{0}$, as the maximum effect of of the ADF with Tukey’s biweight robust estimator happens for $\gamma^{0}=\sigma_{G_{s}}$.

A second contribution is to revise the $\gamma^{t}$ reduction strategy proposed in [13] based on 1/4 of the adjacent pixels as explained in Section III as it was a very conservative choice. In this paper, we evaluated faster and slower $\gamma^{t}$ reductions based on values in range $[0.16\gamma^{t-1},0.96\gamma^{t-1}]$.

With respect to the stopping criteria a simple alternative is to stop filtering the image as $\gamma^{t}\leq 0.01I^{M}$ ($I^{M}$ is the maximum intensity of an image $I$), as 1$\%$ of maximum intensity is imperceptible by humans eyes in medical images.

Finally, diffusion of adjacent pixels should also consider the distance between the source and target pixels. Therefore, we propose the weighted ADF (WADF) based on Equation 4:

[TABLE]

where $d(s,t)$ is the distance between pixels $s$ and $t$. Note that, according to the study proposed in [18, 13] with the WADF for a monotonic variation of the intensities at each iteration we may set $\lambda=1$, eliminating this parameter from the original ADF Equation 1.

IV-B Non-Local Anisotropic Diffusion Function

An interesting question overlooked in previous publications is the weight or distance of non-local patches. Most of the papers, including the previous NL-ADF [21, 22], position the non-local patches in a discrete way in an additional third or fourth image dimension. That is, the closest patch is one-pixel-wide distant from the patch of interest. In a similar way, the original NLM [11] uses the spatial distance as a weighting factor, but this is an arbitrary solution. BM3D [20] and BM4D [26] also build a structure with an additional discrete dimension.

As the ADF does not depend on the distance of patches itself, we decided to evaluate patch distance and quantity over the filtering result. The $n$ more similar patches will have the same distance $d$ to the filtered pixel $s$ as a third (or fourth) dimension to a 2D (or 3D) image. Note that $d$ may assume any float value. We will evaluate multiple patch distance variation in a future work.

The choice of the distance $d$ also impacts on the computation of constant $\lambda$. Still, considering the WADF, we may set $\lambda=1$ even for the non-local weighted anisotropic diffusion filter (NL-WADF) in Equation 5:

[TABLE]

$H_{s}$ comprehends all adjacent pixels to $s$ both local and non-local.

In the case of NL-WADF, the number of pixels with similar intensity to noise and edge pixels should stay proportionally the same. Therefore, we propose also to evaluate the $\gamma^{t}$ decrease in every iteration using the same range of $[0.32\gamma^{t-1},0.96\gamma^{t-1}]$.

Also, we determine the similar patches based on filtered images by a median filter [27]. After this, the patches are used by NL-WADF over the original image.

V Experiments

The first experiment involves estimating the best parameter values for local and non-local WADF. For this purpose, we used BrainWeb Phantom dataset [28] changing the noise intensity between 1 and 9$\%$ of the image maximum intensity. Figure 1 contains slices of the ground-truth image and of an image corrupted by a noise of 7$\%$ of the maximum intensity. The second experiment is a qualitative evaluation in the context of radiotherapy imaging as shown in Figure 2. Images were acquired with Varian TrueBeam and Varian Clinac IX Radiotherapy Systems.

V-A Quantitative Evaluation

We used the following metrics to evaluate the quality of the filtered image in comparison with the ground-truth image with no noise: structural similarity index(SSIM) Equation 6, mean square error(MSE) Equation 7, peak signal-to-noise ratio(PSNR) Equation 8, and image quality index(IQI) Equation 9. We will not present here results for SSIM, since all images achieved accuracy equal to $1.0$. Lower scores are better in terms of MSE and higher scores are better with respect to SSIM and PSNR.

[TABLE]

[TABLE]

[TABLE]

$I$ and $J$ are the filtered and ground-truth images, respectively, $\mu_{I}$ and $\sigma_{I}$ are the mean and standard deviation of pixel intensity in $I$, respectively, and $|I|$ is number of pixels in $I$. $\sigma_{I}J$ is the covariance between images $I$ and $J$.

[TABLE]

[TABLE]

$\max(x,y)$ is the maximum value between x and y. $I^{M}$ is the maximum pixel intensity in $I$.

[TABLE]

For the WADF, we tested the $\gamma^{t}$ reduction rate of $\gamma R=\{0.16\gamma^{t-1},0.32\gamma^{t-1},0.48\gamma^{t-1},0.64\gamma^{t-1},0.80\gamma^{t-1},0.96\gamma^{t-1}\}$. We also evaluated a conservativeness parameter which sets the initial $\gamma_{0}=0.2\sigma_{G_{s}},0.4\sigma_{G_{s}},0.6\sigma_{G_{s}},0.8\sigma_{G_{s}},and\sigma_{G_{s}}$.

We notice in Table II that the best parameter set for WADF with respect to $\gamma^{t}$ reduction rate in average to all noise levels is between 0.64$\gamma^{t-1}$ and 0.80$\gamma^{t-1}$ which is much higher than the value adopted in [13]. It is important to mention that the accuracy given by the proposed metrics are not always favorable to 0.64$\gamma^{t-1}$ and 0.80$\gamma^{t-1}$ reduction, but 0.32$\gamma^{t-1}$ provides better results for more conservative filtering(i.e. $\gamma_{0}\leq 0.4\sigma_{G_{s}}$).

For the NL-WADF, we tested the following parameters: $\gamma R=\{0.32\gamma^{t-1},0.48\gamma^{t-1},0.64\gamma^{t-1},0.80\gamma^{t-1},0.96\gamma^{t-1}\}$; $\gamma_{0}=0.2\sigma_{G_{s}},0.4\sigma_{G_{s}},0.6\sigma_{G_{s}},0.8\sigma_{G_{s}},and\sigma_{G_{s}}$; search radius (SR = 2.0, 3.0 4.0 pixels), non-local patch radius (PR = 1.1, 1.5, 1.8 pixels), the distance from the central pixel of non-local patches (PD = 0.5, 1.0, 2.0 pixels), and the number of non-local patches ($\#$P = 1, 2).

Figure 3 shows sample slices of image showed in Figure 1(b) filtered conservatively($\gamma^{0}=0.4\sigma_{G_{s}}$) and aggressively($\gamma^{0}=0.8\sigma_{G_{s}}$) with parameters SR=3.0, PR=1.9, PD=1.0, $\#$P=2, and $\gamma$R=0.80$\gamma^{t-1}$.

Table III contains the best results for NL-WADF in terms of the proposed metrics. The best set of parameters in average according to MSE and PSNR is to use SR=4.0, PR=1.1, PD=0.5, $\#$P=2, and a $\gamma$R= 0.32$\gamma^{t-1}$. Note that $\gamma_{0}$ is not related to the quality of the result, but to the desired conservativeness, which is more related to the goal application. According to IQI, the best set of parameters is SR=3.0, PR=1.9, PD=0.5, $\#$P=1, and $\gamma$R=0.64$\%$.

Note though, that this is by no means the best parameter set for all noise levels and desired conservativeness. Based on the results in Tables II and III we can also notice that NL-WADF may generate better results as compared to WADF, but this improvement is not statistically significant.

Other important observations are: using one non-local patch produces in 65.7$\%$ of cases a better accuracy than employing two of them and in 69.9$\%$ of the times for more aggressive filtering (i.e. $\gamma_{0}\geq 0.6$). In 59.1$\%$ of the times, for a noise level of 7$\%$, using two non-local patches produces a better accuracy. That is the only noise level with better results using 2 non-local patches.

In terms of IQI, $\gamma^{t}$ reduction of 0.96$\gamma^{t-1}$ produces the higher accuracy in 99$\%$ of the tests, except while dealing with a noise level of 7$\%$.

The analysis of the size of the non-local patch radius is very similar to the number of patches. A radius of 1.1 pixels generates higher accuracy in 52.3$\%$ of times. Smaller radius are a better choice in 62.2$\%$ of times for aggressive filters (i.e $\gamma_{0}=0.6\sigma_{G_{s}}$).

Setting non-local patches at 0.5 pixels from the source patch generates the best results in 55$\%$ of the times and in 61.8$\%$ of the times for more aggressive filtering (i.e. $\gamma_{0}\geq 0.6$). This means that non-local patches should have a stronger impact in removing the noise.

With respect to the size of the similar patches search region, a radius of 2.0 pixels is usually sufficient and produces the best results. Conservative parameters are an exception (i.e. $\gamma_{0}\leq 0.4$) over high-intensity noise in which a radius of 4.0 pixels produces better results. This indicates that a larger search radius is an important parameter to conserve image edges.

V-B Qualitative Evaluation

Using the best parameter set for NL-WADF, we performed a qualitative test with cone beam computed tomography (CBCT) and kV planar radiotherapy images of the low pelvis region provided by hospital AC Camargo in São Paulo/SP, Brazil. A physicist specialist evaluated the quality of the filtered images and their utility for image visualization purposes.

Figures 4 and 5 shows sample slices of filtered planar kV and 3D CBCT images with conservative and aggressive NL-WADF.

Five filtering levels for conservative to aggressive of ten images (3 planar and 7 3D CBCT volumes) were evaluated. The specialist assigned a score from 1 (poor) to 5 (excellent) to each of the filtered images. Table IV show the scores. As noticed for radiotherapy visualization purposes, the best filter consists of applying $\gamma^{0}=0.4\sigma_{G_{s}}$. More conservative parameters do not remove significantly the noise. Using $\gamma^{0}>0.4\sigma_{G_{s}}$ on the other hand on more aggressive filter removes important details in the opinion of the specialist.

VI Conclusions

In this paper, we propose an operational non-local anisotropic diffusion filter for additive noises present in medical images such as MRI and CT. The main contributions of the filter were the distance factor never considered before for non-local patches filtering for ADFs and the weighted contribution of adjacent pixels. Also, we present an extensive evaluation of the parameter values in order to define the best parameter set in terms of SSIM, PSNR, MSE, and DE. Results show that the new non-local version of the ADF is significantly better than the local filter. Finally, we also present a qualitative evaluation of the filter in the context of visual quality of radiotherapy imaging which showed that filtering the image is an important procedure, but the procedure must be conservative avoiding removing important anatomical details.

Future works involve proposing more sophisticated strategies for setting the distance between non-local adjacent patches depending on the source pixel neighbourhood and also based on the neighbourhood similarity.

Acknowledgment

We would like to thank AC Camargo and Dr Luiz Juliano Neto for the support with the radiotherapy images. We also thank the training provided by Varian Medical Systems in agreement with the Brazilian Ministry of Health.