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.

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)(6)(7)(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)(15)(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. 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 higherorder 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 noisereduction effectiveness was obtained. However, the HOSVDbased 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 ∈R J 1 ×J 2 ×...×J N signified an arbitrary Nth-order tensor over the multidimensional space of size J 1 × J 2 × . . . × 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 ...j N using (j 1 j 2 . . . 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 nth 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 ∈R J 1 ×J 2 ×...×J N would be represented as X (n) ∈ R J n ×(J 1 ×...×J n−1 ×J n+1 ×...×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 ∈R J 1 ×J 2 ×...×J N and a matrix V ∈ R K×J n along the nth dimension is a tensor of size J 1 × . . . × J n−1 × K × J n+1 × . . . × J N , denoted as X × n V, with its element expressed by following equation: (1) Definition 4 (Multilinear tensor rank). The multilinear tensor rank of X ∈R J 1 ×J 2 ×...×J N is defined as the rank of the mode-n unfolding matrix X (n) , denoted as rank n (X ).

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: where rank(·) represented the rank of matrix, referring to the number of nonzero singular values of matrix. λ was the tradeoff 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: where ||X|| * = i δ i (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.,

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 ∈ R p×p×p×m . Then the problem estimating the noise-free version X according to The bold values in represented the maximum values across the methods.
where rank 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: where S ∈R r 1 ×r 2 ×r 3 ×r 4 was the core tensor and U (i) ∈ R p×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: where rank(X (n) ) was the number of nonzero singular values of the mode-n unfolding matrix X (n) , and | |Y − X | | 2 F = | Y (n) − X (n) | 2 F 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: 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. 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 h(δ i (X), a). Since both the linear function tr(XY T ) and trace norm | |X| | * were convex function and the function tr(Y T Y) was independent on X, the J(X) was strictly convex when the function J 1 (X) was strictly convex. According to Parekh and Selesnick (36), we can conclude that the function J 1 (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: and the equivalent threshold function of the proximal operator (8) can be defined as (34,38) 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: where Y (n) = U (n) (n) (V (n) ) T . Therefore, the multilinear tensor rank parameter r n can be adaptively estimated as the number of the nonzero elements in θ (n) ; λ, a , and the estimated factorization matrixÛ (n) was selected as the first r n columns of the U (n) . Thus, the estimated core tensor can be formulated as: (1) ) T × 2(Û (2) ) T × 3(Û (3) ) T × 4(Û (4) ) T , (11) and the restored tensorX can be obtained by: ) T × 2(Û (2) ) T × 3(Û (3) ) T × 4(Û (4) )T. (12) 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 −1 ) 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: where Y is the 3D MR images contaminated by Riciandistributed 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 3 × 1 mm 3 × 1 mm 3 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 Algorithm 1 : 3D MR image denoising Input: 3D MR noise images Y noise Output: Denoised MRI image X Initialization: X =Y for t = 1,2,. . . T do Iterative regularization: Search for similar patches and group them to form fourthorder tensor Y for each reference patch; Decompose each tensor by Eq. (5) and compute θ (n) ; λ, a by Eq. (9); Update the factorization matrix U (i) (i = 1, 2, 3, 4); Compute the core tensor S by Eq. (11); Compute the denoised tensorX by Eq. (12); Aggregating allX to obtain denoised image; End for the denoising performance. PSNR was a metric measuring the ratio of the maximum possible signal power to the noise power, which is defined as: 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: where µ x and σ x denote the mean and standard deviation of the ground truth image, respectively, while the µ x * and σ x * denote the mean and standard deviation of the denoised image, respectively. c 1 and c 2 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 3 m), respectively.

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 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. 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.
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.

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 3 × 1.0 mm 3 × 1.25 mm 3 . 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.

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.oasisbrains.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.