Two-dimensional local Fourier image reconstruction via domain decomposition Fourier continuation method

The MRI image is obtained in the spatial domain from the given Fourier coefficients in the frequency domain. It is costly to obtain the high resolution image because it requires higher frequency Fourier data while the lower frequency Fourier data is less costly and effective if the image is smooth. However, the Gibbs ringing, if existent, prevails with the lower frequency Fourier data. We propose an efficient and accurate local reconstruction method with the lower frequency Fourier data that yields sharp image profile near the local edge. The proposed method utilizes only the small number of image data in the local area. Thus the method is efficient. Furthermore the method is accurate because it minimizes the global effects on the reconstruction near the weak edges shown in many other global methods for which all the image data is used for the reconstruction. To utilize the Fourier method locally based on the local non-periodic data, the proposed method is based on the Fourier continuation method. This work is an extension of our previous 1D Fourier domain decomposition method to 2D Fourier data. The proposed method first divides the MRI image in the spatial domain into many subdomains and applies the Fourier continuation method for the smooth periodic extension of the subdomain of interest. Then the proposed method reconstructs the local image based on L2 minimization regularized by the L1 norm of edge sparsity to sharpen the image near edges. Our numerical results suggest that the proposed method should be utilized in dimension-by-dimension manner instead of in a global manner for both the quality of the reconstruction and computational efficiency. The numerical results show that the proposed method is effective when the local reconstruction is sought and that the solution is free of Gibbs oscillations.


Introduction
Magnetic resonance imaging (MRI) is a commonly used medical imaging technique. MRI image reconstruction is based on the inverse Fourier transform of a frequency-limited acquired Fourier spectrum of the object. In this work, we are interested in the 2D Fourier reconstruction of given MRI data. Before we continue to the 2D MRI image reconstruction, we should start from the basic 1D Fourier reconstruction. For the 1D Fourier reconstruction, we consider the problem of how to reconstruct a piecewise smooth function f ðxÞ : O ! R, O = [−1,1] on a uniform grid x j ¼ À 1 þ jDx; j ¼ 0; 1; � � � ; 2N; Dx ¼ 1 N when its Fourier coefficients (also known as k-space data in MRI), ff k g, are given. The Fourier coefficients of f(x) are defined byf and the Fourier partial sum f N (x), which is the image space reconstruction commonly obtained in MRI, is Since the function that we want to find, f(x), is a piecewise smooth function, the approximation by Eq (2) may yield the Gibbs phenomenon. When the patient's MRI image is obtained, usually we want to reconstruct the image so as to reveal the detailed structure of a particular region with the given Fourier coefficients. However, most reconstruction methods, such as the Fourier reconstruction and filtered Fourier reconstruction, are carried out in a global manner [1,2], so they are computationally expensive. And since global methods provide the reconstruction in one piece, so the quality of the reconstruction near different edges is not equal ( [3] shows this limitation). In this work, we propose a local method that focuses on and yields the local reconstruction of the subdomain such that an accurate and non-oscillatory sharp image reconstruction is achieved in that subdomain. If we patch all these local reconstructions together to constitute the whole image, the result will be more accurate than the one obtained by the global method while the subdomain that we are interested in is much enhanced.
As an example, we proposed a local method [3] based on the L 1 regularization method proposed in [4]. The method by [4], known as the sparse polynomial annihilation (PA) method, is an L 1 regularization method based on sparsity of edges for the Fourier reconstruction, for obtaining the sharp image profiles near edges. By using the edge sparsity in the L 1 regularization term, this approach can reduce the Gibbs oscillations near edges and provide a sharp reconstruction. However, as other reconstruction methods, this approach is a global method, which intrinsically puts more weight on strong edges (jumps with large magnitudes) than on weak edges (jumps with small magnitudes). This can sharpen the reconstruction near strong edges, but the reconstruction near weak edges often becomes similar to or even worse than the original Fourier reconstruction. Thus, if the area we are interested in is near weak edges, this approach may not be effective. If the L 1 regularization can be done locally, the reconstruction near weak edges can be improved. Such local operation could be done by our previous 1D method of domain decomposition Fourier L 2 and L 1 minimization based on the Fourier continuation [3]. The key element of our local method is to split the given domain into multiple subdomains first and then carry out the L 1 minimization individually in each subdomain. If necessary, we patch all the reconstruction in each subdomain together.
In this paper, we extend our previous 1D domain decomposition Fourier reconstruction method to 2D image reconstruction. We propose the following 2D methods, 1) the global 2D Fourier continuation sparse PA method and 2) the dimension-by-dimension Fourier continuation sparse PA method. The dimension-by-dimension is basically a 1D Fourier continuation applied in x− and y−directions separately. By comparing numerically the global 2D Fourier continuation sparse PA method with the dimension-by-dimension Fourier continuation sparse PA method, we show that the dimension-by-dimension method should be used for both accuracy and computational efficiency. In this paper, we also provide empirically the best range of the control value of the L 1 regularization, λ, to use. Our method is efficient and flexible in the sense that we can stitch all the local reconstructions from each subdomain to form the accurate whole complete image and that the constituted one is more accurate than the image by the global method. In this paper we provide detailed numerical results using the Shepp-Logan phantom image and real MRI data.
This paper is composed of the following sections. First, a brief explanation of the 1D domain decomposition Fourier continuation method based on the edge sparsity proposed in [3] is given. Then the extension of the 1D method to the 2D Fourier reconstruction will be introduced. Various numerical examples are given immediately after the 2D method. Finally, a concluding remark will be provided.
Global Fourier continuation method. Suppose that we have a non-periodic function f : The point values of f, f(x j ), are given at x j = a + jΔx, where Dx ¼ L N and j = 0, � � �, N. With the global 1D Fourier continuation method, the nonperiodic function f(x) is extended to a periodic function g(z) over I 2 = [a, b + d] (where d = γΔx, γ is a positive integer) with periodicity of L + d. The positive value of γ is arbitrary and the optimal value of γ is function dependent. Practically the value of γ is chosen empirically in consideration of computational efficiency.
The periodic extension g(z) can be written as the following Fourier sum since g(z) is periodic over I 2 where M is a nonnegative integer. The unknown coefficientsĝ k can be obtained by matching where z j = a + jΔx = x j , j = 0, 1, � � �, N. Let the element of the coefficient matrix A be defined by If M = N/2,ĝ k are uniquely determined. If M 6 ¼ N/2, the system can be solved by using the pseudo inverse based on the singular value decomposition (SVD), where A + denotes the pseudo inverse of A. The matrix A easily becomes ill-posed with the large values of N and M, which leads us to the method in the following section. Fourier continuation method using the boundary values. The global Fourier continuation method is simple, but since this method uses every f(x i ) to find the extended periodic function, the reconstruction is easily affected by round-off errors and the Runge phenomenon still exists if f(x) is a Runge function. To minimize such issues, the local Fourier continuation method [10,11] has been used which utilizes f(x i ) near the domain boundaries only. We briefly explain this local Fourier continuation method below.
The extended periodic function g(z) over [a, b + d] is defined as where f match (z) is called the matching function. Here, g(z) is defined different from g(z) in Eq (3). The local Fourier continuation method uses only local values near the boundaries to find the matching function. Define the intervals With the values of f(x) over I left and I right , the matching function f match (z) over these intervals is defined by f match (z) is a periodic function defined over I 3 = [b − δ, b + 2d + δ] with the periodicity of 2(d + δ). By finding f match , the extended periodic function g(z) in Eq (6) is found. To find f match (z) over the whole interval I 3 , we use following formula wheref m k are the unknown Fourier coefficients of f match (z). We findf m k by matching the matching function with the given function. That is, and at Q distinct points in each interval. K is an integer that K < Q. Here p(z) and q(z) are found by using (β + 1) grid points such that the degree of p(z) and q(z) is less than or equal to β. If the degree is β, then p(z) and q(z) are the local interpolation of f(x) over I left and I right , respectively. If the degree is less than β, then p(z) and q(z) can be found by solving the over-determined system in the least-squares sense. To minimize round-off errors, people usually use the Gram polynomial for p(z) and q(z) and construct the linear system [10]. By solving Eqs (9) and (10) in the least-squares sense using SVD we findf m k . Then using Eqs (8) and (6), we find the matching function f match (z j ), j = N, � � �, N + γ and obtain the desired extended function g(z j ), j = 1, � � �, N + γ.

Convex optimization: L 1 regularization based on sparsity of edges
The following convex optimization problem using L 1 regularization based on sparsity of edges was proposed in [4], which yields sharp Fourier reconstruction near edges, where k�k 2 and k�k 1 denote the vector L 2 and L 1 norms, respectively. The first term is the fidelity term. f R is the reconstruction we want to find andf is the given Fourier coefficient vector. The edge operator E p is the sparse polynomial annihilation transform and the superscript p denotes the order of the derivative of the interpolation [12,13]. The polynomial annihilation (PA) is basically higher order derivative of the interpolation. Thus E p f R has large values at the discontinuities but vanishes in the smooth regions of the function. In such a way, E p f R represents sparsity. The constant λ > 0 is a free parameter whose optimal value is chosen empirically [4]. In [4] it was shown that the L 1 regularization with the sparse PA method yields a better reconstruction than the filtering or TV regularization.

1D domain decomposition Fourier continuation sparse PA method
In [3] the domain decomposition sparse PA method was proposed, with which the given domain is split into multiple subdomains and the sparse PA method is applied separately in each individual domain. To make the minimization of Eq (11) done locally, we used the Fourier continuation method.

Remark 1
Here we should note that our domain decomposition Fourier continuation method is not limited to the sparse PA method. The L 1 regularization term in Eq (11) can be replaced with other terms such as TV norm.
For the domain decomposition Fourier continuation method, we split the domain into multiple subdomains f½x 0 ; x j 1 �; ½x j 1 ; x j 2 �; � � � ; ½x j kÀ 1 ; x j k �; ½x j k ; x 2N �g, 0 < j 1 < j 2 < � � � < j k < 2N and carry out the minimization for each subdomain separately. For example, suppose that we reconstruct the function on the subdomain ½x j 1 ; x j 2 �. The function is not necessarily periodic in this subdomain. To deal with this non-periodic function with the Fourier method, we use the Fourier continuation method. However, as shown in [3], the direct application of the Fourier continuation method to the subdomain ½x j 1 ; x j 2 � causes some oscillations near the domain boundaries when trying to stitch the all the reconstructions from each subdomain. To avoid such oscillations and apply the Fourier continuation successfully, we first extend the subdomain so that the extended domain contains ½x j 1 ; x j 2 �. The easiest way to find this extension is to add a fixed number of extra grid points to both boundaries. A more sophisticated way is presented in our previous work [3] based on the TV analysis. Let n m be a fixed integer that determines how many extra grid points are added at each boundary. Thus the extension becomes ½x j 1 À n m ; x j 2 þn m � :¼ ½a; b�, which contains K 0 = (j 2 + n m ) − (j 1 − n m ) + 1 uniform grid points. Thus we have the non-periodic data f N (x)(obtained from all Fourier coefficients via Eq 2) on K 0 uniform grid points over [a, b].
Then we use the following convex optimization method in this subdomain using L 1 regularization based on the sparsity of edges F is a transform matrix explained below.f c is a Fourier coefficient vector which needs to be found before solving Eq (12). Thisf c is not the one defined by Eq (1), but is obtained from For the extended subdomain, [a, b] we apply the Fourier continuation method to f N (x) over [a, b] [3] and obtain the extended periodic function γ is a fixed integer that determines the number of points on the extended interval as defined in the previous section.
, the corresponding Fourier coefficients,f c k , can be found by using the discrete Fourier transform [14], Þ T and F be the transform matrix whose elements are The optimization problem Eq (12) is solved using the Matlab CVX package [15]. f R is the reconstruction we get over the extended interval [a, b + d]. The reconstruction by Eq (12) is over the domain [a, b + d]. We only take the reconstruction over the interval

2D domain decomposition Fourier continuation method
Suppose that we have a periodic function f 0 : Further suppose that the values of f 0 are given at a set of uniform grid points, (x i , y j ), x i = a 0 + iΔ and y j = a 0 + jΔ, where D ¼ L 0 N and i = 0, � � �, N, j = 0, � � �, N. We consider the problem of how to reconstruct the 2D image on n × n uniform grid points over In practice, the subdomain size is determined by the size of the region where the local Gibbs oscillations are dominant. For example, for the top right in Fig 1, we can see clearly the prominent oscillations near (40, 70) in the horizontal direction over about 10 pixels. So the subdomain size will be at least 10 × 10. If there are multiple regions with prominent oscillations but different spatial sizes, then we choose the largest size among them as the subdomain size. Since the function f 0 over the subdomain J is non-periodic in general and the Fourier coefficients computed based on the function values within J are Gibbs-contaminated, we apply the Fourier continuation method to J. Before we use the Fourier continuation method to reduce the Gibbs oscillations near the boundaries of J we first need the function values on a larger subdomain than J as in the 1D case. Here we choose the same number of extra grid points to add to all Thus over this new subdomain J 1 we have a new non-periodic function f : Then we can apply the following two 2D Fourier continuation sparse PA methods to the function f on (n + 2n m ) × (n + 2n m ) uniform grid points over J 1 .

Global 2D Fourier continuation sparse PA method
Now we have a non-periodic function f : And the values of f are given at a set of uniform grid points, ( , with periodicity of L + d on x and y directions. The periodic extension g(z (1) , z (2) ) can be obtained as below where M = (n + 2n m + γ)/2. The unknown coefficientsĝ ðk; lÞ are determined by matching Then we seek gg (a higher number of grid points is used here) by solving the convex optimization problem Heref is zero padding ofĝ ðk; lÞ, that ( F is the analogous 2D Fourier transform matrix of (5). E p x corresponds to E p in (12) in the x direction for each fixed y j , j = 0, � � �, N, and E p y is similarly calculated for the y direction. Finally we obtain the 2D Fourier reconstruction image, f(x i , y j ) = f R (x i , y j ), i = i 1 , i 1 + 1/2, i 1 + 1� � �, i 2 , j = j 1 , j 1 + 1/2, j 1 + 1, � � �, j 2 over the smaller domain J.

Dimension-by-dimension Fourier continuation sparse PA method
The global approach in the previous subsection is computationally slow. Thus we propose to use the dimension-by-dimension approach. That is, we apply the 1D Fourier continuation method using subdomain boundary values for each fixed x i , i = i 1 − n m , � � �, i 2 + n m in y-direction and vice versa. The non-periodic function f(x i , y) on {(x i , y j )jj = j 1 − n m , j 1 − n m + 1, � � �, j 2 + n m } is extended to a periodic function g R (x i , z) on a high resolution grid set {(x i , z I )jI = j 1 − n m , j 1 − n m + 1/2, j 1 − n m + 1, � � �, j 2 + n m + γ−1/2, j 2 + n m + γ, γ is a positive integer}, where if I = j 1 − n m + 1/2 then we have z I ¼ z j 1 À n m þ D 2 . Then the corresponding Fourier coefficients,ĝ R i ðkÞ, to fg R ðx i ; y j Þg j 2 þn m þgÀ 1 j¼j 1 À n m can be found using the discrete Fourier transform where k = −M, � � �, M, M = (n + 2n m + γ)/2. Then we use the convex optimization method using L 1 regularization based on the sparsity of edges to find the reconstruction of the function, g i R ¼ ðg R ðx i ; y j 1 À n m Þ; g R ðx i ; y j 1 À n m þ1=2 Þ; x i 2 þn m � � ½y j 1 À n m ; y j 2 þn m þg �. E p is the sparse PA transform matrix as (12). F is the transform matrix whose elements are We update the values of f on (n + 2n m ) × (2n + 4n m ) grid points over J 1 with these reconstruction data by letting We repeat the same procedure for x-direction and update the values of f over J 1 again. We have f(x, y) on (2n + 4n m ) × (2n + 4n m ) uniform grid points, (x i , y j ), where i = i 1 − n m , i 1 − n m + 1/2, i 1 − n m + 1, � � �, i 2 + n m , j = j 1 − n m , j 1 − n m + 1/2, j 1 − n m + 1, � � �, j 2 + n m . Finally we obtain the 2D Fourier reconstruction, f(x i , y j ), i = i 1 , i 1 + 1/2, i 1 + 1, � � �, i 2 , j = j 1 , j 1 + 1/2, j 1 + 1, � � �, j 2 over the smaller domain J.

Numerical results
In this section, we provide various numerical examples.

Example 1: Shepp-Logan phantom
With this example, we apply both the global 2D Fourier continuation sparse PA method and the dimension-by-dimension Fourier continuation sparse PA method for the reconstruction of the Shepp-Logan phantom image. This reconstructed Shepp-Logan phantom image is found by the following steps. First we use the Shepp-Logan phantom image on the 801 × 801 grid to find its 2D Fourier coefficients,f k x k y , k x , k y = −400, � � �, 400. We use these Fourier coefficients as the exact Fourier coefficients of the Shepp-Logan phantom. Using only a partial number of Fourier coefficients, e.g. k x , k y = −N, � � �, N, N = 42 in this work, we reconstruct the Shepp-Logan phantom via DFT, f N , with which the 2D Fourier continuation method is applied. This experiment mimics an MRI-like acquisition with limited spatial resolution. We compare the results obtained by both the global and dimension-by-dimension methods.
Global method. In Fig 1, Fig 1(a) is the illustration of the Shepp-Logan image on 85 × 85 grids. In Fig 1(b) we choose two sample regions of the reconstructed Shepp-Logan, f N sample region 1 and sample region 2 in the left and right boxes respectively. We apply the global 2D Fourier continuation sparse PA method on these regions. Here we choose λ = 0.01 for both regions. The best value of λ, λ = 0.01 was chosen based on the experiments. The regions in those rectangles with the solid border, J, are the ones where we want to find the reconstruction and the regions in the rectangles with the dashed border, J 1 , are the extended regions with the fixed number of extra grid points to each boundary for the Fourier continuation. The Fig 1(c) and 1(g) are obtained by the following procedure. First we have the 85 × 85 Fourier coefficients that were used to generate the top right subfigure, then apply zero-padding to these Fourier coefficients to generate a new set of 169 × 169 Fourier coefficients. By applying the inverse Fourier transform of these new Fourier coefficients, we create a 169 × 169 image. Finally we take the corresponding regions of the sample regions 1 and 2 in top right. These "double resolution" images have the same matrix size as the reconstruction.
We can see that for both sample regions, J, we eliminated Gibbs oscillations shown in the figures on Fig 1(c) and 1(g) by using the global 2D Fourier continuation sparse PA method.
For both sample regions, the reconstructions (d) and (h) have distinct oscillations (red arrows) near the boundaries of the region J 1 . Since we only consider the smaller region J, these oscillations can be ignored. We can also see that the reconstructions (d) and (h) are blurry near the edges. As we will see in the following section, the reconstruction near edges is too smooth compared to the results by using the dimension-by-dimension method even though the Gibbs oscillations are much reduced.
We also show that the RMSE (root mean squared error) for sample regions 1 and 2, J, are 1.423 and 4.2011 separately.
It took about 250 seconds for the reconstruction with the global approach to be completed with Intel Core i7-3610QM and 2.30GHz.
Dimension-by-dimension method. In Figs 2 and 3: we choose the same sample regions as in Fig 1 and apply the dimension-by-dimension Fourier continuation sparse PA method to these two regions. From top to bottom row is the sample region with 2, 5 and 10 extra grid points added to each boundary. We choose λ = 0.02 for both regions. The value of λ = 0.02 is also chosen by experiments. The images in the first column are the reconstruction of the sample region 1 or 2 of Shepp-Logan phantom image by applying the 1D Fourier continuation sparse PA method in y-direction. In these images, we eliminated the Gibbs oscillations along the vertical direction. The images in the second column are the reconstructed images of the images in the first column by applying the 1D Fourier continuation sparse PA method in x-direction. These images are the reconstructions of the sample region 1 or 2 by applying the dimension-by-dimension Fourier continuation sparse PA method. The images in the third column are the corresponding images to the 169 × 169 exact Shepp-Logan phantom image (ground truth). The images in the forth column are the difference between the reconstructions in the second column and the ground truth in the third column. In the subfigures (d), (h) and (l) in Figs 2 and 3, we observe that the reconstructions near both weak and strong edges are Gibbs-free and sharp. Similar as the global method, the reconstructed images may have oscillations near the subdomain boundaries of the extended domain J 1 . Since we only consider the region J, these oscillations can be ignored. We also observe that the reconstructed images (b), (f) and (j) are much sharper than the reconstructed images (d), (h) and (l) in Fig 1 by the global method near the edges.
We can measure the steepness by visual inspection of the slice image for fixed x or y. In Fig  4, we show the slice image for x = 12 and y = 21 for the sample region 1. From this figure, we can observe that the reconstructed image by the dimension-by-dimension method are much sharper than the reconstructed image by the global method.
By comparing the RMSE in Figs 2 and 3 we observe that if we have 5 or 10 extra points, the RMSE is smaller than the one with 2 extra point. So the reconstruction for 5 or 10 extra points is less noisy in the smooth regions than when we have 2 extra points. By comparing the RMSE in Figs 1 and 2 for sample region 1, we observe that the RMSE for the difference by the dimension-by-dimension method is much smaller than the RMSE for the difference by the global method. Similar conclusion can be found by comparing the RMSE in Figs 1 and 3 for sample region 2. So the reconstruction by the dimension-by-dimension method is more accurate than the global method.
In Fig 5 we observe that the average TV norms for the reconstructions with 2, 5 and 10 extra points added to each boundary are much smaller than the average TV norm for the   2D local Fourier image reconstruction each subdomain and finally stitch all the reconstructions from each subdomain together. By comparing the right image in Fig 6 and left column in Fig 7, we can see that the proposed method eliminates the Gibbs oscillations. By comparing the RMSE, we find that the reconstructions with 5 or 10 extra points are a little better than the reconstruction with 2 extra points.
From Fig 8 we observe that the log 10 errors for 5 or 10 extra points are smaller and the reconstruction is less oscillatory in the smooth region than the case with 2 extra points. Thus when we have 5 or 10 extra points the reconstruction is sharper near the edges and less noisy in the smooth regions than we have 2 extra points.

Examle 2: MRI images
For this example, we provide the MRI image reconstructions with both the global and dimension-by-dimension approaches as in Example 1.
Global method. In Fig 10 we choose two sample regions of the MRI image and apply the global 2D Fourier continuation sparse PA method to these regions. Here we choose λ = 0.00002 for both regions. The chosen value of λ is obtained by experiments. This value is much smaller than the one used for the Shepp-Logan phantom image with the global FC   method. We observe that for both sample regions, J, the reconstructions on the right are blurry near the jumps and that there are oscillations in the smooth part as well. It took about 370 seconds to complete the reconstruction for the right image.
Dimension-by-dimension method. For the dimension-by-dimension method, we first show, in Table 1, the range of λ that provides the best reconstruction results for the given value of n m , the number of extra grid points added to each boundary. The best reconstruction is the reconstruction that yields a sharp reconstruction near jumps while the errors in the smooth region are small. This experiment suggests the range of λ values that can be used when the Fourier continuation method is applied. The lambda values are similar for both the Shepp-Logan phantom image and MRI image when the dimension-by-dimension approach is used. And we already show that the lambda values are very different for the Shepp-Logan phantom image and MRI image while we use the global method. This implies that the dimension-bydimension approach is much more consistent than the global method in terms of choosing the optimal value of λ.
In Figs 11 and 12 we chose two different sample regions, respectively, and applied the dimension-by-dimension Fourier continuation sparse PA method to these regions. For the MRI image, we use the conclusion from the Shepp-Logan example that the reconstructions with 5 or 10 extra points are similar, so we choose the one with 5 extra points for computational efficiency.
The regions in those rectangles with the solid border, J, are the ones where we want to find reconstruction. The images in the left column are the zoomed images of the extended sample regions of the MRI image, which, from top to bottom, are: (1) zoomed images of the extended sample region with 2 (2 is obtained by n/10) extra grid points added to each boundary, (2) zoomed images of the extended sample region with 5 (5 is obtained by n/4) extra grid points added to each boundary. The images in the middle column are the reconstructed images of the images in the left column with the dimension-by-dimension Fourier continuation sparse PA method on the region, J. Here we choose the value of λ according to Table 1. That is λ = 0.0175 for sample with 2 extra grid points and λ = 0.015 for sample with 5 extra grid points. The images in the right column are the zoomed images of the images in the solid rectangle in  2D local Fourier image reconstruction the middle column. From these images we observe that if we have 5 extra grid points the reconstruction is sharper near the edges and less noisy in the smooth regions than when we have 2 extra grid points. For these extended sample regions, the reconstructed images may have oscillations near the boundaries of the extended region J 1 . Since we only need the images in the original region J, these oscillations can be ignored. We also observe that the reconstructed images on the right column in Figs 11 and 12 are much better than the reconstructions on the right column in Fig 10 near both the smooth regions and the edges. By comparing the reconstructions in Figs 10, 11 and 12, we find that the reconstructions with the dimension-by-dimension Fourier continuation sparse PA method are much better than the reconstructions obtained by the global 2D Fourier continuation sparse PA method. Furthermore the dimension-by-dimension Fourier continuation sparse PA method needs about 10 times less computing time than the global 2D Fourier continuation sparse PA method. Thus it is suggested that the dimension-by-dimension approach be used for both accuracy and computational efficiency.
In Fig 13 we split the given MRI image to 8 × 8 subdomains, use the dimension-by-dimension Fourier continuation sparse PA method on each subdomain and finally stitch the reconstructed subdomains together. In Fig 13 we show two reconstruction images: (1) 2 extra grid points (2 is obtained by n/10,n is the number of points on x direction of subdomain) added to each boundary with λ = 0.0175, (2) 5 extra grid points (5 is obtained by n/4) added to each boundary with λ = 0.015. By comparing the regions in the rectangle in Fig 13, we can see that in all the two reconstructions the Gibbs oscillations are eliminated. From these images we observe that if we have 5 extra grid points the reconstruction is sharper near the edges and less noisy in the smooth regions than the case with 2 extra grid points.

Conclusion
In this paper, we extend our 1D domain decomposition Fourier reconstruction method [3] to 2D Fourier image reconstruction. We propose the global 2D Fourier continuation sparse PA method and the dimension-by-dimension Fourier continuation sparse PA method. The global 2D Fourier continuation sparse PA method first divides the 2D image into multiple subdomains and applies the global 2D Fourier continuation to find the 2D periodic extension of the subdomain we are interested in. By finding new Fourier coefficients based on the periodic extension we apply the 2D sparse PA method to obtain the reconstruction. The dimension-bydimension Fourier continuation sparse PA method first divides the 2D image into multiple subdomains and we apply the method in x-and y-directions separately. Finding new Fourier coefficients based on the periodic extension and applying the 1D sparse PA method on these coefficients yield the reconstruction. By splitting the 2D image into multiple subdomains, we obtain sharper reconstruction near both strong and weak edges.
The numerical results in this paper show that the dimension-by-dimension method yields more accurate reconstruction and more efficient than the global method. The dimension-bydimension Fourier continuation sparse PA method for 2D Fourier image reconstruction can be extended to three-dimensional problems by repeating the same procedure in z-direction after applying the method in both x and y-directions.
For our future research we will consider an efficient way to reduce the computational complexity of the proposed method. In this paper we found a range of λ by experiments. We will consider a more systematic way of finding such values. We will also investigate the stability of the global method further and try to devise an efficient 2D global method.