Synopsis

Magnetic Resonance (MR) images are obtained from the measured k-space Fourier coefficients. Due to acquisition time limitations and noise typically only a limited part of the k-space is acquired, resulting in Gibbs ringing on the images. We propose an efficient and accurate local reconstruction method that removes the Gibbs ringing and yields sharp image profiles near local edges.

Introduction

Magnetic Resonance Imaging (MRI) inherently suffers from Gibbs artifacts because only a limited fraction of the k-space can be acquired. Conventionally, Gibbs ringing has been addressed with image filters such as Lanczos, which are associated with a loss of spatial resolution. More advanced techniques showed varying efficacy for strong and soft edges¹² and were computationally expensive because they were carried out in a global manner.^{4, 5, 15}

In this work, we propose a reconstruction method that focuses on only a local subdomain of the image at a time and, thus, yields an accurate and non-oscillatory sharp image. Patching all subdomains together yields the whole image.

Methods

We extended our 1D Domain Decomposition Fourier reconstruction method¹⁴ to 2D Fourier image reconstruction.

Two approaches were compared:

• The global 2D Fourier continuation sparse polynomial annihilation (PA) method and
• the dimension-by-dimension Fourier continuation sparse PA method.

The global 2D Fourier continuation sparse PA method first divided the 2D image into multiple sub-images and, then, applied the global 2D Fourier continuation to the extended sub-images to find the 2D periodic extension of the extended sub-images we are interested in. After finding the new Fourier coefficients based on the periodic extension, we applied the 2D sparse PA regularization to obtain the reconstruction.

For the dimension-by-dimension Fourier continuation sparse PA method, we first divided the 2D image into multiple sub-images and, then, applied the 1D Fourier continuation in y- and x-directions separately to find a periodic extension of the 1D data. Then, we found new Fourier coefficients based on the periodic extension and applied the 1D sparse PA method on these coefficients, which yielded the reconstruction.

We applied the technique in a Shepp-Logan phantom and in rat brain images acquired at 9.4 Tesla.

Results

Figure 1 illustrates sharper reconstructions by the dimension-by-dimension method near structural edges than by the global method. Root mean square error (RMSE) values were smaller for the dimension-by-dimension method, indicating a more accurate reconstruction.

Figure 2 shows the reconstruction obtained by first dividing the whole image into 16 equal sub-images then applying the dimension-by-dimension method on each extended sub-image and finally stitching the 16 reconstructed sub-images together to obtain the reconstruction in the original whole domain. Both reconstructions were Gibbs free and there was no oscillation near the fusion sites. RMSE values for the reconstruction with 5 and 10 extra points were smaller than for the reconstruction with 2 extra points.

Figure 3 shows the reconstruction by the dimension-by-dimension method in the MRI data, removing the Gibbs oscillation and yielding sharper images than the reconstruction by the global method near the edges.

It took about 370 seconds to complete the reconstruction with the global method with 2 extra points and about 16 seconds to use parallel computing with four cores to complete the reconstruction by the dimension-by-dimension method with 2 extra points (20 seconds to complete it with 5 extra points), showing that the dimension-by-dimension method is computationally more efficient than the global method.

Discussion

The dimension-by-dimension method can be extended to three-dimensional problems by repeating the same procedure in the z-direction after applying the method in both x and y-directions.

Most reconstruction methods are carried out in a global manner, in this work, we proposed a reconstruction method that focuses on only a local subdomain of the image at a time and, thus it is computational efficiency. A limitation of our proposed technique is that the best λ in the sparse PA method is problem-dependent. the parameter needs to be determined once for a certain image contrast, for example, in our work we found the best λ for the rat brain image, then we could use the same λ for other rat brain image. For our future research, we will consider an efficient way to reduce the computational complexity of the proposed method. In this paper, we found the λ-values empirically and a more systematic procedure is needed. Future research will also investigate the stability of the global method further and try to devise an efficient 2D global method.

Conclusion

Acknowledgements

References

[1] CVX: Matlab software for disciplined convex programming. CVX Research Inc. http://cvxr.com/cvx/; 2017.

[2] Archibald R, Gelb A, Platte RB. Image reconstruction from undersampled Fourier data using the polynomial annihilation transform. J. Sci. Comput. 2016;67(2): 432–452.

[3] Archibald R, Gelb A, Yoon J. Polynomial fitting for edge detection in irregularly sampled signals and images. SIAM J. Numer. Anal. 2005;43(1):259–279.

[4] Abdul JJ. Advances in the Gibbs phenomenon. Sampling Publishing. 2011.

[5] Abdul JJ. The Gibbs phenomenon in Fourier analysis, splines and wavelet approximations. Springer US. 1998.

[6] Boyd JP. A comparison of numerical algorithms for Fourier extension of the first, second, and third kinds. J. Comput. Phys. 2002;178:118–160.

[7] Bruno OP, Lyon M. High-order unconditionally stable FC-AD solvers for general smooth domains I. Basic elements. J. Comput. Phys. 2010;229:2009–2033.

[8] Bruno OP, Han Y, Pohlman MM. Accurate, high-order representation of complex three-dimensional surfaces via Fourier continuation analysis. J. Comput. Phys. 2007;227:1094–1125.

[9] Bruno OP, Elling T, Sen A. A Fourier continuation method for the solution of elliptic eigenvalue problems in general domains. Mathematical Problems in Engineering. 2015;1–15.

[10] Bruno OP. Fast, high-order, high-frequency integral methods for computational acoustics and electromagnetics, in: Topics in Computational Wave Propagation. Lecture Notes in Computer Science Engineering, Springer, Berlin. 2003;31:43–82. October 31, 2018 15/26

[11] Hesthaven J, Gottlieb S, Gottlieb D. Spectral methods for time-dependent problems. Cambridge University Press, Cambridge; 2007.

[12] Li P, Gao Z, Don W-S, Xie S. Hybrid Fourier-continuation method and weighted essentially non-oscillatory finite difference scheme for hyperbolic conservation laws in a single-domain framework. J. Sci. Comput. 2015;64:670–695.

[13] Shahbazi K, Albin N, Bruno OP, Hesthaven JS. Multi-domain Fourier-continuation/WENO hybrid solver for conservation laws. J. Comput. Phys. 2011;230(24):8779–8796.

[14] Shi R, Jung J. A domain decomposition Fourier continuation method for enhanced L 1 regularization using the sparsity of edges in reconstructing Fourier data. J. Sci. Comput. 2018;74(2):851-871.

[15] Wasserman G, Archibald R, Gelb A. Image reconstruction from Fourier data using sparsity of edges. J. Sci. Comput. 2015;65(2):533–552. October 31, 2018