Adaptive wavelet tight frame construction for accelerating MRI reconstruction

The sparsity regularization approach, which assumes that the image of interest is likely to have sparse representation in some transform domain, has been an active research area in image processing and medical image reconstruction. Although various sparsifying transforms have been used in medical image reconstruction such as wavelet, contourlet, and total variation (TV) etc., the efficiency of these transforms typically rely on the special structure of the underlying image. A better way to address this issue is to develop an overcomplete dictionary from the input data in order to get a better sparsifying transform for the underlying image. However, the general overcomplete dictionaries do not satisfy the so-called perfect reconstruction property which ensures that the given signal can be perfectly represented by its canonical coefficients in a manner similar to orthonormal bases, resulting in time consuming in the iterative image reconstruction. This work is to develop an adaptive wavelet tight frame method for magnetic resonance image reconstruction. The proposed scheme incorporates the adaptive wavelet tight frame approach into the magnetic resonance image reconstruction by solving a l0 regularized minimization problem. Numerical results show that the proposed approach provides significant time savings as compared to the over-complete dictionary based methods with comparable performance in terms of both peak signal-tonoise ratio and subjective visual quality.


Introduction
Compressed sensing (CS) [1,2] aims to reconstruct a signal, which is assumed to be sparse in a particular transform domain, from significantly fewer measurements than mandated by traditional Nyquist sampling. Magnetic resonance imaging (MRI) is a versatile and advanced technique for non-invasive assessment of human tissue with a comparatively slow data acquisition process. Applying CS to MRI is an effective way to overcome the slow data acquisition problem, and some known results from it have demonstrated that CS has great potential in both scan time reduction and reconstruction accuracy improvement. The CS-based methods reconstruct a magnetic resonance (MR) image by solving a nonlinear optimization problem that imposes sparsity constraints on the MR image. The sparsity is achieved under a carefully selected transform system which is crucial for successful recovery in CS.
Although various sparsifying transforms have been used in MRI reconstruction such as wavelet [3], contourlet [4],and total variation (TV) [5] etc, the efficiency of these transforms typically rely on the special structure of the underlying image. For example, the TV transform favors piecewise constant images, while contourlet transform

Preliminaries
In this section, we first give a brief introduction to wavelet tight frame and its construction derived from a multiresolution analysis. We only consider the un-decimated wavelet tight frames, as they are known to perform better than the decimated wavelet systems in most cases [16,17].

Tight frame
A frame is a generalization of Riesz basis and tight frame can be viewed as a redundant system that generalizes orthonormal basis. Let H be a Hilbert space holds for any f ∈ H.This condition is equivalent to the so-called perfect reconstruction property f = ∑ n ⟨f, x n ⟩ x n , f or all f ∈ H (2) There are two operators associated with the tight frame, i.e., analysis operator and its adjoint called a synthesis operator. The analysis operator T decomposes a vector f ∈ H into the canonical coefficients under the tight frame, 202 ADAPTIVE WAVELET TIGHT FRAME CONSTRUCTION FOR ACCELERATING MRI RECONSTRUCTION which is defined by And the synthesis operator satisfies The operator W = T * T is called the frame operator. Then the condition (2) shows that W = I holds if and only if {x n } ⊂ H is a tight frame, where I : H → H is the identical operator. It should be noted that T T * is usually not identical operator.

Wavelet tight frame for finite dimension
The wavelet tight frame is a tight frame generated by a bank of where h 0 is the so-called low-pass filter and h i , i = 1, 2, · · · , r − 1 are the so-called high-pass filters. Given a signal g ∈ R N , the one-level un-decimated wavelet transform with respect to the filters h is defined as follows where * denotes the filtering procedure. With a slight abuse of notation, we also use W (h) to denote the corresponding operator matrix. It is easy to show that when the symmetric boundary condition of is employed, then the wavelet transform is a Toeplitz-plus-Hankel matrix.
In order to make this wavelet transform W (h) to be a tight frame, the equation W (h) * · W (h) = I should be hold, where (·) * denotes the adjoint operator. Following the previous work, we know that this perfect reconstruction condition can be equivalently stated as one of the full Unitary Extension Principle (UEP) condition [18,19], which has proven to be a highly useful methodology for constructing a wavelet tight frames with a structure of multiresolution analysis. The UEP condition is stated as follows.
where δ k = 1 if k = 0 and δ k = 0 otherwise. A family of wavelet tight frame is constructed by using the above UEP condition [20,21], e.g., the piecewise linear B-spline framelet can be constructed by the associated filters It can be easily shown that h = {h 0 , h 1 , h 2 } satisfies the condition (4), so that W (h) forms a tight frame with wavelet structure. For 2D images, the tensor product method is often used to form a 2D wavelet tight frames. Suppose that h = {h i } r−1 i=0 satisfies the condition (4), then W (h) forms a tight frame for 1D signal in R N . The 2D wavelet tight frames can be constructed by W (h) ⊗ W (h) , which is a tight frame for signal in R N 2 generated by concatenating all columns of a 2D signal in R N ×N . The tensor products between every pair of filters in can also equivalently construct the transform operator W (h) , where To get a multi-level wavelet transform, we can recursively using the one-level wavelet transform W (h) to the low-pass coefficients. We omit the detailed discussion, and interested reader may refer to [22,23].

Adaptive wavelet tight frame construction
As mentioned above, due to the various structure of practical images, one predefined 2D/1D sparsifying transform working well for one type of images may not work for another. To address this issue, adaptive transform such as learned dictionary has been introduced to provide a better sparse approximation to the underlying image to be recovered. A representative dictionary learning method is K-SVD [10], which gives the learned dictionary by solving following minimization problem Therein, ∥·∥ 2 denotes ℓ 2 norm and the ℓ 0 pseudo-norm ∥·∥ 0 accounts for the number of nonzero entries. The training data y i is may be acquired by some available images or overlap patches extracted from a previous reconstructed image. Solving problem (5) gives the learned (overcomplete) dictionary under which each example y i has sparse representation α i . The regularization parameter λ balances a trade-off between data fidelity and sparsity. The K-SVD method is a very popular method that has been applied to a wide variety of image processing applications and medical image reconstruction [11,12,24,25]. But the problem (5) is very challenging due to the extreme ill-posedness of this problem. An alternating scheme introduced in [13] is proposed by imposing additional orthonormal constraint on the designed dictionary. Owing to this constraint, the orthogonal dictionary learning-based algorithm for image processing [14] and MRI reconstruction [13]provides significant time saving and competitive result as compared K-SVD-based methods. Different from these local patch-based transform, the wavelet system is usually considered to be a global transformation which has been shown disadvantage in a variety of applications such as image denoising, image inpainting and image reconstruction, etc. In fact, if the wavelet filters have finite support, which is the case in many practical situation, the global wavelet transform can also be considered as a local patch-based transform due to the Toeplitz structure of the resulting transform matrix. Therefore, the worse performance of the wavelet transform-based methods compared other patch-based approaches arise from that the transform is predefined and fixed, rather than globality of the transform system. Motivated by this fact, Cai et. al [15]constructed an adaptive tight wavelet frame by solving the following balanced minimization problem with orthogonal constraints: where W (h) is the wavelet tight frame associated with a set of filters h , (·) T denotes the transpose operator and is replaced by (·) H (conjugate transpose) for complex-valued matrix, y ∈ R N is the measured data, and α is a sparse vector which is close to the canonical coefficients of y under the system Through an iterative scheme for solving (6), the wavelet tight frame can be obtained by alternatively updating the coefficient vector α and the filters h = {h i } r−1 i=0 which generates a tight frame as in (3), resulting two subproblems. The first subproblem, i.e., solving (6) with respect to α given the frame filter h , can be efficiently solved by soft hard thresholding. But the second subproblem, i.e., solving (6) with fixed , in general is a complex minimization problem. By assuming a special case, where ⟨h i , h j ⟩ = 1 / r 2δ i−j ,i, j = 0, 1, · · · , r − 1 , then h can be easily obtained by the singular value decomposition.

MRI reconstruction using adaptive wavelet tight frame
There are two ways to employ an adaptive wavelet tight frame for specific task in image processing such as image denoising and inpainting, and medical image reconstruction. One is that, firstly, constructing the adaptive wavelet tight frame from a reference image, and then using it as a predefined transform to the specific task (e.g., [15,26,27]). For single-coil MRI reconstruction, under the adaptive wavelet tight frame obtained by solving (6), the underlying MR image can be reconstructed from undersampled k-space data by solving the following formulation: where x is the MR image to be reconstructed, y denotes the single-coil undersampled k-space data with at some undersampling factor, F u is the undersampling Fourier operator, and W is the learned adaptive tight frame with wavelet structures. Another way for utilizing the adaptive wavelet tight frame is that incorporating the model (6) and (7). In the case of MRI reconstruction, this resulting in the following balanced sparse optimization problem with ℓ 0 regularization: where we use (·) H instead of (·) T because the MR data are complex-valued.
In this paper, due to the difficulty in choosing the reference MR image, we consider the second way to reconstruct the MR image by solving the problem (8). Starting with the initial r 2 filters h = {h i } r 2 −1 i=0 with support on Z 2 ∩ [0, r − 1] 2 for some positive integer , we take an alternating iterative scheme to alternatively update the estimations of the sparse vector , the underlying image initialized by zero-filling reconstruction, and the filters h . Specifically, we solve the ℓ 0 regularization problem (8) by alternatively minimizing over one variable with others fixed, resulting three subproblems: and The solution to the first subproblem is known to have a unique solution α * obtained by applying a hard thresholding operator on the canonical frame coefficients W (h)x: where H √ λ (·) is an element-wise hard thresholding operator defined as The subproblem (10) is a least squares problem whose normal equation is equivalent to Solving the above equation directly is impractical because it requires to invert a N 2 × N 2 matrix premultiplying x for an image size of N × N . Fortunately, following the way described in [11], we can effectively obtain the Fourier transform of the image x instead of x itself. We omit the details here and directly present the solution of equation (14). Let S 0 (k), where k is the position index in the k-space, S 0 (k) denotes the zero-filled version of the measurements y , and is the Fourier transform of the reconstructed image W H α with the modified coefficients α , i.e., S 0 (k) = F F u H y and S 1 (k) = F W H α where F is the standard discrete Fourier transform matrix. Then the Fourier transform of x, denoted as S(k) , can be obtained by where Ω denotes indices subset of k-space corresponding to sampled locations. The subproblem (11) is a constrained minimization problem with the quadratic constraints W H W = I . Due to the special wavelet structure of W , this minimization problem is rather complex in terms of implementation and solving it could be computationally expensive. To efficiently solve (11), we employ sufficient condition for the orhogonality of the wavelet transform which was proposed and proved in [15]. Let D i denotes a column vector obtained by column-wise concatenating h i , followed by flipping it in the up-down direction, and D is a matrix generated by collecting all D i -s as columns. Then in order that W (h) H W (h) = I holds, we need only On the other side, the wavelet transform of an image can be implemented by first extracting overlapping patches of size r × r in the image, and second premultipling by the matrix D . Therefore, the subproblem (11) can be restated as where X is a matrix whose i-th column X i is the vectorized version of the patches extracted from ,x and A is a rearranged version of α where each column is correlated to DX i (9). The following proposing gives the exact solution to (16), or equivalently to (11).
is the singular value matrix of AX H , U andV are the collection of left and right singular vectors, respectively. Then, D * = 1 / r U V H is the optimal solution of the constrained optimization problem (16). This solution is unique if and only if all singular values of AX H are nondegenerate and nonzero.

Algorithm 1. Adaptive wavelet tight frame construction for MRI reconstruction (ATFMRI)
Input: k-space measurements y Output: Reconstructed image x Initialization:

1) Initialize tight frame filters
using some existing tight frame; 2) Initialize the reconstruction using zero-filled measurements: x (0) = F u H y . Iteration: For the k-th iteration 2) Apply hard-thresholding to the wavelet coefficients via (12) and (13), obtaining α (k) ; 3) Update k-space data using Eq. (15), and then apply inverse FFT to the updated k-space data, resulting in updated reconstruction x (k) ;

4) Compute the tight frame filters
according to the Proposition 1.
4a) Construct matrix X whose columns consist of all concatenated patches extracted from x (k) , and construct A by rearranging the entries of α (k) according to X ; 4b) Compute the SVD of AX H , i.e., AX H = U ∑ V H ; 4c) Set h i (k) to be the i-th column vector of the matrix D = 1 / r U V H .
The complete procedure of the numerical solver for solving (8) is summarized in Algorithm 1. In our implementation we stop the algorithm when the relative change of the norm of the reconstructed image between two 206 ADAPTIVE WAVELET TIGHT FRAME CONSTRUCTION FOR ACCELERATING MRI RECONSTRUCTION consecutive iterations is smaller than a certain threshold. As each subproblem in the alternating iteration above is exactly solved with closed-form optimal solutions, the cost functional in (8) is bound below and monotonically decreasing, resulting in convergence of the sequence of the cost functional. However, due to the nonconvex objective function in (8), the convergence of iterates in the algorithm cannot be theoretically guaranteed. To yield a convergence result, several strategies have been adopted such as replacing the ℓ 0 pseudo-norm by ℓ 1 norm [28,29,30], and gradually decreasing the regularization parameters during iterations [13], etc.

Numerical experiments
In this section, we present experiments on several in-vivo MR scans of size 512 × 512 which are available at American Radiology Services [31]. The proposed ATFMRI is compared with two representative algorithms on recovering MR images from the undersampled k-space data. One of these was proposed by Lustig et al. [3], denoted as sparseMRI, which utilizes non-adaptive wavelet and TV as global sparsifying transform. Another is patch-based overcomplete dictionary learning method for MRI reconstruction, denoted as DLMRI, proposed by Ravishankar et al. [11]. In our experiments, we retrospectively undersampled data with varying undersampling factors. Two metrics are used for the quantitative measurements: the peak signal-to-noise ratio (PSNR) and a high frequency error norm (HFEN). The PSNR in decibels (dB) is calculated by: where MSE is the mean squared error between the reference and the reconstructed images. The HFEN is used to quantify the quality of reconstruction of edges and fine features and pixels, and is computed as the ℓ 2 -norm of the result obtained by filtering the difference between the reconstruction and reference image with Laplacian of Gaussian filter.
For fair comparison, we use the publicly available code for both sparseMRI and DLMRI, and all methods use the same simulated noise and sampling pattern for the k-space. Furthermore, for each compared method, we set the main parameters by sweeping through a set of values and selecting the best ones for comparison in terms of PSNR. For the proposed method, the main parameters are fixed and empirically set as: the filter size is 5 × 5 , the regularization parameters are set to λ = 0.08 and µ = 20/σ respectively, where σ denotes the standard deviation of the measurement noise. In our implementation, we considers the measurement noise has exact noise level or can be estimated when the exact value is unknown. All the experiments were implemented using Matlab 2011b on an operation system with 32-bit Windows 7, with a core 3.7 GHz Intel processor and 4 GB of memory.

Simulations on different MR scans
Two experimental evaluations are carried out to study the performance of our proposed algorithms. In this subsection, the four reference images (see Fig. 1) are transformed to k-space and then retrospectively undersampled at a moderate acceleration factors of 4 with Cartesian sampling patterns. The noiseless scenario (i.e., σ =0 in the experiments) is considered in this evaluation. We will study the noise case with various undersampling factors in next subsection. Table 1 provides the associated quantitative results in terms of PSNR, HFNE and CPU time of the compared methods. It can be shown that the adaptive system-based methods yields significant benefits in improving reconstruction quality in comparison with the non-adaptive ones. Although the PSNR and HFEN values of the results from the two adaptive system-based methods are comparable, the proposed ATFMRI provides significant time saving as compared to DLMRI.
In Fig. 2, we show qualitative comparisons between these schemes using the spine image with Cartesian sampling (the sampling mask is showed in Fig. 2(b)) at acceleration factor 4. The zero-filled Fourier reconstruction (Fig. 2(c)) has significant undesirable artifacts due to undersampling. Some of which persist in the sparseMRI reconstruction (see Fig. 2(d)). This phenomenon can also be found in later Fig. 3. In contrast, the results with adaptive transformation (Fig. 2 (e) and (f)), including learned dictionary and learned tight frame, are seen to be substantially less artifacts and more close to the reference image ( Fig. 2(a)). Fig. 2 (g) and (h) give quantitative performance of the proposed method and experimental evidence that the proposed ATFMRI algorithm is convergence for this dataset. The regions of interest (ROI), as indicated by a white rectangular in Fig. 2(a), are enlarged and displayed at the bottom-right corner of each image. It can also be seen from the enlarged images that the DLMRI and the ATFMRI scheme are more effective than the non-adaptive ones in reducing the over-smoothing of the reconstructed images.

Performance with various sampling rate
In this subsection, we study experimental evaluation under k-space Gaussian noise of standard deviation 16. Table 2 lists the reconstruction results of the compared methods under various acceleration factors. As a special case, when the acceleration factor is 2, the sparseMRI has similar performance with the adaptive transform-based methods. However, when the acceleration factor is more than 2, the adaptive transform-based methods have major advantages compared to the non-adaptive one. On the other hand, it can be seen that the sparseMRI and ATFMRI methods greatly improved the computational efficiency compared with DLMRI. Fig. 3(a-c) are the example reconstruction of k-space data with 4-fold acceleration and noise level using sparseMRI, DLMRI and ATFMRI respectively. Fig. 3(d-f) are their difference images from the reference image respectively. It can be seen from the difference images that the DLMRI and the ATFMRI scheme are more effective than the non-adaptive ones in improving PSNR and HFEN of the reconstructed images. The corresponding quantitative values are listed in Table 2

Conclusion and future work
In this paper, we first present an optimization model incorporating adaptive wavelet tight frame into MRI reconstruction, and then design an alternating scheme to solve the problem. Different from previous methods, which first drive a discrete wavelet tight frame system from the reference images and then use this constructed  adaptive tight frame for various image processing tasks, we iteratively update the tight frame using previous reconstruction. Based on the unitary extension principle condition, we have developed an efficient method to simultaneously construct an adaptive tight frame and reconstruct MR image from the undersampled k-space data. Numerical evaluations suggest that the adaptive transform-based methods have major advantages compared to the non-adaptive one. On the other hand, although the quantitative performance of the two adaptive system-based methods are comparable, the proposed ATFMRI provides significant time saving as compared to DLMRI.
In the future work, we would like to put efforts on two directions. One is that constructing an adaptive tight frame under more general conditions using the UEP. Also, we will develop a more efficient and convergent algorithm for solving the proposed model.