A Modified Higher-Order Singular Value Decomposition Framework With Adaptive Multilinear Tensor Rank Approximation for Three-Dimensional Magnetic Resonance Rician Noise Removal

The magnetic resonance (MR) images are acknowledged to be inevitably corrupted by Rician distributed noise, which adversely affected the image quality for diagnosis purpose. However, the traditional denoising methods may recover the images from corruptions with severe loss of detailed structure and edge information, which would affect the lesion detections and diagnostic decision making. In this study, we challenged improving the Rician noise removal from three-dimensional (3D) MR volumetric data through a modified higher-order singular value decomposition (MHOSVD) method. The proposed framework of MHOSVD involved a parameterized logarithmic nonconvex penalty function for low-rank tensor approximation (LRTA) algorithm optimization to suppress the image noise in MR dataset. Reference cubes were extracted from the noisy image volume, and block matching was performed according to nonlocal similarity for a fourth-order tensor construction. Then the LRTA problem was implemented by tensor factorization approaches, and the ranks of unfolding matrices along different modes of the tensor were estimated utilizing an adaptive nonconvex low-rank method. The denoised MR images were finally restored through aggregating all recovered cubes. We investigated the proposed algorithm MHOSVD on both the synthetic and real clinic 3D MR images for Rician noise removal, and relative results demonstrated that the MHOSVD can recover images with fine structures and detailed edge preservation with heavy noise even as high as 15% of the maximum intensity. The experimental results were also compared along with several classical denoising methods; the MHOSVD exhibited a sufficient improvement in noise-removal performance at various noise conditions in terms of different measurement indices such as peak signal-to-noise ratio and structural similarity index metrics. Based upon the comparison, the proposed MHOSVD has proved a relative state-of-the-art performance with excellent detailed structure reservation.

Introduction

Magnetic resonance imaging (MRI), as a widely accepted noninvasive imaging modality, was acknowledged as a useful diagnostic tool with high resolution and excellent contrast sensitivity to anatomical properties (1). However, despite the significant improvements in instrument technique during recent years, the detection or assessment of specific diagnostic information through referring physicians or computer-aided analysis still frequently suffered from serious random noise. As is known, the assessment of medical image quality is usually described by various physical measures including the contrast between different tissues, the detailed representation relative to imaging spatial resolution, and the image noise characteristics. The introduction of noise into MR images can usually cause imaging blur or fine structure coverage/distortion that may severely degrade the image quality and affect the diagnostic accuracy. As reported, a high level of noise may directly affect the signals resulting in anatomical inconsistencies that especially correspond to cardiac and brain images or may bias the orientations of tensors in functional MRI. Meanwhile, computer-aided diagnosis (CAD) approaches recently have been developed as assistive tools to help radiologists to detect and clarify abnormalities. Typical CAD systems based on medical images for lesion separation and recognition generally employed techniques including segmentation, feature extraction, and classification. However, one key problem that challenges the accuracy of CAD tasks was sensitivity to artifacts such as noise. Hence, to improve the performance of CAD approaches, image preprocessing was necessary wherein proper denoising algorithms were essential for the reliability and robustness of computer-aided detection/diagnosis. An ideal denoising method should restore the images by removing the stochastic corruptions of noise with preservation of the shapes and detailed structures of the tissue against abnormities as well. Previous studies have considered the complex raw data of MRI images be corrupted by white Gaussian noise with zero mean and same variance in both the real and imaginary parts; thus, the noise characters were transformed into Rician distribution in magnitude images (2, 3). Therefore, an effective denoising algorithm specific for Rician noise based upon three-dimensional (3D) MR image datasets can play a fundamental role to improve the diagnostic accuracy for both the radiologists and the CAD tools.

To date, several techniques have generally been implemented in commercial MRI devices to improve the signal-to-noise ratio (SNR) for image quality control (4). Along with the hardware evolutions such as magnetic field intensity that dramatically increased, efficient imaging processing algorithms were also utilized to recover the undesirable random variations that obscured the diagnostic information. One applicable way of denoising was to take the average from multiple repeated acquisitions, with the price of screening time extension, which was difficult for patients to keep static. The more popular way was to involve the reduction of noise power level during post-imaging processing. As MRI intrinsically posed a paradox between SNR and resolution, the ideal noise removal procedure should be able to perform SNR improvements without harming the fine structures and detailed features in images. Numerous methods for MR image denoising have been developed based upon different characters of noise. Filter-based approaches constructed the noise-reduction schemes with linear or nonlinear filters, such as anisotropic diffusion and total variation (TV) techniques, to improve image quality with edge preservation considerations (5–8). A few other studies utilized the statistics/estimation of noise properties to conduct denoising; examples included maximum likelihood (9, 10), linear minimum mean square error (11, 12), and phase error (13), which were generally used as estimators for Rician noise. Recently, the transform approaches were demonstrated to be powerful by exploiting the different representations of noise against signal either spatially or spectrally in transform domain (14–16). For example, wavelet transforms were used to decompose the signal and noise into multiresolution subspaces, which facilitated the Rician noise removal while preserving edge and fine details (17). Besides, curvelet transform and contourlet transforms were also used for denoising MR images (18, 19). Moreover, deep learning methods were also involved in the Rician noise removal in MR images. For example, Jiang et al. (20) developed a multichannel conventional neural network (CNN) method for noise removal in synthetic and clinical MR volumetric image. Manjón et al. proposed a two-step 3D MR image denoising algorithm by combining CNN with rotation-invariant nonlocal means (NLM) filtering. The relative results were encouraging, which improve the visual quality significantly; however, the deep leaning frameworks still need to meet the high burdens of huge training dataset and time-consuming computation (21).

More recent studies have intended to combine the strengths of multiple approaches (e.g., filter and transform approaches) to better regulate denoising in MRI (22, 23). The well-known NLM filter that Buades et al. (24) proposed was one classical method that considered exploiting nonlocal patch similarity for noise removal (24). The redundancy patterns within images ensured useful information to be restored by NLM with outstanding edge sharpness maintained. Kervrann et al. (25) soon extended the NLM algorithm into variations combined with transform approaches with the purpose of taking advantage of both similarity and sparsity. One famous NLM variation was block matching 3D (BM3D), which constructed its scheme by grouping image fragments according to similarity, filtering with the sparse representation in transform domain, and transforming back by aggregation procedure (26). Maggioni et al. (27) then modified the BM3D method using NLM and cosine/wavelet transform into BM4D implementation, and they demonstrated state-of-the-art in volumetric MRI data recovery. Meanwhile, the development in tensor factorization, such as the higher-order singular value decomposition (SVD) (HOSVD) (28), enabled better sparse representations using learned basis rather than the fixed orthogonal basis that cosine/wavelet transforms generally used in BM3D/BM4D. Rajwade et al. (29) derived the HOSVD denoising algorithm, which manipulated the NLM filtering with HOSVD transform instead. As the learning basis was adaptable to image contents, redundant procedures of patches rearrangements were saved, and further performance improvements were promising. In Zhang et al. (30) have extended the HOSVD algorithm for MRI, and excellent noise-reduction effectiveness was obtained. However, the HOSVD-based denoising method generally computed the transform coefficients using hard threshold function with nonconvex and nonlinear properties, which was considered to constrain the denoising performance somehow. Further, the hard threshold function would also cause additional shock and pseudo-Gibbs effect, which may limit the retention of image structure.

The present study aimed to improve the performance of HOSVD application for denoising 3D MR volumetric data. To address the drawback mentioned above, we modified the HOSVD algorithm [named as modified HOSVD (MHOSVD)] with an adaptive multilinear tensor rank approximation method utilizing nonlocal similarity and nonconvex logarithmic regularization. Image cubes from volumetric data were stacked based on similarity to construct a fourth-order tensor. Then tensor factorization was performed. A low-rank tensor was derived to approximate the sparse representations by exploiting the image structural redundancy. We estimated the ranks of unfolding matrices along different modes of the tensor using an adaptive nonconvex low-rank method. Experiments were performed based on various MR images, obtained from either computer synthesis or clinical screening, to investigate denoising improvements of proposed framework against several advanced algorithms with the following novelties:

1. Considered the nonlocal similarity in 3D MR datasets and formulated the grouped similar cubes into a low-rank tensor approximation (LRTA) problem.

2. Applied a nonconvex low-rank function to adaptively estimate the rank of unfolding matrices along different modes.

3. Used the proximal operator to obtain the global closed-form solution for the constructed low-rank approximation problem.

We structured the rest of this manuscript as follows: Background reviews necessary backgrounds about the tensor. Proposed Model describes our algorithm proposed for Rician noise removal based upon adaptive HOSVD framework. We validated the proposed schemes to both synthetic and clinical 3D MR images, and the effectiveness of noise reduction was evaluated and compared with that of existing denoising algorithms in Experiments and Results. Last but not the least, we draw our conclusion in Conclusions.

Background

A tensor can be considered as a generalization of vector and matrix for the representation of multidimensional quantities. Throughout the article, we represented the scalars as lowercase letters (e.g., x), vectors as bold lowercase letters (e.g., x), matrices as bold uppercase letters (e.g., X), and tensors as bold calligraphic letters (e.g., ${X}$). For example, the symbol ${X}\in R^{J_1\times J_2\times \dots \times J_N}$ signified an arbitrary Nth-order tensor over the multidimensional space of size J₁ × J₂ × … × J_N. We also used the subscript lowercase letters to present the element tensor; e.g., an arbitrary element of the tensor ${X}$ can be denoted as ${{X}}_{j_1{j}_2\dots j_N}$ using (j₁j₂…j_N).

Definition 1 (Tensor fiber). A tensor generally consists of fibers that can be thought as the high-dimensional generalization of matrix rows and columns. The mode-n fiber of tensor ${X}$ is defined by fixing the index in all dimensions but the n-th dimension.

Definition 2 (Tensor unfolding). Tensor unfolding, also called tensor matricization, is to rearrange the elements of the tensor into the matrix format under predefined order. The mode-n matricization of an Nth-order tensor ${X}\in R^{J_1\times J_2\times \dots \times J_N}$ would be represented as ${\text{X}}_{(n)}\in R^{J_n\times (J_1\times \dots \times J_{n-1}\times J_{n+1}\times \dots \times J_N)}$, which ordered the mode-n fibers of tensor ${X}$ into its columns.

Definition 3 (n-mode matrix product). The product of a tensor ${X}\in R^{J_1\times J_2\times \dots \times J_N}$ and a matrix $\text{V}\in R^{K\times J_n}$ along the n-th dimension is a tensor of size J₁ × … × J_n−1 × K × J_n+1 × … × J_N, denoted as ${X}{\times }_n\text{V}$, with its element expressed by following equation:

$\begin{array}{ll}{({X}{\times }_n\text{V})}_{j_1\dots j_{n-1}k{j}_{n+1}\dots j_N}={\displaystyle \sum _{j_n=1}^{J_n}}{x}_{j_1{j}_2\dots j_N}{v}_{k{j}_n}.& (1)\end{array}$

Definition 4 (Multilinear tensor rank). The multilinear tensor rank of ${X}\in R^{J_1\times J_2\times \dots \times J_N}$ is defined as the rank of the mode-n unfolding matrix X_(n), denoted as $ran{k}_n({X})$.

Proposed Model

Low-Rank Matrix Approximation

Low-rank matrix approximation (LRMA) aimed to solve the minimization problem through optimizing the cost function between the given data and an approximation with reduced rank. The method was helpful in various image processing such as 2D image denoising (31, 32). The LRMA usually applied for image noise removal following steps, as follows: searching and grouping nonlocal similar patches, performing collaborative filtering, and aggregating the recovered patches to restore the denoised images. Mathematically, the LRMA method can be modeled as follows:

$\begin{array}{ll}{}_X^{\min }\frac{1}{2}\Vert Y-X{\Vert }_F^2+\lambda rank(X),& (2)\end{array}$

where rank(·) represented the rank of matrix, referring to the number of nonzero singular values of matrix. λ was the trade-off parameter balancing the contribution between fidelity and regularization term. However, due to the high nonconvex and nonlinear properties, problem (Eq. 2) was practically intractable. Therefore, convex relaxation nuclear norm was used instead, and problem (Eq. 2) can be transformed into the following LRMA problem:

$\begin{array}{ll}{}_X^{\min }\frac{1}{2}{\Vert Y-X\Vert }_F^2+\lambda {\Vert X\Vert }_{\ast },& (3)\end{array}$

where ${\left|\right|\text{X}\left|\right|}_{*}=\sum _i{\delta }_i(\text{X})$ denoted the nuclear norm that was defined as the sum of a singular value from matrix X. Cai et al. (33) have proved that the solution of the LRMA problem (Eq. 3) can be conveniently obtained by singular value thresholding; i.e., $\hat{X}=U{S}_{\lambda }(\Sigma )V^{\mathrm{\text{T}}}$, where Y = UΣV^T was the SVD of the noisy image Y, and S_λ(Σ) = max(Σ−λ, 0).

Low-Rank Tensor Approximation

LRTA can extend the LRMA algorithm for multidimensional image processing. In this study, we employed the LRTA for 3D MR image denoising. First, reference cube sized of p × p × p was extracted from the MR noisy image. Block matching was performed with a local search window of size k × k × k, and m similar cubes (including the reference cube) were found and stacked to a fourth-order tensor ${Y}\in R^{p\times p\times p\times m}$. Then the problem estimating the noise-free version ${X}$ according to noisy image ${Y}$ can be regarded as an LRTA mathematically modeled as

$\begin{array}{ll}{{X}}^{*}{=}_{{\text{X}}_{}}^{\text{argmin}}\frac{1}{2}{\Vert {Y}-{X}{\Vert }_F^2\text{ s}.\text{t}.\text{ rank}}_i({X})\le r_i(i=1,2,3,4),& (4)\end{array}$

where $ran{k}_i({X})$ denoted the rank of the mode-i unfolding matrix of the fourth-order tensor ${X}$. Based upon the HOSVD (28), tensor ${X}$ can be decomposed into the following:

$\begin{array}{ll}{X}={S}{\times }_1{\text{U}}^{(1)}{\times }_2{\text{U}}^{(2)}{\times }_3{\text{U}}^{(3)}{\times }_4{\text{U}}^{(4)},& (5)\end{array}$

where ${S}\in R^{r_1\times r_2\times r_3\times r_4}$ was the core tensor and ${\text{U}}^{(i)}\in R^{p\times r_i}(i=1,2,3,4)$ was the factorization matrix. The key for solving problem (Eq. 5) was to estimate the multilinear tensor rank parameter r_i(i = 1, 2, 3, 4). Therefore, the LRTA problem (Eq. 4) can be independently decomposed into four relative minimization problems:

$\begin{array}{ll}{X}_{(n)}^{\ast }{=}_{X_{\left(n\right)}}^{argmin}\frac{1}{2}|\left|Y_{\left(n\right)}-X_{\left(n\right)}\right|{|}_F^2+\lambda rank(X_{(n)})\left(n=1,2,3,4\right),& (6)\end{array}$

where rank(X_(n)) was the number of nonzero singular values of the mode-n unfolding matrix X_(n), and ${\left|\right|{Y}-{X}\left|\right|}_F^2={\left|\right|{\text{Y}}_{(n)}-{\text{X}}_{(n)}\left|\right|}_F^2$ were the Frobenius norm.

Adaptive Multilinear Tensor Rank Algorithm

As the solution of problem (Eq. 6) was intractable, the nuclear norm with nonconvex properties have been demonstrated an effective tool for sparsity of singular values compared with the convex rank function. Motivated by Selesnick and Bayram (34), Chen and Selesnick (35), and Parekh and Selesnick (36), we propose a parameterized logarithmic nonconvex penalty function on singular values, formulated as follows:

$\begin{array}{ll}{X}_{(n)}^{\ast }{=}_{X_{\left(n\right)}}^{argmin}\frac{1}{2}|\left|Y_{\left(n\right)}-X_{\left(n\right)}\right|{|}_F^2+\lambda {\displaystyle \sum _{i=1}^{r_n}\frac{1}{a}\log \left(1+a{\delta }_i\left(X_{\left(n\right)}\right)\right),}& (7)\end{array}$

where δ_i(X_(n)) was the i-th singular value of the mode-n unfolding matrix X_(n). Although the nonconvex penalty function was employed, the proposed LRTA problem (Eq. 7) was strictly convex when 0 ≤ a < 1/λ, which could be proven by Lemma 1, as follows.

Lemma 1(37): The nonconvex function $\psi (x,a)=\frac{1}{a}log(1+ax)$ satisfied Assumption 1 as mentioned in literature (36). The function h: R→R represented as h(x, a) = ψ(x, a)−|x| was continuously differentiable and concave and satisfied −a ≤ h″(x, a) ≤ 0.

Eq. (7), which was considered to be strictly convex when 0 ≤ a < 1/λ, was proved, as follows.

Proof: Define the function J:R^{m × n}→R as

$\begin{array}{l}\begin{array}{c}J(X)=\frac{1}{2}\Vert Y-X{\Vert }_F^2+\lambda {\displaystyle \sum _{i=1}^{r_n}}(h({\delta }_i(X),a)+|{\delta }_i(X)|)\\ =\frac{1}{2}tr(Y^{T}Y)-tr(X{Y}^T)+\frac{1}{2}tr(X^{T}X)\hfill \\ +\lambda {\displaystyle \sum _{i=1}^{r_n}}(h({\delta }_i(X),a)+|{\delta }_i(X)|)\hfill \\ =\frac{1}{2}tr(Y^{T}Y)-tr(X{Y}^T)+\frac{1}{2}tr(X^{T}X)\hfill \\ +\lambda {\displaystyle \sum _{i=1}^{r_n}}h({\delta }_i(X),a)+\lambda \Vert X{\Vert }_{*}\hfill \\ =J_1(X)+\frac{1}{2}tr(Y^{T}Y)-tr(X{Y}^T)+\lambda \Vert X{\Vert }_{*},\hfill \end{array}\end{array}$

where $J_1(X)=\frac{1}{2}tr(X^{T}X)+\lambda {\sum }_{i=1}^{r_n}h({\delta }_i(X),a)$. Since both the linear function tr(XY^T) and trace norm ||X||_* were convex function and the function tr(Y^TY) was independent on X, the J(X) was strictly convex when the function J₁(X) was strictly convex. According to Parekh and Selesnick (36), we can conclude that the function J₁(X) was strictly convex when 0 ≤ a < 1/λ.

According to Selesnick and Bayram (34), the proximal operator of the proposed penalty function (Eq. 7) was defined as:

$\begin{array}{ll}\theta (y;\lambda ,a){=}_{x\in R}^{argmin}\frac{1}{2}{(y-x)}^2+\lambda \frac{1}{a}log(1+ax),& (8)\end{array}$

and the equivalent threshold function of the proximal operator (8) can be defined as (34, 38)

$\begin{array}{ll}\theta \left(y;\lambda ,a\right)=\{\begin{array}{c}\frac{\left|y\right|}{2}-\frac{1}{2a}+\sqrt{{\left(\frac{\left|y\right|}{2}+\frac{1}{2a}\right)}^2-\frac{\lambda }{a}}\left|y\right|\ge \lambda .\\ 0\left|y\right|<\lambda \end{array}& (9)\end{array}$

When the proximal operator θ(•) was applied to matrix X, the element-wise would need to be treated. Then the proposed LRTA problem (Eq. 7) can globally solved with optimal solution:

$\begin{array}{ll}{\text{X}}_{(n)}^{*}={\text{U}}^{(n)}\theta ({\Sigma }^{(n)};\lambda ,a){({\text{V}}^{(n)})}^T,& (10)\end{array}$

where ${\text{Y}}_{(n)}={\text{U}}^{(n)}{\Sigma }^{(n)}{({\text{V}}^{(n)})}^T$. Therefore, the multilinear tensor rank parameter r_n can be adaptively estimated as the number of the nonzero elements in θ(Σ⁽ⁿ⁾; λ, a), and the estimated factorization matrix ${\hat{\text{U}}}^{(\mathrm{\text{n}})}$was selected as the first r_n columns of the U⁽ⁿ⁾. Thus, the estimated core tensor can be formulated as:

$\begin{array}{ll}\hat{{S}}={Y}\times 1{({\hat{U}}^{(1)})}^{\text{T}}\times 2{({\hat{U}}^{(2)})}^{\text{T}}\times 3{({\hat{U}}^{(3)})}^{\text{T}}\times 4{({\hat{U}}^{(4)})}^{\text{T}},& (11)\end{array}$

and the restored tensor $\hat{{X}}$ can be obtained by:

$\begin{array}{ll}\hat{{X}}=\hat{{S}}\times 1{({\hat{U}}^{(1)})}^{\text{T}}\times 2{({\hat{U}}^{(2)})}^{\text{T}}\times 3{({\hat{U}}^{(3)})}^{\text{T}}\times 4({\hat{U}}^{(4)})\text{T}.& (12)\end{array}$

Variance-Stabilizing Transform

MR images were mainly contaminated by Rician-distributed noise (3). For Rician noise removal in 3D MR image, we implemented the forward and inverse variance-stabilizing transform (VST) (38) into our proposed denoising framework. The forward VST was used to convert the 3D MR images with signal-dependent Rician-distribution noise into the new volumetric data with homoscedastic Gaussian-distributed noise. Once the proposed denoising algorithm was applied, the inverse VST (VST⁻¹) counteracted the recovered data to obtain the denoised images. Therefore, the whole scheme of the proposed denoising algorithm in 3D MR images can be formulated as:

$\begin{array}{ll}\hat{{X}}=\text{VS}{\text{T}}^{-1}(\text{Denoising}(\text{VST}({Y}\text{,}{\sigma }_{\text{n}})\text{,}{\sigma }_{\text{VST}})\text{,}{\sigma }_{\text{n}}),& (13)\end{array}$

where ${Y}$ is the 3D MR images contaminated by Rician-distributed noise, σ_n represents the Rician noise level, and σ_VST is the standard deviation of noise after VST. We describe the proposed denoising method in Algorithm 1.

Experiments and Results

In this section, we validated the performance of the proposed algorithm (refer as MHOSVD) on both the synthetic and clinical 3D MR images. The noise-free MR images for noise synthesis, including T1-weighted (T1w) data, T2-weighted (T2w) data, and proton density-weighted (PDw) data, were downloaded from BrainWeb database (39, 40). The raw data were formatted in a size of 181 × 217 × 181 with 1 mm³ × 1 mm³ × 1 mm³ voxel resolution. The noise data were simulated by adding different levels of spatially invariant Rician noise (from 1% to 15% of the maximum intensity with an increase of 2%).

The peak SNR (PSNR) and structural similarity index measurement (SSIM) were calculated to quantitatively evaluate the denoising performance. PSNR was a metric measuring the ratio of the maximum possible signal power to the noise power, which is defined as:

$\begin{array}{ll}PSNR=10{log}_{10}\frac{MA{X}^2}{MSE},& (14)\end{array}$

where MSE is the mean squared error (MSE) between the denoised image and ground truth and MAX represents the maximum image intensity. According to the definition, a higher PSNR value would indicate lower noise content that related to a better image quality.

Meanwhile, the metric of SSIM was defined according to human eye perception:

$\begin{array}{ll}SSIM=\frac{(2{\mu }_x{\mu }_{x^{*}}+c_1)(2{\sigma }_x{\sigma }_{x^{*}}+c_2)}{({\mu }_x^2+{\mu }_{x^{*}}^2+c_1)({\sigma }_x^2+{\sigma }_{x^{*}}^2+c_2)},& (15)\end{array}$

where μ_x and σ_x denote the mean and standard deviation of the ground truth image, respectively, while the ${\mu }_{x^{*}}$ and ${\sigma }_{x^{*}}$ denote the mean and standard deviation of the denoised image, respectively. c₁ and c₂ are constants. Therefore, the range of SSIM was from −1 to 1. An SSIM value close to 1 referred to a perfect similarity comparing the denoised image with the ground truth, which resulted in a higher performance of image restoration.

Parameter Sets

The parameters that may determine the denoising performance of the proposed algorithm (MHOSVD) included the following: the size of cube p, the number of similar cubes in a group m, the size of searching volume L, and noise-feedback parameters β and γ. A large value of p and m tended to benefit the removal of image data corruption at high-level noise. However, a large value of parameter m that denoted the number of similar cubes would also result in intragroup dissimilar cubes as a trade-off. Efforts were paid to balance the advantages and disadvantages; the setting of parameters p and m at different noise levels was chosen and shown in Table 1. Similarly, the increase in searching volume parameter L can help better in cube grouping according to image similarity for higher PSNR and also with a cost of the computational burden. The size of the search window (L) was set as 13 according to previous studies (27). The noise-feedback parameters β and γ, which were critical to the feedback of the residual image and residual noise, were set to 0.65 and 0.2 on the basis of previous research (30), respectively. In addition, the parameter a from parameterized logarithmic nonconvex penalty function (Eq. 7) and threshold τ from the proximal operator (9) were set as 0.5/τ and 2 × log(p³m), respectively.

TABLE 1

Table 1. The setting of the cube size p and numbers of similar cubes m at varying noise levels.

Denoising on Synthetic Three-Dimensional Magnetic Resonance Images

The performance of proposed algorithm (MHOSVD) was first tested on synthetic brain 3D MR images, against several existing classical algorithms including RSNLMMSE (11), BM4D (27), PRI-NLM3D (22), and HOSVD-R (30). The corresponding PSNR performance and SSIM performance on T1w, T2w, and PDw images were compared across algorithms. As illustrated by Figure 1, the proposed MHOSVD exhibited encouraging improvements over the other three methods, which had been previously verified to be state-of-the-art. Further, we also tabulated the PSNR (Table 2) and SSIM (Table 3) obtained by these four denoising algorithms (RSNLMMSE, BM4D, PRI-NLM3D, HOSVD-R, and MHOSVD) on MR images (T1w, T2w, and PDw) adding Rician noise under different noise levels (1, 3, 5, 7, 9, 11, 13, and 15%, respectively). The proposed MHOSVD outperformed RSNLMMSE (ranged from 1.9337 to 3.8249 dB), PRI-NLM3D (1.091 to 2.2149 dB), BM4D (0.4305 to 1.3077 dB), and HOSVD-R (0.0737 to 0.46 dB) with regard to PSNR. The same superiority was also observed according to the metric SSIM, indicating that the proposed MHOSVD can better remove Rician noise with structure preservation against corruptions across a board range of noise level.

FIGURE 1

Figure 1. Peak signal-to-noise ratio (PSNR) and structural similarity index measurement (SSIM) comparisons across four denoising methods (RSNLMMSE, BM4D, PRI-NLM3D, HOSVD-R, and MHOSVD) for Rician noise removal in three-dimensional magnetic resonance (3D MR) images synthesized with corruptions under different noise levels.

TABLE 2

Table 2. Peak signal-to-noise ratio (PSNR) comparisons across RSNLMMSE, BM4D, PRI-NLM3D, HOSVD-R, and MHOSVD for Rician noise removal in synthetic three-dimensional magnetic resonance (3D MR) images.

TABLE 3

Table 3. Structural similarity index measurement (SSIM) comparisons across BM4D, PRI-NLM3D, HOSVD-R, and MHOSVD for Rician noise removal in synthetic three-dimensional magnetic resonance (3D MR) images.

Examples of denoising results of all four algorithms on synthesized MR images (T1w, T2w, and PDw) corrupted by the Rician noise at level of 15% are visually shown in Figures 2–4. Through visible observation, the RSNLMMSE and the PRI-NLM3D exhibited a relatively degraded performance as undesirable oversmoothing and higher detailed structure loss were found to be accompanied with noise reduction, than did the other three methods (BM4D, HOSVD-R, and MHOSVD). Though all four algorithms showed intensity oscillations in homogenous regions according to the results of residual images (Figures 2–4) and enlarged regions of interest (ROIs) (Figure 5), the proposed MHOSVD decreased the unpleasant fluctuation in performance somehow. Nevertheless, the figures illustrate that our algorithm MHOSVD retained fine textures pretty well in 3D volumetric images.

FIGURE 2

Figure 2. Denoising T1-weighted (T1w) images by different methods (RSNLMMSE, BM4D, PRI-NLM3D, HOSVD-R, and MHOSVD) under 15% Rician noise. First line: T2w image without noise and with 15% Rician noise. Second line: the denoised images with four different methods. Third line: the corresponding residual images (the absolute difference between the denoised and noise-free images).

FIGURE 3

Figure 3. Denoising T2-weighted (T2w) images by different methods under 15% Rician noise. First line: T2w image without noise and with 15% Rician noise. Second line: the denoised images with four different methods. Third line: the corresponding error images.

FIGURE 4

Figure 4. Denoising proton density-weighted (PDw) images by different methods under 15% Rician noise. First line: T2w image without noise and with 15% Rician noise. Second line: the denoised images with four different methods. Third line: the corresponding error images.

FIGURE 5

Figure 5. Performance comparison with enlarged region of interest (ROI) visualization using four different methods under 15% Rician noise. ROIs selected from above (red rectangles in Figures 2–4). First line: T1-weighted (T1w) images. Second line: T2-weighted (T2w) images. Third line: proton density-weighted (PDw) images.

Compared with hard threshold function, the proposed logarithmic penalty function was closer to the rank function. The singular values on unfolding matrices along different modes also were compared across ground truth and restored images (noisy level set as 15%) via HOSVD-R and MHOSVD method, respectively. As shown in Figure 6, the singular values extracted from recovered images using our proposed method were obviously closer to those of the noise-free images than HOSVD-R, especially in the small singular value domain.

FIGURE 6

Figure 6. The singular value decomposition (SVD) comparison of different mode unfolding matrices about ground truth and denoised image, destroyed by 15% Rician noise, through the HOSVD-R and the proposed MHOSVD method. (A) Mode 1 matrix; (B) mode 2 matrix; (C) mode 3 matrix.

In addition, we further applied the proposed algorithm to verify its reliability and effectiveness on the basis of MR images with various abnormities. The ground truth was obtained from the Multimodal Brain Tumor Image Segmentation Challenge (BraTS) 2013, which provided MRI datasets that have already been optimized by several image preprocessing approaches for further application simulations such as CAD segmentations (41). Figure 7 shows some examples randomly selected from the BraTS, as Figure 7A focuses on T1w contrast-enhanced modes from three different patients, while Figure 7B involves a specific case imaged via different modes including T1w contrast enhanced, T2w, and fluid-attenuated inversion recovery (FLAIR). In this study, Rician noise with a level of 11% was introduced into the BraTS data to synthesize the noisy images. As illustrated by the second column of Figure 7, the boundary of lesions can be severely blurred with the impact of noise, and some small structures were even buried by the stochastic variances that noise introduced. The proposed algorithm (MHOSVD) was applied to remove the image corruption and restore the image qualities, and the relative results were compared with those of several classical denoising method including RSNLMMSE, PRI-NLM3D, BM4D, and HOSVD-R. Figure 7 illustrates that all algorithms can repair images with different noise removal abilities. Comparing columns 3–7 with column 2, it can be visibly observed that the blurry images were degraded after the denoising procedure (Figure 7), the boundary of pathological context was re-highlighted, and fine structures covered by the additive noise were restored as a result of noise removal. Further, compared with other classical methods (RSNLMMSE, PRI-NLM3D, BM4D, and HOSVD-R), our proposed algorithm can emphasize the lesion edge more effectively and preserve more small structures, which may assist the physician in decision making under complex conditions when images were acquired with heavy noise and play a key role for accuracy improvements in computer-aided detection/diagnostic tools such as tumor segmentation and recognition.

FIGURE 7

Figure 7. Denoising performance based upon synthesized brain images with abnormities acquired via different MRI modes. The ground truth images (1st column) were downloaded from the open source of BraTS, which included clinical datasets after several image preprocessing approaches. A level of 11% Rician noise was added into the ground truth to synthesize the noised images (second column). The proposed algorithm (MHOSVD) was applied to remove the corruptions against four classical denoising methods (RSNLMMSE, PRI-NLM3D, BM4D, and HOSVD-R), and the results were compared. (A) Denoising performances based on T1-weighted contrast-enhanced images acquired from three different patients, and each row corresponds to one patient. (B) Denoising performances based on MR datasets acquired by different imaging modes [from top to bottom: T1-weighted contrast-enhanced, T2, and fluid-attenuated inversion recovery (FLAIR)] of the same patient.

Denoising on Clinical Three-Dimensional Magnetic Resonance Images

We also validated the denoising performance of the proposed algorithm on clinical MR datasets downloaded from the publicly available Open Access Series of Imaging Studies (OASIS) database (http://www.oasis-brains.org) (42). The clinical datasets (brain images, T1w) were of size 256 × 256 × 128 with a resolution of 1.0 mm³ × 1.0 mm³ × 1.25 mm³. The obtained MR images were originally corrupted by a certain level of clinical noise through acquisition. The Rician noise level was investigated according to (38), and the noise estimates for these MR images (OAS1_0092 and OAS1_0112) were ~4.5 and 3% of the maximum intensity, respectively.

The restored results of our proposed MHOSVD for T1w brain images in sagittal, coronal, and transverse planes are presented in Figure 8. Further comparisons across four different algorithms (RSNLMMSE, BM4D, PRI-NLM3D, HOSVD-R, and MHOSVD) are shown in Figure 9. According to performance comparison (Figure 9), the statistical approach RSNLMMSE can restore the images but also severely erased structural details, and the PRI-NLM3D also showed some blurry edge and several tiny structures loss after noise removal. Again, based upon clinical images, the outstanding noise-reduction performance with excellent detailed structure reservation was practically observed using MHOSVD.

FIGURE 8

Figure 8. Example of clinical image denoising using the proposed method. First line: clinic noise image. Second line: the denoised images with the proposed method. Third line: the corresponding error images (the absolute difference between the denoised and noise images). Regions of interest (ROIs) (red rectangle) were selected and enlarged for better detail representation in this figure.

FIGURE 9

Figure 9. Denoising performances comparison based on regions of interest (ROIs) of top images (red rectangle in Figure 8). (A) Enlarged ROI from OAS1_0092. (B) Enlarged ROI from OAS1_0112. The proposed algorithm (MHOSVD) was compared with existing state-of-the-art methods including RSNLMMSE, BM4D, PRI-NLM3D, and HOSVD-R.

Conclusions

In this study, we proposed an adaptive multilinear tensor rank approximation algorithm based on parameterized nonconvex logarithmic function for Rician noise removal in 3D MR images. The framework extracted 3D cube from noise images and searched similar cubes according to the Euclidean distance between cubes to construct the corresponding fourth-order tensor. The nonconvex logarithmic function was applied for the rank estimation of unfolding matrices along different modes of the tensor. Then the fourth-order tensors can be effectively recovered through the adaptive multilinear tensor rank approximation. Afterwards, the recovered cubes were aggregated to obtain the recovered volumetric images. Finally, we verified our algorithm on synthetic and clinical MR images. The proposed method exhibited a state-of-the-art performance in Rician noise removal with excellent fine detail preservation, which even outperformed several existing classical denoising methods (RSNLMMSE, BM4D, PRI-NLM3D, and HOSVD-R). Further, the capability to effectively capture the sparsity of multidimensional dataset would enable our work to be extended to several other domains, such as hyperspectral images (43, 44), color images (45), and video (46).

Data Availability Statement

The image dataset testing the algorithm effectiveness and efficiency was originally accessed from Open Access Series of Imaging Studies (OASIS) database (http://www.oasis-brains.org).

Author Contributions

LW designed the experiment, collected and analyzed the data, and drafted the manuscript. DX and XW participated in experiment setup and interpreted the results partly. LC and WH supervised the project and revised the manuscript. All listed authors approved the final manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (31800824), the Chongqing Science and Technology Program (cstc2018jcyjAX0390), and the project of Advanced Scientific Research Cooperation and high-level personnel training between China and America-Oceania Area, funded by Ministry of Education of the People's Republic of China.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling editor declared a shared affiliation with the authors at the time of review.

References

1. Brosnan T, Wright G, Nishimura D, Cao Q, Macovski A, Sommer FG. Noise reduction in magnetic resonance imaging. Magn Reson Med. (1988) 8:394–409. doi: 10.1002/mrm.1910080404

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Henkelman RM. Measurement of signal intensities in the presence of noise in MR images. Med Phys. (1985) 12:232–3. doi: 10.1118/1.595711

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Gudbjartsson H, Patz S. The Rician distribution of noisy MRI data. Magn Reson Med. (1995) 34:910–4. doi: 10.1002/mrm.1910340618

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Mohan J, Krishnaveni V, Guo Y. A survey on the magnetic resonance image denoising methods. Biomed Signal Process Control. (2014) 9:56–69. doi: 10.1016/j.bspc.2013.10.007

CrossRef Full Text | Google Scholar

5. Tang J, Sun Q, Liu J, Cao Y editors. An adaptive anisotropic diffusion filter for noise reduction in MR images. In: 2007 International Conference on Mechatronics Automation. Harbin: IEEE (2007).

Google Scholar

6. Yuan J. An improved variational model for denoising magnetic resonance images. Comput Math Appl. (2018) 76:2212–22. doi: 10.1016/j.camwa.2018.05.044

CrossRef Full Text | Google Scholar

7. Liu RW, Shi L, Huang W, Xu J, Yu SCH, Wang D. Generalized total variation-based MRI Rician denoising model with spatially adaptive regularization parameters. Magn Reson Imaging. (2014) 32:702–20. doi: 10.1016/j.mri.2014.03.004

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Rudin LI, Osher S, Fatemi E. Nonlinear total variation based noise removal algorithms. Phys D Nonlinear Phenom. (1992) 60:259–68. doi: 10.1016/0167-2789(92)90242-F

CrossRef Full Text | Google Scholar

9. Sijbers J, den Dekker AJ, Van Audekerke J, Verhoye M, Van Dyck D. Estimation of the noise in magnitude MR images. Magn Reson imaging. (1998) 16:87–90. doi: 10.1016/S0730-725X(97)00199-9

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Sijbers J, den Dekker AJ, Scheunders P, Van Dyck D. Maximum-likelihood estimation of Rician distribution parameters. IEEE Trans Med Imaging. (1998) 17:357–61. doi: 10.1109/42.712125

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Aja-Fernández S, Alberola-López C, Westin C-F. Noise and signal estimation in magnitude MRI and Rician distributed images: a LMMSE approach. IEEE Trans Image Process. (2008) 17:1383–98. doi: 10.1109/TIP.2008.925382

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Golshan HM, Hasanzadeh RP, Yousefzadeh SC. An MRI denoising method using image data redundancy and local SNR estimation. Magn Reson Imaging. (2013) 31:1206–17. doi: 10.1016/j.mri.2013.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Tisdall D, Atkins MS editors. MRI denoising via phase error estimation. In: Medical Imaging 2005: Image Processing. San Diego, CA: International Society for Optics and Photonics (2005).

Google Scholar

14. Bao P, Zhang L. Noise reduction for magnetic resonance images via adaptive multiscale products thresholding. IEEE Trans Med Imaging. (2003) 22:1089–99. doi: 10.1109/TMI.2003.816958

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Anand CS, Sahambi J editors. MRI denoising using bilateral filter in redundant wavelet domain. In: TENCON 2008-2008 IEEE Region 10 Conference. Hyderabad: IEEE (2008).

Google Scholar

16. Wood JC, Johnson KM. Wavelet packet denoising of magnetic resonance images: importance of Rician noise at low SNR. Magn Reson Med. (1999) 41:631–5. doi: 10.1002/(SICI)1522-2594(199903)41:3<631::AID-MRM29>3.0.CO;2-Q

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Nowak RD. Wavelet-based Rician noise removal for magnetic resonance imaging. IEEE Trans Image Process. (1999) 8:1408–19. doi: 10.1109/83.791966

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Ashamol V, Sreelekha G, Sathidevi P editors. Diffusion-based image denoising combining curvelet wavelet. In: 2008 15th International Conference on Systems, Signals Image Processing. Bratislava: IEEE (2008).

Google Scholar

19. Do MN, Vetterli M. The contourlet transform: an efficient directional multiresolution image representation. IEEE Trans Image Process. (2005) 14:2091–106. doi: 10.1109/TIP.2005.859376

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Jiang D, Dou W, Vosters L, Xu X, Sun Y, Tan T. Denoising of 3D magnetic resonance images with multi-channel residual learning of convolutional neural network. Jpn J Radiol. (2018) 36:566–74. doi: 10.1007/s11604-018-0758-8

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Manjón JV, Coupé P editors. MRI denoising using deep learning. In: International Workshop on Patch-based Techniques in Medical Imaging. Granada: Springer (2018).

Google Scholar

22. Manjón JV, Coupé P, Buades A, Collins DL, Robles M. New methods for MRI denoising based on sparseness and self-similarity. Med Image Anal. (2012) 16:18–27. doi: 10.1016/j.media.2011.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Coupé P, Yger P, Prima S, Hellier P, Kervrann C, Barillot C. An optimized blockwise nonlocal means denoising filter for 3-D magnetic resonance images. IEEE Trans Med Imaging. (2008) 27:425–41. doi: 10.1109/TMI.2007.906087

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Buades A, Coll B, Morel J-M. A review of image denoising algorithms, with a new one. Multiscale Model Simul. (2005) 4:490–530. doi: 10.1137/040616024

CrossRef Full Text | Google Scholar

25. Kervrann C, Boulanger J, Coupé P editors. Bayesian non-local means filter, image redundancy and adaptive dictionaries for noise removal. In: International Conference on Scale Space and Variational Methods in Computer Vision. Ischia: Springer (2007).

Google Scholar

26. Dabov K, Foi A, Katkovnik V, Egiazarian K. Image denoising by sparse 3-D transform-domain collaborative filtering. IEEE Trans Image Process. (2007) 16:2080–95. doi: 10.1109/TIP.2007.901238

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Maggioni M, Katkovnik V, Egiazarian K, Foi A. Nonlocal transform-domain filter for volumetric data denoising and reconstruction. IEEE Trans Image Process. (2012) 22:119–33. doi: 10.1109/TIP.2012.2210725

PubMed Abstract | CrossRef Full Text | Google Scholar

28. De Lathauwer L, De Moor B, Vandewalle J. A multilinear singular value decomposition. SIAM J Matrix Anal Appl. (2000) 21:1253–78. doi: 10.1137/S0895479896305696

CrossRef Full Text | Google Scholar

29. Rajwade A, Rangarajan A, Banerjee A. Image denoising using the higher order singular value decomposition. IEEE Trans Pattern Anal Mach Intell. (2012) 35:849–62. doi: 10.1109/TPAMI.2012.140

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Zhang X, Xu Z, Jia N, Yang W, Feng Q, Chen W, et al. Denoising of 3D magnetic resonance images by using higher-order singular value decomposition. Med Image Anal. (2015) 19:75–86. doi: 10.1016/j.media.2014.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Zha Z, Yuan X, Li B, Zhang X, Liu X, Tang L, et al. Analyzing the weighted nuclear norm minimization and nuclear norm minimization based on group sparse representation. arXiv:1702.04463. (2017).

Google Scholar

32. Xie Y, Gu S, Liu Y, Zuo W, Zhang W, Zhang L. Weighted Schatten $ p $-norm minimization for image denoising and background subtraction. IEEE Trans Image Process. (2016) 25:4842–57. doi: 10.1109/TIP.2016.2599290

CrossRef Full Text | Google Scholar

33. Cai J-F, Candès EJ, Shen Z. A singular value thresholding algorithm for matrix completion. SIAM J Optim. (2010) 20:1956–82. doi: 10.1137/080738970

CrossRef Full Text | Google Scholar

34. Selesnick IW, Bayram I. Sparse signal estimation by maximally sparse convex optimization. IEEE Trans Signal Process. (2014) 62:1078–92. doi: 10.1109/TSP.2014.2298839

CrossRef Full Text | Google Scholar

35. Chen P-Y, Selesnick IW. Group-sparse signal denoising: non-convex regularization, convex optimization. IEEE Trans Signal Process. (2014) 62:3464–78. doi: 10.1109/TSP.2014.2329274

CrossRef Full Text | Google Scholar

36. Parekh A, Selesnick IW. Enhanced low-rank matrix approximation. IEEE Signal Process Lett. (2016) 23:493–7. doi: 10.1109/LSP.2016.2535227

CrossRef Full Text | Google Scholar

37. Parekh A, Selesnick IW. Convex denoising using non-convex tight frame regularization. IEEE Signal Process Lett. (2015) 22:1786–90. doi: 10.1109/LSP.2015.2432095

CrossRef Full Text | Google Scholar

38. Foi A editor Noise estimation removal in MR imaging: the variance-stabilization approach. In: 2011 IEEE International Symposium on Biomedical Imaging: from Nano to Macro. Chicago, IL: IEEE (2011).

Google Scholar

39. Collins DL, Zijdenbos AP, Kollokian V, Sled JG, Kabani NJ, Holmes CJ, et al. Design and construction of a realistic digital brain phantom. IEEE Trans Med Imaging. (1998) 17:463–8. doi: 10.1109/42.712135

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Kwan R-S, Evans AC, Pike GB. MRI simulation-based evaluation of image-processing and classification methods. IEEE Trans Med Imaging. (1999) 18:1085–97. doi: 10.1109/42.816072

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Menze BH, Jakab A, Bauer S, Kalpathy-Cramer J, Farahani K, Kirby J, et al. The multimodal brain tumor image segmentation benchmark (BRATS). IEEE Trans Med Imaging. (2014) 34:1993–2024. doi: 10.1109/TMI.2014.2377694

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Marcus DS, Wang TH, Parker J, Csernansky JG, Morris JC, Buckner RL. Open Access Series of Imaging Studies (OASIS): cross-sectional MRI data in young, middle aged, nondemented, and demented older adults. J Cogn Neurosci. (2007) 19:1498–507. doi: 10.1162/jocn.2007.19.9.1498

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Luo F, Zhang L, Du B, Zhang L. Dimensionality reduction with enhanced hybrid-graph discriminant learning for hyperspectral image classification. IEEE Trans Geosci Remote Sens. (2020) 58:5336–53. doi: 10.1109/TGRS.2020.2963848

CrossRef Full Text | Google Scholar

44. Luo F, Zhang L, Zhou X, Guo T, Cheng Y, Yin T. Sparse-adaptive hypergraph discriminant analysis for hyperspectral image classification. IEEE Geosci Remote Sens Lett. (2019) 17:1082–6. doi: 10.1109/LGRS.2019.2936652

CrossRef Full Text | Google Scholar

45. Kong Z, Yang X. Color image and multispectral image denoising using block diagonal representation. IEEE Trans Image Process. (2019) 28:4247–59. doi: 10.1109/TIP.2019.2907478

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Bengua JA, Phien HN, Tuan HD, Do MN. Efficient tensor completion for color image and video recovery: low-rank tensor train. IEEE Trans Image Process. (2017) 26:2466–79. doi: 10.1109/TIP.2017.2672439

PubMed Abstract | CrossRef Full Text | Google Scholar