Improving Non-Cartesian MRI Reconstruction through Discontinuity Subtraction

Non-Cartesian sampling is widely used for fast magnetic resonance imaging (MRI). Accurate and fast image reconstruction from non-Cartesian k-space data becomes a challenge and gains a lot of attention. Images provided by conventional direct reconstruction methods usually bear ringing, streaking, and other leakage artifacts caused by discontinuous structures. In this paper, we tackle these problems by analyzing the principal point spread function (PSF) of non-Cartesian reconstruction and propose a leakage reduction reconstruction scheme based on discontinuity subtraction. Data fidelity in k-space is enforced during each iteration. Multidimensional nonuniform fast Fourier transform (NUFFT) algorithms are utilized to simulate the k-space samples as well as to reconstruct images. The proposed method is compared to the direct reconstruction method on computer-simulated phantoms and physical scans. Non-Cartesian sampling trajectories including 2D spiral, 2D and 3D radial trajectories are studied. The proposed method is found useful on reducing artifacts due to high image discontinuities. It also improves the quality of images reconstructed from undersampled data.


INTRODUCTION
Non-Cartesian data sampling, such as radial and spiral encoding, is often performed for fast data acquisition in MRI to capture the rapidly decaying signals. Therefore, non-Cartesian trajectories have gained increased attention in various applications such as functional brain imaging, hyperpolarized gas imaging, contrast-enhanced MR angiography, and cardiac imaging. However, their disadvantages lie primarily in more complicated image reconstruction. The most widely used MRI reconstruction scheme generally consists of three steps: preweighting, gridding by convolving the original data with interpolation kernels, and d-dimensional inverse Fourier transform [1]. Conventional gridding methods have used predefined gridding kernel functions, such as triangle, Gaussian, sinc, and Kaiser-Bessel functions, and so forth [1,2].
The recent introduction and development of NUFFT algorithms [3][4][5][6][7] provide an efficient tool to evaluate the nonuniform discrete Fourier transform (NUDFT). The NUFFTs calculate optimal kernels with selected scaling function, so that the interpolation error is minimized in a least-square sense. Some of these NUFFT algorithms have been applied to spiral and radial MRI reconstructions [8,9] and achieved better accuracy-speed trade off. These NUFFT-based direct reconstructions are referred to as direct reconstruction methods in the rest of the paper.
Unfortunately, NUFFT reconstruction, under the same gridding scheme, provides improvement only on the interpolation accuracies. In fact, image quality is still challenged by k-space truncation and density compensation for the non-Cartesian sampling, and as a result, the point spread function (PSF) is not as sharp. For an imaging system, PSF is one of the most important characteristics and ideally should be a single delta function. Lauzon and Rutt [10] studied the "principal PSF," the Fourier transform of the sampling function, for polar (2D radial) sampling trajectory and observed intensity leakage (existence of sidelobes and broadening of the mainlobe). Since the final image is a convolution of the PSF with the true magnetization distribution, the leakage of PSF then results in intensity leakage in the reconstructed image. Furthermore, the intensity of artifacts is proportional to the intensity of the region being convolved. As a result, the artifacts leaked by high intensity regions into their surrounding 2 International Journal of Biomedical Imaging low intensity regions become most visible. For example, the skull in brain images usually has very high intensity and its leakage into inside of the brain significantly degrades the image quality in various regions of interest.
A straightforward approach to improve the reconstruction result is to sharpen the PSF, which directly corresponds to an optimized selection of density compensation function (DCF). For non-Cartesian sampling, DCF is applied to the raw data to compensate for the sampling density. However, the mathematically ideal DCF [1,[11][12][13] may not provide highest image quality. Some recent works [14][15][16][17] formulate the DCF computation as an optimization problem to emphasize certain PSF characteristics, for example, narrow mainlobe and low sidelobes. Nevertheless, the optimized solutions are rarely close to being ideal and the computational cost of optimization becomes inevitably high when large data sets are processed.
Instead of correcting PSF itself, this paper develops an alternative method with discontinuity subtraction, based on the fact that intensity leakage is most severe at high discontinuities. Here we define "high discontinuity" as the structural region that has high intensity difference (contrast) with the surrounding regions, such as the skull. Similar ideas can be found in some other applications such as computerized tomography [18] and seismic data processing [19]. However, the proposed method does not require any a priori knowledge about the image, which is needed in [18]. The computational cost is slightly higher than direct reconstruction but much lower than iterative methods.
The paper is organized as follows. Section 2 develops multidimensional NUFFTs that are used for reconstruction and k-space simulation. These NUFFTs (NUFFT-1 and NUFFT-2) are extensions of Liu and Nguyen's idea of 1D least-square NUFFT [5,6]. In Section 3, the principal PSF of the direct NUFFT reconstruction method [9] is analyzed and the scheme of the discontinuity subtraction method is described. The proposed method is then compared with the direct NUFFT reconstruction method on both spiral and radial trajectories through computer simulation and physical scan studies in Section 4, followed by discussions in Section 5.

k-space simulation and reconstruction
In MRI, the spatially encoded MRI signal s(k) in k-space can be modeled as the spatial-frequency response of the distribution of the initial transverse magnetization M xy in spatial domain x. The acquired raw data are therefore Fourier coefficients of the underlying M xy whose image I x can be reconstructed by taking inverse Fourier transform of the sampled signal in k-space. Generally, in d-dimension, where t is the readout time.
For numerical evaluation, the integrands must be truncated and discretized. Note that in (1), M xy (x) usually has a finite support because the object is physically bounded, and thus no truncation is needed. Equation (1) does not need a density compensation either since the image is uniformly sampled. However, (2) would require both truncation and density compensation because of the infinite k-space support and nonuniform sampling. Discretizing (1) in Cartesian spatial domain and (2) in non-Cartesian k-space gives which assumes the image voxel size in any pth dimension Δx p = 1, and where w(k) is the preweighting factor (or DCF). In this work, we use Voronoi DCF [13] for spiral trajectory and the Jacobian DCF for radial trajectory. Equation (3) is named as "k-space simulation formula" that simulates non-Cartesian k-space samples from a known image, and (4) as "reconstruction formula" that recovers the image from known non-Cartesian k-space samples. Since (4) only involves the inverse Fourier transform once, we refer to this reconstruction as direct reconstruction.
In this work, (3) and (4) are evaluated by the multidimensional NUFFTs, which will be briefly summarized in Section 2.2. Comparison of these NUFFT reconstruction and NUFFT k-space signal simulation to conventional gridding methods can be found in detail in our previous work [9]. Again, one should note that (3) provides high-fidelity estimation, while (4) is limited in accuracy and usually consists of leakage artifact, which will be discussed in later sections.

NUFFT-1 and NUFFT-2
Multidimensional NUFFT-1 (NUFFT of the first type) and NUFFT-2 (NUFFT of the second type) [9] are developed for fast evaluation of the NUDFTs between the non-Cartesian kspace and the Cartesian image space. NUFFT-1 evaluates with the following procedures: (1) calculate a set of ddimensional kernels Ψ m with least-square criterion for each of the non-Cartesian samples G m at v m ; (2) take the ddimensional regular FFT on g n ; (3) convolve the FFT results with Ψ m and find G m . Take NUFFT-1 as an example. Each kernel Ψ m is computed as a tensor product of 1D least-square kernels [5,6] where μ is the oversampling factor and q is the kernel size.
The kernel varies from sample to sample, but it is only dependent on the location of each sample. NUFFT-1 is utilized to evaluate (4) and NUFFT-2 is used for (3). Note when NUFFT-2 is used to evaluate (3), the kernel is in fact the conjugate of Ψ m due to the negative sign in the exponent. As a result, when NUFFT-1 and NUFFT-2 are used together, the two kernels only need to be calculated once and stored to be called at every NUFFT operation.
Our numerical studies find that the NUFFTs compute (6) and (5) at a computational cost which asymptotically approaches that of a regular FFT with L 2 approximation error of 10 −4 . In this work, we use μ = 2 and q = 4.

The point spread function and intensity leakage
System impulse response, also known as point spread function (PSF), is one of the most important performance measurements of an imaging system. The resulting image can be interpreted as a convolution of the object and the system PSF. If we assume the system is ideal before image reconstruction, that is, the raw data are exact spatial frequency responses at the sampling locations, the system PSF becomes simply the Fourier transform of the sampling function (9) where the sampling function and W(k) is the DCF. We define (9) as principal PSF, to be consistent with the one of the polar (radial) MR imaging system given by [10]. It is shown that the PSF is the DFT of DCF. The desired PSF function should have narrow mainlobe width and low sidelobe level.
In Figure 1, we show PSFs of radial and spiral sampling trajectories. The radial trajectory has 400 rays and each ray is oversampled with 183 points within a k-region of [−0.5, 0.5].
The spiral trajectory has 16 interleaves and 2048 samples on each interleave. The PSFs are calculated on a 128×128 Cartesian grid. Both PSFs, as shown in Figure 1, have significant sidelobes and the mainlobe is broadened. When convolved with the ideal object, these PSFs introduce significant intensity leakages. The caused ringing artifacts are most visible at high intensity discontinuities.
In practice, undersampled radial or spiral acquisitions are widely used by MRI researchers to capture rapid physiological images. In these functional studies, temporal resolution usually trades off spatial resolution, so that undersampling is often used to improve the temporal resolution. For non-Cartesian sampling, wellsampling refers to a sampling scenario whose lowest sampling density still satisfies Nyquist criterion, and any scenarios with lower sampling density than that are considered as undersampling. Figure 2 shows that undersampling introduces replica sidelobes besides the normal leakage. These PSF imperfections will be translated to ringing artifacts and streaking artifacts in the images.

Reconstruction and leakage reduction
To reduce these artifacts, one approach is to optimize DCF for the desired properties of PSF, that is, narrow mainlobe and low sidelobes. However, this approach is quite complicated and minimizing the sidelobes is usually at the cost of widening the mainlobe. Alternatively, we consider the problem from a new perspective and propose a discontinuity subtraction method to reduce the artifacts caused by the high discontinuities. Furthermore, unlike the subtraction method suggested by [18], we use direct reconstruction result to find the discontinuities so that no a priori information about the image is needed. During each iteration thereafter, everything is transformed back to k-space as a projection onto convex sets (POCS). Subtraction is applied in k-space and the data fidelity is ensured. The detailed procedures of our NUFFTbased leakage reduction method are the following.
(1) Given the non-Cartesian k-space data s 0 (k), apply (4) to find the estimated image I 0 (x). Here, gradient edge detection and thresholding segmentation techniques are used for simplicity. (2) Detect the highest discontinuity I d1 in the estimated image with existing edge detection and image segmentation techniques.  NUFFT-1 and NUFFT-2 needs to be calculated only once. For example, an image with m high discontinuities will need only one pre-processing, m + 1 NUFFT-1 and m NUFFT-2. It takes less than 2m + 1 times as long as direct reconstruction. Given the number of high discontinuities m in a practical medical image is usually small, the method is faster than most iterative reconstruction methods. For example, m is typically 2-3, which corresponds to 2-3 iterations and takes 5-7 times as long as direct reconstruction.

Computer simulation study
The Shepp-Logan phantom is utilized here for computer simulation study. It is known as the phantom of "head," which is usually used to simulate brain image with the skull. The non-Cartesian k-space signal samples on both radial and spiral trajectories are accurately simulated using (3). The Acquired data in k-space s 0 (k) Enforce data fidelity in k-space s 0 (k) = s 0 (k) s di (k) (4) Direct reconstruction into image space Still discontinuity in image space?

Yes No
Final image in image space Extract discontinuity in image space I di (x) Figure 3: Flowchart of the discontinuity subtraction reconstruction algorithm.
proposed leakage reduction method is compared with the direct reconstruction on both 2D and 3D phantoms. Figure 4 shows reconstruction results of a 128 × 128 phantom simulated with a 16-interleaved spiral trajectory. Two discontinuities are subtracted during the reconstruction. For illustration purposes, the first high discontinuity image, the corresponding subtraction remainder, and the smoothed remainder are shown in the figure. The direct reconstruction in Figure 4(a) suffers from ringing artifacts caused by the high discontinuity at the "skull." Figure 4 Figure 4(i). The oscillations around the "skull" are visibly mitigated by discontinuity subtraction. Reconstructions of the same 128 × 128 phantom simulated with 400 radial rays are shown in Figure 5. Again, the intensity leakages are well improved by discontinuity subtraction. Table 1 compares the performance of the two reconstruction methods quantitatively and lists the L 2 reconstruction errors defined as e 2 = I recon − I 2 / I 2 , where I recon is the reconstructed image and I is the original phantom. Both errors decay with increased spatial resolution, but the proposed leakage reduction method performs consistently better than the direct reconstruction method.
In Section 3.1, we discussed the influences of undersampling on the principal PSF. Reconstructing these undersampled data sets becomes more challenging. In Figure 6, we simulated the k-space samples of a 128 × 128 phantom on only 120 rays (70% undersampled) and applied both direct and leakage reduction methods to reconstruct the image. The streaking artifacts are much less significant when discontinuity subtraction (Figure 6(b)) is applied. The three tiny  ellipses near the bottom of the phantom are not resolvable in Figure 6(a) but can be clearly differentiated in Figure 6(b). This is better shown by the line profiles across those three ellipses in Figure 6(c). This suggests an improvement of resolution. In fact, the L 2 reconstruction errors are, respectively, 15.91% and 4.33% for Figures 6(a) and 6(b). A 64 3 Shepp-Logan phantom is also simulated using 3D isotropic radial trajectory to study the performance of the proposed method in 3D. 12868 rays are simulated to make it a well-sampled case just satisfying the Nyquist sampling rate. Figure 7 shows the three cross planes of the original phantom and reconstructed images using direct reconstruction and leakage reduction reconstruction. The ringing artifacts around the "skull" are significant in Figure 7(b) but are greatly suppressed in Figure 7(c) where discontinuity subtraction is applied. The L 2 errors are 25.36% for Figure 7(b) and 11.39% for Figure 7(c), respectively. Table 2 summarizes the reconstruction processing time of the above experiments using direct method (T 1 ) and the proposed method (T 2 ). The 2D experiments in the table used image size 128 2 and the 3D case used 64 3 . All of the reconstruction methods are implemented with MATLAB 7.0 on the same computer (CPU: AMD ATHLON 2600+; RAM: 2MB).

Physical scan study
Physical MR phantom scans were conducted using spiral and radial samplings as well. A phantom with sharp edges was scanned with a 6-interleaved spiral trajectory. Figure 8 shows the reconstructed images using the two reconstruction methods and center line profiles. Ringing artifacts were generated by the high discontinuities in the phantom structure. Applying the leakage reduction method partially suppressed the artifacts and improved image quality.
Radial scans were also conducted to study the performance of the two methods on undersampled data sets. A cylindrical water tube phantom was scanned by running a 2D radially encoded spin echo sequences on a 2-Tesla Oxford horizontal magnet. We first collected an 800-ray data set to support the Nyquist sampling rate for a 256 × 256 image, then reduced the number of rays to 200 and collected another data set. Direct reconstruction was applied to reconstruct both data sets. It is found that the image of the wellsampled case (Figure 9(a)) looks decent but the one of 3/4 undersampled case in Figure 9(b) has evident-streaking artifacts and low signal-to-noise ratio. These artifacts are mitigated in Figure 9(c) where the proposed leakage reduction reconstruction was applied, and the image quality is visually close to the well-sampled case.

DISCUSSIONS
The studies described show that the proposed leakage reduction method is helpful on reducing the ringing artifacts near the edges and the streaking artifacts due to undersampling. In this section, we further the discussion on the challenges left to the current algorithm, its applicability, and some potential research directions for improvements.  We found that discontinuity detection, via image segmentation, is quite important to the algorithm. Although the method does not require any a priori knowledge of the object, it should be noted that it does rely on detecting the discontinuity from the direct reconstruction using postprocessing techniques. The minimum requirement is that the discontinuities in the direct reconstruction result should be well detectable. It is found that the algorithm may not improve over the direct reconstruction when segmentation does not perform well (such as cases with low discontinuities). We find some cases, for example, those with low signal-to-noise ratio or significantly undersampled k-space, extremely challenging. However, it should also be pointed out that the algorithm merely requires a rough segmentation of the current highest discontinuities, and this is usually not as challenging as image segmentation tasks in a normal sense. Second, due to the POCS process that ensures data fidelity in the k-space, the algorithm at least should not perform worse than direct reconstruction. Given the above challenges, the current algorithm is most applicable to imaging objects that involve high intensity discontinuities, for example, the imaging of brain with skull and cardiac imaging. In these applications, the lower intensity anatomies are of more interest so that the reduction of leakages from the higher anatomies into the low intensity regions becomes very useful.
Based on the analysis, we suggest several research directions. First of all, although segmentation is not the emphasis in the work, there is a large body published work on this topic. We refer the reader to this literature for more effective segmentation methods, and we believe that integrating them in the current framework will result in even more improved results. Second, the proposed method is a good alternative to trade off image quality and computational complexity, but it does not guarantee optimal suppression of artifacts. Some recent works [20,21] suggested encouraging results using constrained optimization methods that minimize certain edge preserving priors (e.g., total variation or wavelet coefficients) with the constraint of data fidelity. The concept can be integrated into our scheme as well, for example, by updating the image using the gradient of the prior(s) before applying the discontinuity subtraction.

CONCLUSIONS
A leakage reduction reconstruction method has been proposed for improving non-Cartesian MR image reconstruction. The study explored two major non-Cartesian sampling trajectories, spiral and radial. Studies have shown that principal PSF imperfections cause significant image artifacts when discontinuous structures exist or when the k-space signal is undersampled. The proposed method uses discontinuity subtraction and the discontinuity is found in the image space by post-processing techniques. The subtraction is taken in the k-space, in order to enforce data fidelity. In this scheme, the multidimensional NUFFTs are exploited for k-space simulation and direct image reconstruction. Both computer-simulated Shepp-Logan phantoms and physical scanned data were reconstructed to compare the proposed method with the direct reconstruction method. Results suggested that the proposed method successfully reduced artifacts in images with discontinuities or ones reconstructed from incomplete data. Qing Huo Liu is a Professor of electrical and computer engineering at Duke University. He received the Ph.D. degree in electrical engineering from the University of Illinois at Urbana-Champaign, in 1989. His research interests include computational electromagnetics and acoustics, biomedical imaging, inverse problems, geophysical subsurface sensing, electronic packaging, and the simulation of photonic and nanodevices. He is a Fellow of the IEEE and a Fellow of the Acoustical Society of America. Currently, he serves as an Associate Editor for Radio Science, and an Associate Editor for IEEE Transactions on Geoscience and Remote Sensing. He received the 1996 Presidential Early Career Award for Scientists and Engineers (PECASE) from the White House, the 1996 Early Career Research Award from the Environmental Protection Agency, and the 1997 CAREER Award from the National Science Foundation.