A spatially variant high-order variational model for Rician noise removal

Rician noise removal is an important problem in magnetic resonance (MR) imaging. Among the existing approaches, the variational method is an essential mathematical technique for Rician noise reduction. The previous variational methods mainly employ the total variation (TV) regularizer, which is a first-order term. Although the TV regularizer is able to remove noise while preserving object edges, it suffers the staircase effect. Besides, the adaptability has received little research attention. To this end, we propose a spatially variant high-order variational model (SVHOVM) for Rician noise reduction. We introduce a spatially variant TV regularizer, which can adjust the smoothing strength for each pixel depending on its characteristics. Furthermore, SVHOVM utilizes the bounded Hessian (BH) regularizer to diminish the staircase effect generated by the TV term. We develop a split Bregman algorithm to solve the proposed minimization problem. Extensive experiments are performed to demonstrate the superiority of SVHOVM over some existing variational models for Rician noise removal.


INTRODUCTION
Magnetic resonance (MR) images have been widely used in medical imaging.Due to the thermal noise caused by patients during the scan process (Nowak, 1999;Aja-Fernández & Vegas-Sánchez-Ferrero, 2016), the MR images are inevitably degraded.The noises in the MR images have negative impacts on various tasks of medical image processing and analysis, such as classification, segmentation, visualization (Aja-Fernández & Vegas-Sánchez-Ferrero, 2016).Hence, noise removal is the fundamental task for processing MR images.
It was shown that the noises in the MR images can be modeled by the Rician distribution (Henkelman, 1985;Aja-Fernández & Vegas-Sánchez-Ferrero, 2016).The task of Rician noise removal refers to estimating the clean MR image from a noisy one.Since the Rician noise is signal-dependent, it is a great challenge to denoise the clean MR image.Next, the previous works that address the problem of Rician noise removal are reviewed.First, several denoising methods based on statistics are presented.In (Henkelman, 1985;McGibney & Smith, 1993;Bernstein, Thomasson & Perman, 1989), the first and second moments of the Rician distribution were employed to estimate the clean MR images.Using the local sample statistics, Aja-Fernández & Vegas-Sánchez-Ferrero (2016) derived a closed-form staircase effect generated by TV.An efficient split Bregman algorithm is developed to solve the proposed model.The proposed model is evaluated on a large dataset in comparison with several existing variational models for Rician noise removal.Experimental results show the superiority of SVHOVM in Rician denoising.
The main contributions include the following: • A novel variational model for Rician noise removal is proposed.Particularly, the common TV regularizer is modified such that it becomes spatially variant according to the characteristics of pixels.The BH regularizer, which is a high-order term, is utilized to enhance the denoising results.
• An efficient split Bregman algorithm is developed to solve the proposed problem.
• Extensive experiments are conducted to discuss the effects of the parameters of SVHOVM and to evaluate its performance.Experimental results show that the proposed model outperforms some existing variational models for Rician noise reduction.
The rest of the article is organized as follows.In Section 'Preliminary and related works', some preliminaries and a brief overview of related works are presented.The proposed model is described in Section 'Proposed model'.In Section 'Numerical implementation', a split Bregman algorithm for solving the proposed problem is presented.Section 'Experimental results' discusses experimental results.

PRELIMINARY AND RELATED WORKS Preliminary
The MR imaging systems use quadrature detectors to produce two-or three-dimensional complex data.The raw MR data is always perturbed by Gaussian noise.The complex representation of the raw MR data is given by where F R and F I are the real and imaginary parts of the raw MR data F; u ∈ R p×q is the true amplitude of the noise-free image; η 1 and η 2 ∈ R p×q are Gaussian noise with zero mean and standard deviation σ .For clinical analysis, the magnitude MR images are often used.Mathematically, the magnitude MR image is computed by Since the magnitude MR images are obtained by the non-linear transformation, the distribution of the overall noises for the magnitude MR image is no longer Gaussian.It was shown in Henkelman (1985), Aja-Fernández, Alberola-López & Westin (2008) that the noises in the magnitude MR images have the Rician distribution which is given by where I 0 is the modified Bessel function of the first kind with order zero (Bowman, 2012).
The form of the modified Bessel function of the first kind with real order ν are given by with (n) = (n − 1)! is the gamma function; ν ∈ R.

Related works
The goal of Rician noise removal is to estimate the noise-free image u from the noisy magnitude MR image f .Most of variational models utilize the MAP approach to estimate u by maximizing a posterior given f , that is ũ = arg max u P(u|f ).In Getreuer, Tong & Vese (2011), they applied the Bayes's rule to derive the MAP model as where the first two terms form the data fidelity, which is derived from (3) using the MAP framework; the last term is the TV of u; α is a non-negative regularization parameter; ∇ is the gradient operator; is the image domain.
Since the MAP model is non-convex, Getreuer, Tong & Vese (2011) approximated its data fidelity term by a convex function as follows where with A(•) is the cubic rational polynomial approximation of I 1 (•)/I 0 (•) with I 1 (•) being the modified Bessel function of the first kind with first order; c = 0.8426.By exploring the statistical property of the Rician distribution, Chen & Zeng (2015) added the quadratic term into the MAP model to obtain a strictly convex model as Following the similar idea, Yuan (2018) added a convex gradient data fidelity into the MAP model as For brevity, the following notations are used: GTV for the convexified model derived by Getreuer, Tong & Vese (2011) (Eqs.( 5)-( 8)); CZ for the model of Chen & Zeng (2015) (Eq.( 9)); Yuan (2018) for the model proposed by Yuan (Eq.( 10)).

PROPOSED MODEL
As described in 'Introduction', existing variational methods for Rician noise removal mainly utilize TV.Besides, the adaptability received less research attention.To this end, the author proposes a spatially variant high-order variational model (SVHOVM) for Rician denoising.The spatially variant TV (SVTV) and the bounded Hessian (BH) regularizers are introduced by minimizing the following minimization functional where < •,• > denotes the Euclidean inner product; • 1 and • 2 stand for the 1 -and 2 -norms, respectively; α(f )∇u 1 is the SVTV regularizer with α(•) being the weighting function; ∇ 2 u 1 is the BH regularizer where ∇ 2 denotes the Hessian operators; β are non-negative regularization parameters.
The BH regularizer is exploited to remove the side effect produced by the TV term.In Papafitsoros & Schönlieb (2014) showed that the BH regularizer is able to remedy the staircase effect and to preserve structural details.The SVTV term is the common TV regularizer weighted by the function α(•) for each pixel.The weighting function is defined as where α 0 is a non-negative parameter; G ω stands for the Gaussian filter with zero mean and standard deviation ω; κ denotes a contrast parameter; '' * '' represents the convolution operator.
The weighting function can adaptively manipulate the smoothing strength of the TV regularizer.Its values vary depending on the image gradients of pixels.Namely, for a fixed κ, in flat regions where |∇G ω * f | < κ, the weighting function is large, which means a strong noise reduction.In contrast, at object edges where |∇G ω * f | > κ, the weighting function is small, which indicates the edge preservation.Thus, the weighting function is effective in reducing noise while maintaining object edges.

NUMERICAL IMPLEMENTATION
In this section, a split Bregman algorithm is developed to solve the proposed problem (11).Following (Goldstein & Osher, 2009), two auxiliary variables are introduced to obtain the following constrained problem: where d and z are auxiliary variables.Applying the Bregman iteration gives the following unconstrained problem: where b 1 and b 2 are the Bregman iteration variables; θ 1 and θ 2 are the penalty parameters.The problem Eq. ( 14) is solved by an alternating direction method (Gabay & Mercier, 1976;Bertsekas, 2014).In each step, either u, d or z is minimized while keeping other variables fixed.With u and z fixed, the d-subproblem is obtained as: which has the following solution: Similarly, the z-subproblem and its solution are given by By fixing d and z, the u-subproblem is obtained as Let E(u) denote the functional of ( 19).The Newton's method is applied to solve the u-subproblem (19) as follows: with where the variable (uf /σ 2 ) of I 0 and I 1 is omitted for brevity.Finally, the b 1 and b 2 variables are updated by: In summary, the denoised image is found by iteratively computing u, d and z via the sequence of Eqs. ( 16), ( 18), ( 20), (21), and ( 22).The number of iterations is used as the stopping criterion.It is worthy of note that the weighting function is iteratively refined by computing (12) using the restored image of the previous iteration.This refinement offers an enhanced weighting function, resulting in better denoised images.The overall split Bregman algorithm for solving the proposed problem (11) is summarized in Fig. 1.

EXPERIMENTAL RESULTS
In this section, numerical experiments are conducted to discuss the affects of the parameters of SVHOVM and to evaluate the proposed model in comparison with existing variational methods for Rician noise removal.The experiments are performed on the IXI and SB datasets.The IXI dataset (Information eXtraction from Images: https://brain-development.org/ixi-dataset/.)contains real MR images.The SB dataset contains simulated MR images generated by BrainWeb (Simulated Brain Database: http://www.bic.mni.mcgill.ca/brainweb/.)(Cocosco, 1997;Kwan, Evans & Pike, 1999;Kwan, Evans & Pike, 1996;Collins et al., 1998).It consists of three-dimensional T1w, T2w, and PDw volumes of 181×217×181 voxels with zero noise.The MR images are perturbed by Rician noise with three noise levels σ = 5,15 and 25.Sample MR images of the datasets are shown in Fig. 2. The PSNR and SSIM (Wang et al., 2004) are used to measure the performance of models.Besides, visual quality is also employed for qualitative evaluation.Let u and ũ denote the noise-free and the denoised images.The PSNR and SSIM indices are defined as: where M and N are the sizes of images; µ ũ, σ ũ and µ u , σ u are the means and standard deviations of ũ and u, respectively; c 1 and c 2 are constants.

Ablation study
In this section, numerical experiments are conducted to discuss the influence of various parameters on the proposed algorithm (Fig. 1).The regularization parameters α and β, the contrast parameter κ, and the inner iteration number n max are considered.Note that the inner iteration corresponds to the loop of the Newton's method (Eqs.20-22).

Regularization parameters
The effects of the parameters α 0 and β on the performance of SVHOVM are investigated.These parameters determine the weights of the SVTV and BH regularizers.One of these parameters is varied while keeping the other fixed.The effects of α 0 are shown in Figs.3C-3E.The parameter β is fixed to a low value (particularly, β = 1) in order to diminish the influence of this parameter.The parameter α 0 is set to 1,10 and 20 to show its effects.One can see that the parameter α 0 controls the smoothness of denoising results.As α 0 increases, more noise is reduced but the staircase effect becomes more obvious.Then, the influence of the parameter β is examined.The parameter α 0 is fixed by a large value (particularly, α 0 = 15) in order to generate the staircase effect.The parameter β is gradually increased to demonstrate its effects.Figures 3F-3H show that as β gets larger, the BH regularizer diminishes the staircase effect more effectively, producing smooth transition between flat regions.However, the large values of β result in blurred images.Thus, the regularization parameter β should be not too large so that the artifacts generated by the TV term are diminished without generating any serious blur in the denoised images.
Next, the dependence of the PSNR measure on the parameters α 0 and β is considered.The parameter α 0 is fixed at different values while β is varied.Figure 4A shows that the PSNR values change in a continuous manner.For a fixed α 0 , the PSNR measures initially increase with the value of β, reaching the maximum value and then decreasing.When α 0 increases, the PSNR measure rises to the global maximum, followed by a decrease.One can see from Fig. 4A that an optimal denoised result can be attained by alternatively adjusting the two regularization parameters α 0 and β.

Contrast parameter
Figure 4B shows the dependence of the weighting function α(•) on the gradient magnitude for different values of the contrast parameter κ.As can be seen, the weighting function is monotonically decreasing with the increase of the gradient magnitude.As the gradient magnitude becomes larger, α(•) goes to 0 and the strength of the TV term is down-weighted.Thus, the weighting function controls the regularization strength of the TV term.Figure 4B also demonstrates that the parameter κ adjusts the range in which the weighting function receives low values.As κ declines, this range is extended.It means that when κ decreases, the more image details as well as noise are preserved.Therefore, an optimal value of κ needs to be determined to achieve a balance between noise reduction and image detail preservation.The parameter κ is set to 0.8 empirically.

Inner iteration number
The impacts of the inner iteration number n max on the performance of the proposed algorithm are examined.Figure 5A shows that the proposed algorithm is less sensitive to the number of inner iterations.Meanwhiles, the computational time per iteration of the proposed algorithm increases by about 25% on average when the number of inner loops is increased by one (Fig. 5B).Thus, one inner iteration is used in order to reduce the computational complexity of the proposed algorithm without affecting the quality of the denoised images.

Comparative study
Next, SVHOVM is compared with some existing variational models for Rician noise removal.The models GTV, CZ, and Yuan are used as references.The regularization parameters of models are tuned using a simple alternating optimization method to achieve  Glowinski, Pan & Tai (2016), the penalty parameters θ 1 and θ 2 are set by 5 and 5, respectively.The number of iterations is fixed to 500.Beyond this, the performance of models is nearly unchanged.
Table 1 shows the PSNR and SSIM results of different models on the IXI and SB datasets.The figures highlighted in bold represent the best results for each volume and noise level.From Table 1, the following observations are made.First, the proposed model attains the best results for most of the cases.SVHOVM fails to achieve the best performance in SSIM for the T1w volume of the SB dataset with σ = 15 and 25.On average, SVHOVM outperforms GTV, CZ, and Yuan models by 0.76,1.39,0.95 in PSNR and 0.0232,0.0561,0.00346 in SSIM, respectively.Second, as the noise level increases, the improvement of SVHOVM over competitive models declines.The proposed model gives the average gains in PSNR of 1.23,0.99,and 0.87 for σ = 5,15, and 25, respectively; the corresponding figures in SSIM are 0.0349,0.0233,and 0.0561.It can be explained that large noises reduce the effectiveness of the weighting function.Third, SVHOVM yields the best performance on the T1w volume, followed by the PDw volume and then the T2w volume.
The advantages of SVHOVM are further confirmed by Figs.6-8. Figure 6 shows denoising results of different models on an image of the IXI dataset with σ = 5.The parts of the denoised images are enlarged for visual comparison.Besides, the curves of 1D intensity values are shown.One can see from Fig. 6 that SVHOVM yields the best denoised image.For the results of SVHOVM, the flat regions are smoother and the object edges are sharper compared with those of other models.The intensity curve of SVHOVM fits to the ground truth better than those of the competitive models.
Figures 7 and 8 demonstrate the restored images of different models on the images of the IXI dataset for the higher noise levels.The residual images, which are the image difference between the noisy images and the denoised images, are shown.It can be seen that the structural information exists in the residual images of all models.It is due to the fact that the Rician noise is signal-dependent.One can observe that fewer information is left in the residual images of SVHOVM over the cases of the competitive models.In summary, the quantitative and qualitative results show superior performance of SVHOVM compared with other existing variational models for Rician noise removal.

CONCLUSION
In this article, the author presented a spatially variant high-order variational model for Rician noise removal.The SVTV regularizer was proposed in order to adjust the smoothing strength according to the characteristics of pixels.In addition, the proposed model employs the BH regularizer to reduce the staircase effect.The split Bregman algorithm was derived to solve the minimization problem efficiently.Extensive numerical experiments showed that the proposed model outperforms the existing variational models in terms of both quantitative and qualitative criteria.
The author hopes that the proposed method can serve as a tool for clinical analysis.One limitation of the proposed method is that it depends on parameters, especially the regularization parameters.On the one hand, the parameters allow the clinical experts to adjust the level of noise reduction to observe the image details.On the other hand, a parameter-dependent method requires the users to understand the effects of the parameters in order to obtain the optimal results.In the future, the author will investigate a method to automatically search for the optimal parameters of the proposed model.

Figure 1
Figure 1 The flowchart of the split Bregman algorithm for solving the proposed problem (11).The parameters k max and n max denote the outer and inner iteration numbers, respectively.Full-size DOI: 10.7717/peerjcs.1579/fig-1

Figure 6
Figure 6 Denoising results of different models on an image of the IXI dataset with the noise level σ = 5.The parts of the denoised images, which are framed by green boxes, are enlarged for visual comparison.The 1D curves of intensity value are shown below the denoised images.Image source credit: IXI dataset, CC BY-SA 3.0 (https://brain-development.org/ixi-dataset/).Full-size DOI: 10.7717/peerjcs.1579/fig-6

Phan (2023), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.1579 10/17 Table 1 Average PSNR and SSIM results of different methods on the IXI and SB datasets with different noise levels.
The values highlighted in bold represent the best results for each volume and noise level.