Magnetic resonance image reconstruction based on image decomposition constrained by total variation and tight frame

Abstract Objectives Magnetic resonance imaging (MRI) is a commonly used tool in clinical medicine, but it suffers from the disadvantage of slow imaging speed. To address this, we propose a novel MRI reconstruction algorithm based on image decomposition to realize accurate image reconstruction with undersampled k‐space data. Methods In our algorithm, the MR images to be recovered are split into cartoon and texture components utilizing image decomposition theory. Different sparse transform constraints are applied to each component based on their morphological structure characteristics. The total variation transform constraint is used for the smooth cartoon component, while the L 0 norm constraint of tight frame redundant transform is used for the oscillatory texture component. Finally, an alternating iterative minimization approach is adopted to complete the reconstruction. Results Numerous numerical experiments are conducted on several MR images and the results consistently show that, compared with the existing classical compressed sensing algorithm, our algorithm significantly improves the peak signal‐to‐noise ratio of the reconstructed images and preserves more image details. Conclusions Our algorithm harnesses the sparse characteristics of different image components to reconstruct MR images accurately with highly undersampled data. It can greatly accelerate MRI speed and be extended to other imaging reconstruction fields.


INTRODUCTION
MRI has been widely used in clinical medicine due to its excellent imaging capabilities of anatomical structures and physiological functions without ionizing radiation. 1,2owever, because conventional MRI requires sequential acquisition of the entire k-space data, the speed of MRI is slow, which affects clinical throughput and image The design of sparse transforms of MR images is a crucial factor that significantly impacts the quality of CS reconstruction.Researchers have done a lot of work on constructing sparse transforms of MR images.The total variation (TV) transform used the sparsity of images in the gradient domain to perform MR reconstruction, which preserved the edges of images but was prone to resulting in blocky effects. 5,6Orthogonal transforms, such as orthogonal wavelet transforms and curvelet transforms, are widely used in CS because orthogonality is beneficial to theoretical analysis and algorithm design. 7,80][11][12][13] This is because redundant atoms in TF can make the transformed coefficients sparser and better capture different image features than orthogonal dictionaries do.Some representative TFs are patch-based methods [9][10][11] and shiftinvariant wavelet frames. 12,13Conventional patch-based methods learned adaptive dictionaries from example image patches to alleviate the drawback of using the fixed basis but might lose important features due to the subdivision of patches. 9To make up for this deficiency, a convolutional sparse coding method was proposed to represent the image with a summation of a set of convolutions of the feature maps. 10,11With shift-invariant wavelet frames, a projected fast iterative soft-thresholding algorithm (pFISTA) was proposed to quickly solve the analysis model and acquire accurate MR reconstruction images. 12,13However, the studies mentioned above use a single or mixed set of sparse transforms for the entire image or image patches without accounting for the respective sparsity characteristics of different morphological components of the image, such as smooth parts and texture parts.Although the lowrank plus sparse model decomposed a full data matrix into two parts, it relied on the spatiotemporal correlations of multiple sequential images, making it suitable for dynamic MRI reconstruction rather than single-slice MR image reconstruction. 14,15n order to better analyze and process images, image decomposition theory has been employed in image denoising, 16 image fusion, 17 texture removal, 18 etc.It decomposes the image into the cartoon and texture components.The cartoon component comprises the smoother structural features of the image, which is usually constrained by the TV transform.The texture component contains the detailed information and noise of the image, which is constrained by the L 1 norm of sparse transforms. 16By applying distinct sparse transform constraints to different components, image decomposition can effectively represent image features.
Capitalizing on the benefits of image decomposition, the contribution of this paper is that we have proposed a novel MRI reconstruction algorithm based on image decomposition constrained by total varia-tion and L 0 norm of tight frame.It can simultaneously handle both image decomposition into two components and MR image reconstruction.By employing different sparse transforms tailored to each individual component, our method can effectively harness the sparse characteristics of different image components and improve MR image reconstruction performance with highly undersampled data.Specifically, the MR images to be reconstructed were divided into the cartoon and texture components.Since TF transform has superior performance in representing image detail features than orthogonal transforms and L 0 norm can obtain sparser solutions than L 1 norm, 19 we table introduced the L 0 norm constraint of TF transform to the texture component.For the cartoon component, we adopted the TV transform constraint to represent smooth features.The experimental results under various imaging conditions have demonstrated that our proposed algorithm outperforms two state-of -the-art CS MRI reconstruction methods.Our idea could also be extended to other CS reconstruction fields to improve their performance.In the following sections of the paper, we will use the term ''TV+TFL 0 '' to denote our proposed MRI reconstruction algorithm, which is based on image decomposition constrained by the TV and L 0 norm of TF.
The paper is organized as follows: first, the background of CS MRI and image decomposition theory is presented.Then our proposed TV+TFL 0 algorithm and its solving process are described in the Methods section.In the Results section, numerous experiments are conducted and the experimental results demonstrate the performance of our proposed algorithm.Conclusions and future work are given in the last section.

CS MRI
The undersampled k-space measurements y ∈ C M can be briefly described as follows: where x ∈ C N is the MR image to be reconstructed, F u ∈ C M×N (M < N) denotes the undersampled discrete Fourier transform operator, M and N stand for the number of the k-space measurements and the image pixels,  ∈ C M is the additive noise.Since M < N, undersampling improves the imaging speed of MR,but it also causes Equation (1) to be underdetermined.CS MRI leverages the sparsity of image transforms to solve this problem.The reconstruction model is expressed as follows 20 : (2) where Ψ denotes the sparse transform, ‖ ⋅ ‖ 0 is the L 0 norm constraint, λ serves as the regularization parameter to balance sparsity and data fidelity.Since the L 0 norm is nonconvex, which is difficult to solve, it is replaced with the L 1 norm which can be solved efficiently as a convex optimization 20 : arg min (3)

Image decomposition
Image decomposition theory states that images can be decomposed into cartoon and texture components to represent different features of the images 17 : where x c and x t represent the cartoon component and the texture component of the image x respectively.x c is a geometric and smoothly-varying component, while x t is an oscillatory component that contains details and noise.The cartoon-texture decomposition is commonly obtained through the TV+L 1 model which consists of a TV constraint and an L 1 norm constraint of sparse transforms 16 : arg min where ‖∇ x c ‖ 1 is the TV constraint, ‖Ψx t ‖ 1 is the L 1 norm constraint of a sparse transform Ψ,  c , and  t are the regularization coefficients of the cartoon component and the texture component, respectively.

Proposed TV+TFL 0 algorithm
In order to better represent the features of images sparsely and further improve the MR image reconstruction quality, we applied image decomposition theory to CS MRI.Building upon the image decomposition framework in Equation ( 5), we propose a new MR image reconstruction model named TV+TFL 0 : where the last term represents the MRI data fidelity term,  is the regularization coefficient for the data fidelity term, and the preceding three terms represent the image decomposition.Despite its greater computational challenge, the L 0 norm can obtain sparser solutions compared to the L 1 norm. 19Hence the texture component by image decomposition in this paper is constrained using the L 0 norm.Besides, the sparse transform Ψ adopts a tight frame for better capturing image details, satisfying Ψ * Ψ = I, where the superscript * denotes the adjoint operator, and I is the identity matrix. 12

Solving process
To solve the optimization problem of Equation ( 6), we employ an alternating iterative minimization approach, wherein one variable is optimized at a time while keeping other variables fixed.Specifically, the solving process is divided into three steps: Firstly, the cartoon component x c is optimized via Equation ( 7) with the MR image to be reconstructed x and the texture component x t fixed: Equation ( 7) can be solved using the nonlinear conjugate gradient method. 21econdly, we tackle x t with x and x c fixed for solving Equation ( 6) to yield: Equation ( 8) is the L 0 -norm non-convex optimization problem, and the iterative hard thresholding algorithm (IHTA) has been demonstrated to be a highly effective solution method. 22Besides, Ψ in Equation ( 8) is a tight frame, and paper 12 has proposed a fast algorithm called pFISTA to solve the tight frame based L 1 norm problem.Therefore, here we combine IHTA and pFISTA, denoted as projected fast iterative hard-thresholding algorithm (pFIHTA), to solve Equation ( 8): where k represents the number of iterations, is the step size, and H  t (q) denotes the hard thresholding function: (1) Initialization: x 0 is obtained by Equation ( 13), x 0 t = 0 (2) While stop criteria are not met, do: Update x k c by using the nonlinear conjugate gradient method to solve Equation (7); Update x k t by using pFIHTA to solve Equation ( 8); Update x k+1 by using Equation (12) to solve Equation (11); End Thirdly, we get x via the following formulation according to Equation ( 6) with x c and x t fixed: Equation ( 11) is the least squares problem and its analytic solution can be obtained as follows: ) where F represents the full Fourier transform matrix normalized, the superscript H denotes the Hermitian transpose operation, the matrix FF H u F u F H is a diagonal matrix consisting of ones and zeros, and the ones correspond to the sampled locations in k-space. 20n the first step of the above-mentioned process, x needs to be assigned an initial value.The initial value is obtained by a tight frame based L 0 norm optimization problem as shown in Equation (2) when Ψ is a tight frame.Then Equation (2) can also be solved by pFIHTA: In summary, the solving steps of the MRI reconstruction model based on TV+TFL 0 image decomposition proposed in this paper are illustrated in Table 1.

Experiments
The performance of the proposed algorithm has been evaluated by recovering MR images using a variety of sampling schemes at different undersampling factors, with and without noise.Three MR images used in the experiments (the 7th and the 27th slice of the brain MRI 1 23 and the 31st slice of the brain MRI 2 9 ) are shown in Figure 1.The size of all experimental images is 256 × 256.The MR images were normalized to a maximum magnitude of 1. Similar to prior work on CS MRI, 9,20 the CS data acquisition was simulated by undersampling the 2D discrete Fourier transform of the MR images.The sampling schemes used included 2D variable density random sampling, 1D variable density random sampling and pseudo radial sampling, as depicted in Figure 2. 2D variable density random sampling is usually considered as an ideal case for algorithm verification and can also be used for the 2D phase encodings in 3D imaging in practical applications. 12,21In the noisy case, zero-mean complex white Gaussian noise was added in k-space with the standard deviation = 0.01 or 0.02, which is widely used for noisy model in MRI. 12,24Our proposed method TV+TFL 0 was compared with two state-of -theart CS MRI algorithms, MCTV-L 2 6 and pFISTA. 12,13or MRI reconstruction, MCTV-L 2 used a non-convex isotropic TV constraint and pFISTA adopted a tight frame based L 1 norm constraint.
In all MR reconstruction experiments, the parameters used in all algorithms were hand-tuned to achieve the best reconstruction performance in terms of PSNR.For the MCTV-L 2 algorithm, the regularization coefficient was set as  TV = 0.002, and the augmented Lagrangian parameter = 150.For the pFISTA algorithm, the regularization coefficient was  p = 0.001, and the step size was  p = 0.1.For the TV+TFL 0 algorithm proposed in this paper, the parameters were set as = 0.03,  1 = 1.3,  c = 0.004,  t = 0.01, = 1.1,  → +∞ in noiseless experiments and  = 1000 in noisy experiments.To compare fairly, the tight frames in pFISTA and TV+TFL 0 were the same, both using shift-invariant discrete wavelet transform (SIDWT).The shift-invariant property of SIDWT is beneficial in suppressing artifacts caused by the Gibbs phenomenon in discontinuous regions of the images. 12

Performance metrics
The quality of the reconstruction results is evaluated by the peak signal-to-noise ratio (PSNR) and a high frequency error norm (HFEN).PSNR and HFEN are two commonly used image quality metrics in CS MRI. 9,20,25SNR is defined as 25 : where I is the ground truth image, MAX I denotes the maximum possible value of the image, I rec is the reconstructed image, m and n are the number of rows and columns of the image respectively.HFEN is used to assess the reconstruction quality of edges and fine features.It is calculated as the L 2 norm of the result obtained by filtering the difference between the reconstructed and ground truth images with a Laplacian of Gaussian filter. 20esides, the error images, which are the absolute value of the difference between the reconstructed images and their corresponding ground truth images in Figure 1, are employed as a visual metric for evaluating the reconstruction quality. 9

PSNR and HFEN comparison
The PSNR of the reconstruction results of MCTV-L 2 , pFISTA and TV+TFL 0 under different undersampling rates and sampling schemes with and without noise are presented in Table 2.It can be observed from Table 2 that the proposed TV+TFL 0 in this paper greatly improved PSNR compared to the other two algorithms under various imaging conditions for all three test images.Using the imaging condition with the noise stan-dard deviation  = 0.01, the undersampling rate R = 0.2 and pseudo radial sampling as an example, the PSNR of the 7th slice of the brain MRI reconstructed by MCTV-L 2 , pFISTA, and TV+TFL 0 were 29.38, 29.74, and 30.75 dB, respectively (as shown in Table 2).Our method yielded an enhancement in PSNR of more than 1 dB.According to Table 3, TV+TFL 0 exhibits lower HFEN values than MCTV-L 2 and pFISTA in most cases, indicating its superior capability in capturing edges and fine features.Although pFISTA produces slightly higher PSNR than MCTV-L 2 in some cases, it generally has the worst HFEN among the three methods.Overall, TV+TFL 0 outperforms the other two algorithms in terms of PSNR and HFEN.
To further verify the performance of TV+TFL 0 at different undersampling rates, Figure 3  reconstruction performance at various undersampling rates.

Visual inspection comparison
To visually assess the image reconstruction quality, the three experimental MR images recovered by MCTV-   our TV+TFL 0 algorithm shown in Figures 4-6f , indicating that pFISTA lost more features.The reason is that MCTV-L 2 and pFISTA only exploit the sparsity in the finite difference domain and the SIDWT domain respectively, which is insufficient to capture the complex structures within the images.As shown in Figures 4-6c,f , in contrast to MCTV-L 2 and pFISTA, our TV+TFL 0 method successfully suppressed more artifacts and retained more detail information by utilizing two more suitable sparse transform constraints based on image decomposition theory.Under other imaging conditions as listed in Table 2, compared with the other two algorithms, TV+TFL 0 also obtained superior reconstructed images, which is similar to the example shown.

Computation cost comparison
For comparing the computation cost of the three algorithms, we conducted the experiment on reconstructing Figure 1c by using 1D variable density random sampling with the undersampling rate R = 0.2 and the noise standard deviation  = 0.01.The experiment was implemented in Matlab 2019b on a PC with a 3.10 GHz Intel Core i5 CPU and 8.0 GB memory.The CPU running time of MCTV-L 2 ,pFISTA,and TV+TFL 0 are 6.35,64.47,and 79.29 s, respectively.As shown in Tables 2 and 3, although TV+TFL 0 takes a little longer time than that of MCTV-L 2 and pFISTA, TV+TFL 0 can achieve the best reconstruction performance of Figure 1c.

CONCLUSION
In this paper, we have developed a novel MRI reconstruction model based on TV+TFL 0 image decomposition for better sparse representation of MR image features.According to image decomposition theory, we applied distinct sparse transform constraints to the cartoon and texture components according to their respective distinct structural characteristics.For the cartoon component, the TV transform constraint was employed to represent smooth features.For the texture component, the L 0 norm constraint based on tight frame redundant transform was utilized for capturing detailed features.Compared with MCTV-L 2 and pFISTA algorithms, the results of various experiments have demonstrated that our proposed TV+TFL 0 achieved superior performance for reconstructing MR images by both quantitative and qualitative visual evaluations, such as improved PSNR and HFEN, enhanced detail clarity and reduced artifacts.However, due to the cartoontexture image decomposition and their respective different sparse transform constraints, the computational complexity and parameters that need to be tuned increase.In the future, we will study how to speed up sparse transform learning and automatic regularization parameter selection method for the proposed algorithm.

AU T H O R C O N T R I B U T I O N S
Guohe Wang: methodology, software, and writing; Xi Zhang: data analyses; Li Guo: conceptualization and supervision.All authors reviewed and edited the manuscript.

AC K N OW L E D G M E N T S This work was supported by the Tianjin Municipal Education Commission Research Program (2019KJ181).
The authors sincerely thank the editors and reviewers for their constructive comments.

C O N F L I C T O F I N T E R E S T S TAT E M E N T
The authors declare no conflicts of interest.

F I G U R E 1
The experimental data of MRI.The parts of figure (a) and (b) are the 7th and the 27th slice of the brain MRI 1.The part (c) is the 31st slice of the brain MRI 2. F I G U R E 2 Sampling schemes.(a) 2D variable density random sampling.(b) 1D variable density random sampling.(c) Pseudo radial sampling.
displays how the PSNR and HFEN values of the 27th slice of the brain MRI 1 reconstructed by the three algorithms change with the undersampling rate from 0.2 to 0.7 under 1D variable density random sampling scheme without noise.From Figure 3 it can be clearly observed that the PSNR and HFEN values of TV+TFL 0 are better than those of the other two compared algorithms.Hence, we can deduce that TV+TFL 0 has superior MRI F I G U R E 3 The parts of the figure (a) and (b) are PSNR and HFEN curves of the 27th slice of the brain MRI 1 recovered by MCTV-L 2 , pFISTA, and TV+TFL 0 at different undersampling rates.F I G U R E 4 Reconstruction results of the 7th slice of the brain MRI 1.The parts of the figure (a)-(c) are the reconstruction results of MCTV-L 2 , pFISTA, and TV+TFL 0 .The parts (d)-(f) are corresponding error images of (a)-(c).

F I G U R E 6
Reconstruction results of the 31st slice of the brain MRI 2. The parts of the figure (a)-(c) are the reconstruction results of MCTV-L 2 , pFISTA, and TV+TFL 0 .The parts (d)-(f) are corresponding error images of (a)-(c).
The solving steps of the MRI reconstruction model based on TV+TFL 0 image decomposition.
TA B L E 1 Input:y, ,  1 ,  c ,  t , , and  Output: reconstructed MR image x 3 Note:The lowest HFEN value is highlighted in bold.