High Resolution Time Resolved Contrast Enhanced MR Angiography Using k-t FOCUSS

Hong Jung; Eung Yeop Kim; Jong Chul Ye

doi:10.13104/jksmrm.2010.14.1.10

Abstract

Purpose

Recently, the Recon Challenge at the 2009 ISMRM workshop on Data Sampling and Image Reconstruction at Sedona, Arizona was held to evaluate feasibility of highly accelerated acquisition of time resolved contrast enhanced MR angiography. This paper provides the step-by-step description of the winning results of k-t FOCUSS in this competition.

Materials and Methods

In previous works, we proved that k-t FOCUSS algorithm successfully solves the compressed sensing problem even for less sparse cardiac cine applications. Therefore, using k-t FOCUSS, very accurate time resolved contrast enhanced MR angiography can be reconstructed. Accelerated radial trajectory data were synthetized from X-ray cerebral angiography images and provided by the organizing committee, and radiologists double blindly evaluated each reconstruction result with respect to the ground-truth data.

Results

The reconstructed results at various acceleration factors demonstrate that each components of compressed sensing, such as sparsifying transform and incoherent sampling patterns, etc can have profound effects on the final reconstruction results.

Conclusion

From reconstructed results, we see that the compressed sensing dynamic MR imaging algorithm, k-t FOCUSS enables high resolution time resolved contrast enhanced MR angiography.

Introduction

Angiography is an imaging technique to depict structures of vessels and blood flow in arteries and veins. The main interest of this technique is to find the disorder of blood vessels such as stenosis, occlusion, or aneurysms.

X-ray angiography is a standard tool for blood vessel diagnosis. For X-ray angiography, contrast agent must be injected into the patient using a catheter. Then, to obtain blood vessel images, digital subtraction angiography (DSA) is generally used which subtracts a subsequently acquired image without contrast agent from images with contrast agent. However, the procedure is quite invasive.

In order to resolve these issues, there have been intensive researches to exploit other imaging modalities such as CT or MRI in angiography. Currently, CT angiography (CTA) using iodine contrast agent is most widely used as a replacement of the conventional X-ray angiography (1). However, CTA is considered harmful for routine use due to radiation risk (2). In contrast, MR angiography (MRA) is safer (3). Contrast enhanced MRA (CE-MRA) uses gadolinium as a contrast agent. Gadolinium reduces T1 so that the signal intensity of blood vessel filled with gadolinium becomes high. However, it is generally accepted that the resolution of CE-MRA is poorer than that of CTA due to slow acquisition time. Since MR acquisition is based on Fourier transform, the number of measured k-space data should satisfy Nyquist sampling limit to avoid aliasing artifacts, which usually results in trade-off between spatial and temporal resolution in MRA.

Recently, several novel approaches have been developed to reduce scanning time while preserving image qualities. For example, parallel imaging methods such as sensitivity encoding (SENSE) (4), simultaneous acquisition of spatial harmonics (SMASH) (5), generalized autocalibtrating partially parallel acquisitions (GRAPPA) (6) accelerate data acquisition time by obtaining a subset of k-space data and resolve aliasing artifacts by exploiting the diversity of coil sensitivity. However, it was reported in (7, 8) that an acceptable image quality can be generated at only limited acceleration factor due to the amplification of noise.

Our approach is different from conventional algorithms thanks to the use of "compressed sensing" theory (9, 10). According to compressed sensing theory, perfect reconstruction is possible even from sampling rates dramatically smaller than the Nyquist sampling limit by solving an l₁ minimization problem, as long as the non-zero spectral signal is sparse and the samples are obtained with an incoherent basis (10). Therefore, compressed sensing may be one of the most suitable approaches for CE-MRA, since the temporal variation of CE-MRA is usually very sparse in spatio-temporal domain. Recently, we have developed compressed sensing based dynamic MR imaging algorithm called k-t FOCUSS (11, 12). k-t FOCUSS has been successfully applied to cardiac cine imaging. Hence, applying k-t FOCUSS to CE-MRA, great performance can be expected thanks to its spatio-temporal sparsity. Furthermore, by incorporating k-t FOCUSS in radial trajectory using golden trajectory (13), the advantage of radial trajectory can be exploited as well. In order to make k-t FOCUSS more effective, the temporal basis should be chosen judiciously such that the transformed signal can be effectively sparsified. We found that principal component analysis (PCA) provides a suitable basis. Experimental results using the data provided by Recon challenge at 2009 ISMRM workshop for Data Sampling and Image Reconstruction show that an accurate reconstruction can be achieved from very sparse measurements.

Even though the results presented in this communication is based on simulated MRA data from X-ray cerebral angiography, the insight we obtained from this exercise may provide a promising direction in CE-MRA.

Methods

k-t FOCUSS has been developed for high spatiotemporal resolution dynamic MR imaging such as cardiac cine imaging (11, 12). The basic concept for k-t FOCUSS comes from compressed sensing theory which tells us that the accurate reconstruction is possible using only limited number of measurements if the unknown image can be sparsely represented. Here, the essential requirement for the success of compressed sensing is the "sparsity" of images. Time resolved CE-MRA is therefore suitable for compressed sensing framework since the difference between adjacent frames is only restricted inside of vessels which are relatively stationary and sparse. Furthermore, the temporal variation of contrast agent is often slowly varying, which can be effectively sparsified in appropriate basis. According to compressed sensing theory, those very sparse signals can be very accurately reconstructed by solving a computationally feasible l₁ minimization problem as follows:

min ‖ρ‖₁

[1]

subject to ‖v- Fψρ‖₂ ≤ ε,

where ρ, ψ, F, and υ represent a vector stacked by unknown sparse signal, temporal sparsifying transform, Fourier transform from image space to k-space, and a vector stacked by k-t measurements, respectively. Here ε, denots the noise level. In order to solve Eq. [1], k-t FOCUSS finds:

[2]

ρ_n+1 = W_nq_n,

by iteratively solving re-weighted l₂ minimization problem as follows:

min ‖q_n‖₂

[3]

subject to ‖v- FψW_nq_n‖₂ ≤ ε,

where W_l+1 is updated as:

[4]

W_{n + 1} = (\begin{matrix} \begin{matrix} {|ρ_{n + 1} (1)|}^{p} & 0 \\ 0 & {|ρ_{n + 1} (2)|}^{p} \end{matrix} & \begin{matrix} \dots & 0 \\ \dots & 0 \end{matrix} \\ \begin{matrix} ⋮ & ⋮ \\ 0 & 0 \end{matrix} & \begin{matrix} ⋱ & ⋮ \\ \dots & {|ρ_{n + 1} (N)|}^{p} \end{matrix} \end{matrix}), 1 / 2 \leq p \leq 1 .

In (11), it was proved that this algorithm solves l₁ minimization problem by setting p = 0.5. If we set 0.5<p<1, Eq. [3] will solve l_q minimization for q<1, which guarantees a sparser global minimization. However, since l_q (q<1) is not a convex minimization, the solution can be fallen to local minimization. The k-t FOCUSS algorithm can be then represented as:

[5]

ρ_n+1=Θ_nΨ^HF^H(FΨΘ_nΨ^HF^H+λI)^-1υ,

where Θ_n=W_nW_n^H and λ is a regularization parameter. If we increase λ, Eq. [5] obtains sparser solution but the data fidelity term in Eq. [3] becomes less important. Usually, λ≪1 provides good image qualities.

Note that there exist several options to be optimized in k-t FOCUSS. For example, the specific choice of sparsifying transform, sampling patterns and etc. These issues will be discussed in detail in the sequel.

Temporal Sparsifying Transform

As mentioned in previous section, the "sparsity" of unknown signal is an essential condition for compressed sensing reconstruction. In (11, 12), Fourier transform (FT) along temporal direction is used to sparsify cardiac cine images. FT tends to sparsify a signal especially when signal shows sinusoidal pattern. Therefore, when the object has periodic motion (for example moving heart), temporal FT can effectively sparsify images.

However, even if the signal does not have an obvious periodic pattern, the Fourier spectrum is mostly focused at the low frequency when the signal varies slowly. Time resolved CE-MRA often has such smooth temporal variation so that FT along temporal direction works well as a sparsifying transform. To see the sparsity in Fourier spectrum, x-f image and f-axis plot at one x-point are illustrated with original x-t image and temporal variation at the same point in Figure 1. The majority of energy is concentrated at a few low frequencies since the signals smoothly vary.

Basically, k-t FOCUSS is based on compressed sensing and the sparsifying transform ψ does not need to be constrained to FT. If there are other bases to represent images more sparsely, those bases will be more effective for k-t FOCUSS reconstruction. Based on this idea, we already used principal component analysis (PCA) to find the most sparsifying transform for functional MRI (11). The resultant transform is often called Karhunen-Loeve transform (KLT) (14). KLT is well-known as an optimal energy compaction transform for signals (14). Hence, for non-periodic signal such as time resolved CE-MRA, KLT can be a more suitable sparsifying transform than FT.

Unlike the FT, KLT is a data dependent transform. More specifically, let σ_k denote a random vector spanning time series from t₁ to t_N.

[6]

σ_k = [σ(k,t₁), σ(k,t₂),…, σ(k,t_N)]^T

In general, the covariance matrix C of σ_k can be expanded as follows:

[7]

C = E [σ_{k} σ_{k}^{H}] = \sum_{l = 1}^{N} γ_{l} Ψ_{l} Ψ_{l}^{H},

where

{\{γ_{l}\}}_{l = 1}^{N}

and

{\{Ψ_{l}\}}_{l = 1}^{N}

are the eigenvalues and the corresponding orthonormal eigenvectors of C (14). Then, the original signal σ_k can be represented as follows (14):

[8]

σ_{k} = \sum_{l = 1}^{N} ρ_{l} Ψ_{l},

where

[9]

ρ_l = < σ_k, Ψ_l >.

The eigenvectors

{\{Ψ_{l}\}}_{l = 1}^{N}

and coefficients

{\{ρ_{l}\}}_{l = 1}^{N}

correspond to the sparsifying transform Ψ and unknown sparse signal ρ in Eq. [1], respectively.

The KLT provides a separation of the randomness and the time variation in the signal σ_k. In particular, the randomness is summarized in the sequence

{\{ρ_{l}\}}_{l = 1}^{N}

while the time variation in the process is embodied in the deterministic functions

{\{Ψ_{l}\}}_{l = 1}^{N}

(14). Therefore, if each pixel of CE-MRA has similar behavior along temporal direction, we can represent the temporally varying images very sparsely with a small amount of randomness.

However, estimating KLT bases is a very underdetermined problem when the data are highly accelerated. Dealing with this problem, we approximate KLT bases using only a partial data set. For example, in the Cartesian trajectory, as shown in Figure 2 (a), the fully sampled low frequency k-t data can be used to estimate the covariance matrix C. More specifically, let V be a fully sampled k-t data set in low frequency region:

[10]

V = (\begin{matrix} \begin{matrix} v (k_{1}, t_{1}) & v (k_{2}, t_{1}) \\ v (k_{1}, t_{2}) & v (k_{2}, t_{2}) \end{matrix} & \begin{matrix} \dots & v (k_{M}, t_{1}) \\ \dots & v (k_{M}, t_{2}) \end{matrix} \\ \begin{matrix} ⋮ & ⋮ \\ v (k_{1}, t_{N}) & v (k_{2}, t_{N}) \end{matrix} & \begin{matrix} ⋱ & ⋮ \\ \dots & v (k_{M}, t_{N}) \end{matrix} \end{matrix}),

where M and N are the total number of fully sampled low frequency k-data and the total number of time frames, respectively. Then, the covariance matrix C is obtained by:

[11]

C = VV^H.

Finally, the KLT bases

{\{Ψ_{l}\}}_{l = 1}^{N}

can be estimated using eigenvector decomposition of C as in Eq. [7]. Here, it is important to note that even though KLT bases are obtained from low spatial resolution images, the temporal changes are not smoothed at all because the resultant basis is complete in C^N and we can represent any temporal variation in C^N with linear combination of KLT bases. Therefore, the only concern about the approximated KLT bases is if the unknown signal can be represented with sparse coefficients ρ. We confirmed that the resultant coefficients are sparse enough in the following experiments. This is one of the main differences from partially separable function (PSF) approach (15), where only l significant principle components are used for temporal transform by assuming that the rank of the unknown x-t images is l. However, in real cases, it is hard to know the rank of the unknown x-t images a priori. In contrast, in our approach, the bases are complete and the sparse non-zero coefficients are obtained using l minimization.

In radial trajectory, since we do not have fully sampled low frequency measurements on the same coordinates along time, a novel method should be considered. In this paper, initial dynamic images are reconstructed with k-t FOCUSS using FT. Then, we set the mask M that may correspond to the vessel area. The mask M is determined by a simple thresholding:

[12]

M (x) = \{\begin{matrix} 1, if |σ̂ (x)| \geq ϒ \\ 0, otherwise, \end{matrix}

where σ̂ and Υ represent temporal average of initial k-t FOCUSS results and specified thresholding value, respectively. Temporally varying signals on the non-zero mask are then extracted to estimate the temporal covariance matrix. Then, KLT bases are obtained by eigen-decomposition as described above. Figure 1 briefly summarizes these steps and demonstrates the sparse signals after KLT. Through this process, the temporal variation due to streaking artifact patterns can be excluded in determining KLT bases.

Incoherent Sampling

In (11, 12), we achieved successful performance of k-t FOCUSS in Cartesian trajectory in cardiac cine imaging which is usually less sparse than MRA. Meanwhile, several high performance algorithms, such as highly-constrained back-projection (HYPR) (16) or iterative HYPR (I-HYPR) (17), use radial trajectory for MRA. Therefore, one may wonder which trajectory is more suitable for MRA. Since both trajectories have incoherent sampling basis with respect to image space, compressed sensing framework works for both cases. Figure 2 (a) and (b) show the artifact patterns in Cartesian and radial trajectories, respectively. The artifacts are spread over whole image space like noise pattern and this proves the incoherency of sampling basis in both trajectories.

However, due to differences in Cartesian and radial trajectories, there exists pros and cons for each method. For example, in radial trajectory, spatial resolution is mainly determined by imaging field of view (FOV) and there is smaller trade off between spatial resolution and number of views in contrast to the Cartesian trajectory. Furthermore, due to the oversampled k-space region, radial trajectory is more robust to motion artifacts (18) and more accurate contrast reconstruction is possible. Moreover, using golden ratio trajectories (13), we can distribute sampling trajectories without overlap on different time points. More specifically, by specifying the angle increment between consecutive views with 111.25°, the sampling incoherency can be preserved for arbitrary choices of the number of views and time frames for already acquired data set. From a computational complexity point of view, Cartesian trajectory is, however, more advantageous than radial trajectory. Specifically, the reconstruction can be done using fast Fourier transform (FFT), whereas radial trajectory requires computationally expensive projection/back-projection or gridding. Another issue comes from using KLT. As already explained in the perviovs section, radial trajectory require initial reconstruction to obtain KLT bases, whereas Cartesian trajectory does not need the initial reconstruction since fully sampled low frequency k-t data can be used for KLT basis estimation.

Up-Sampling

Increasing imaging FOV makes the resultant image sparser. According to Fourier slice theorem, each read-out data in radial trajectory corresponds to the Fourier transform of a sinogram obtained from projection of an image. Therefore, inverse Fourier transform of one read-out corresponds to the sinogram of original image. Hence, just zero-padding at both ends of each sinogram with arbitrary size, imaging FOV can be extended without any distortion of original image contents. The zero-padding scheme is a quite common technique in signogram-based reconstruction (19). It is usually employed for correct implementation of discretized filtered back projection. However, in our approach, we employed this zero-padding scheme for sparser representation of the unknown images.

Actually, when the number of views is close to Nyquist sampling criteria, this step is not necessary to avoid aliasing artifacts. However, if the number of views is very limited so that the empty k-space is estimated with small number of measurements, estimation error of k-space data can result in aliasing artifact in image domain especially when the size of FOV is similar to the size of image contents. In this case, employing up-sampling scheme mentioned above, image quality can be greatly improved without serious aliasing artifacts, as will be demonstrated later.

Materials

The accelerated simulated data were provided at the "Recon Challenge" of 2009 ISMRM workshop on Data Sampling and Image Reconstruction (http://www.ismrm.org/workshops/Data_09/recon.htm). Data were simulated from X-ray cerebral angiography. X-ray data were collected 3 frames per second, with 2048 in-plane resolution, for a total of 8 seconds (24 collected frames) which span wash-in to wash out of injected bolus. The images were blurred to 512 resolution, and linearly interpolated in time between frames to create a total of 512 time frames. Coil sensitivity maps obtained from an axial slice through a water phantom using an 8 channel head coil were superimposed on the image to create 8 coil images. Independent noise was added to each channel. From each temporal image, one read-out was obtained with size of 512 in radial trajectory so that the given total number of simulated k-data was 512×512. Here, the time interval between each interpolated time frame corresponds to TR in real MR acquisition. Then, reconstruction was performed at 12 time points corresponding to the 39, 79, 118, 157, 197, 236, 275, 314, 354, 393, 432, 472 nd time frames.

In radial acquisition, the view order was based on the golden ratio whose angle increment between consecutive views is 111.25°(13). In this scheme, there was no overlap of trajectories in 512 views and we could choose arbitrary number of views to reconstruct each time frame since any consecutive views are broadly distributed over k-space as described in (13). In this paper, the experiments were conducted at 12- and 51.2-fold acceleration. At 12-fold acceleration, all of the given 512 views were used for reconstruction. Specifically, k-space data at different time were aggregated with 39 or 40 sized windows whose centers were placed at target time frames, such as the 39, 79, 118, 157, 197, 236, 275, 314, 354, 393, 432, 472 nd. The window size 39 or 40 was determined by a time distance between current target frame and previous target frame. For 51.2-fold acceleration, we simply used 10 views for each time frame. From these data sets, the reconstruction was performed for each coil using k-t FOCUSS with various compressed sensing parameters. For up-sampling, we generated 1024 sized read-out by padding zeros at both ends of each 512 sized sinogram and the up-sampling factor was 2. For comparison of the performance in different trajectories, we additionally acquired randomly sampled data in Cartesian trajectory. Here, the acceleration factor was 12. Final results were obtained by root sum of squares of 8-coil images. For the quantitative evaluation of each result, the interpolated X-ray images corresponding to 39, 79, 118, 157, 197, 236, 275, 314, 354, 393, 432, 472 nd time frames are multiplied by coil sensitivity maps, and the root-square sum of squared images of 8-coils are used as a ground truth. Our winning result at the Recon Challenge was the result reconstructed using up-sampling and KLT at 12-fold acceleration factor.

Experimental Results

For the effect of different parameters, we illustrated radial k-t FOCUSS results using FT with no up-sampling, FT with up-sampling, and KLT with up-sampling together with ground truth in Figure 3. For fair comparison, the pixel values were equally scaled with fixed minimum and maximum values for different methods. Here, the acceleration factor was 12. When there was no up-sampling, the overall image quality was poor with streaking and aliasing artifacts as indicated by arrows. In contrast, after up-sampling, FT and KLT results illustrate no such artifacts. The magnified images help to see the detailed structures and contrast changes around vessels. It is observed that KLT results have the finest vessel structures.

For quantitative comparison of FT and KLT methods, we calculated mean square errors (MSE) at 51.2 and 12-fold acceleration factors in Figure 4. When acceleration factor is low, both methods show comparable MSE. However, when acceleration factor is high, the improvement in KLT result is significant. For further analysis about temporal resolution, temporal variations of averaged pixel values within specified 5×5 window were plotted together with true values and corresponding errors in Figure 5 for 12 and 51.2-fold accelerations. Similar to MSE results, when the acceleration factor is low, both results correctly follow the true values, whereas at high acceleration KLT result shows better performance.

In order to compare the performance in Cartesian and radial trajectories, we have conducted comparative studies in Figure 6. Here, the acceleration factor was 12 and KLT was used as a sparsifying transform for both trajectories. Figure 6 (a) shows one time frame. The images of both trajectories look nice without aliasing or streaking artifacts but in magnified version, it is observed that radial trajectory has finer vessel structures. Figure 6 (b) and (c) again confirm that radial results are better than Cartesian results by plotting spatial and temporal variation.

Finally, Table 1 briefly summarizes k-t FOCUSS using FT and KLT in Cartesian and radial trajectory in terms of reconstruction time and MSE. Reconstruction time was measured for each coil by CPU time on Xeon 3 GHz Linux platform with 4 GB RAM using Matlab 7.0.4. Then, MSE was calculated for each time frame. Here, the acceleration factor was 12 and this sampling ratio was enough to directly estimate KLT bases in Cartesian trajectory. Hence, the reconstruction time was similar in using FT and KLT in Cartesian trajectory. Similar to the radial trajectory, KLT result with Cartesian sampling shows smaller MSE than FT result. Next, temporal FT without up-sampling in radial trajectory shows very poor performance. By using up-sampling by factor of 2, the MSE value was significantly reduced. If KLT is used for temporal transform, the results can be further improved. Note that the problem of up-sampling is longer time for reconstruction. This is due to large memory requirement. Furthermore, since KLT basis should be estimated after initial reconstruction, KLT in radial trajectory requires additional computational time.

Conclusion

This paper described the k-t FOCUSS algorithm that won the Recon Challenge at 2009 ISMRM Data Sampling and Image Reconstruction workshop. Since MRA is usually composed of very sparse vessel structures, it is a very suitable application of compressed sensing. In this paper, we proved that k-t FOCUSS that has been successfully applied to cardiac cine imaging is also effective in time resolved CE-MRA. Several sparsifying transforms and various options have been compared. Through various experiments, we confirmed that KLT with radial trajectory with up-sampling outperforms other combinations of sparsifying transforms, sampling patterns, etc.