Inversion of incomplete spectral data using support information with an application to magnetic resonance imaging

In this paper we discuss an imaging method when the object has known support and its spatial Fourier transform is only known on a certain k-space undersampled pattern. The simple conjugate gradient least squares algorithm applied to the corresponding truncated Fourier transform equation produces reconstructions that are basically of a similar quality as reconstructions obtained by solving a standard compressed sensing problem in which support information is not taken into account. Connections with previous one-dimensional approaches are highlighted and the performance of the method for two- and three-dimensional simulated and measured incomplete spectral data sets is illustrated. Possible extensions of the method are also briefly discussed.


Introduction
In many application areas, ranging from geophysics to magnetic resonance imaging (MRI), one is confronted with the problem of reconstructing an object, a function, or an image from incomplete Fourier spectral data (see, e.g. [1][2][3]). This is an ill-posed problem in general and very difficult or impossible to solve without any additional information (support or sparsity information, for example). However, by taking a priori information about the object into account, it may be possible to successfully reconstruct the object of interest based on incomplete Fourier data. In compressed sensing (CS), for example, we take into account that the object or image has a sparse representation in some basis and accurate reconstructions are possible provided that the undersampling artefacts are incoherent [4]. As an illustration, in MRI the prototype CS problem consists of minimising the objective function where Ψ is a sparsifying transform, F the (unitary) discrete Fourier transform (DFT) matrix, d the data vector, and S k a diagonal matrix with ones and zeros on the diagonal representing incoherent k-space measurements, where a diagonal entry equal to one corresponds to a point in k-space for which data is available. Often, an additional total variation functional is added to the above objective function as well and the bases that are used for the sparsifying transform Ψ are typically global bases (wavelets, noiselets, etc.) defined over the complete field of view (FoV). With CS, no information about the support of the object is required to successfully image the object of interest. However, optimisation algorithms that minimise objective functions that consist of a leastsquares objective function describing data fidelity and an ℓ 1 objective function that enforces sparsity are generally more complex than algorithms that minimise least-squares objective functions. In addition, in CS a regularisation parameter needs to be determined to balance the data fidelity and sparsifying objective functions in the total objective function.
In this paper, we consider the problem of reconstructing a bounded object from incomplete spectral data in case the support of the object is known. Our motivation comes from magnetic resonance imaging (MRI) where, Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. at least in principle, the support of the object or body part of interest can be determined before an actual clinical scan takes place. During the so-called pre-scan, for example, an image is produced based on not-fully-sampled k-space data. From this image, the FoV is determined and this scan can also be used to estimate the support of the object and its location within the FoV. In general, this approximate support is not exact, of course, but for single body parts (a head, an arm, or a leg, for example) the support can be determined in a fairly straightforward manner using binary maps. Determining approximate supports in case multiple (possibly disjoint) objects or body parts are present within the FoV may be more challenging, but the overall idea remains the same. In this paper, we essentially follow this approach and determine the approximate support of an object (an apple) based on incomplete data that was obtained with a low-field MR scanner [5,6]. We show that image quality improves when this (approximate) support information is included in the reconstruction algorithm. Further examples in which different undersampling patterns are used, are presented as well to demonstrate that directly including support information in the reconstruction algorithm generally improves image quality.
Having support information available, our objective is to reconstruct the image on the FoV from incomplete spectral data. In this case, we can still use CS techniques, of course. In [7], for example, a CS technique is introduced that takes partially known support information into account, while in [8] and [9] the CS problem is formulated as an optimisation problem that aims at finding the solution to the data equation which has the smallest number of components outside of the support.
In this paper, however, we propose a much simpler reconstruction algorithm in case support information is known. Specifically, we build upon the work presented in [2] and propose an image reconstruction algorithm, which is essentially the conjugate-gradient least-squares (CGLS) algorithm applied to the normal equation that corresponds to a space-and frequency-restricted Fourier transform equation. In [2], this approach has been proposed for one-dimensional space-and band-limited Fourier transform problems, while here we generalise this approach to two-and three-dimensional imaging problems. In particular, we consider two-and threedimensional (random) undersampling patterns in the spectral domain that are typically used in MRI to speed up data acquisition and study the performance of the method. We stress that for the two-and three-dimensional spatial support functions and spectral domain undersampling patterns considered here, we have to resort to numerical methods, since no analytical results concerning the eigenfunctions and singular functions of the corresponding truncated Fourier operators are available (as in Slepian-Pollak theory [3,10]). Finally, we also compare our reconstructions to reconstructions obtained via CS techniques and possible extensions of the method (parallel imaging and total variation regularisation) are briefly discussed as well.

Basic equations
Starting point of our analysis is the discrete Fourier transform equation where x is a discrete (vectorised) image function defined on the FoV, F is the unitary discrete Fourier transform (DFT) in one, two, or three dimensions, and d is a vector containing Fourier transform data. In MRI, for example, (2) follows from discretising the signal representation [11] ( ) ( ) ( ) i z is the gradient vector due to the application of the gradient coils of the MR scanner. We assume that the support of an object located within the FoV is known and introduce its indicator or support matrix as a diagonal matrix S x , where a diagonal element is equal to the indicator function of the object applied to the pixel (2D) or voxel (3D) that corresponds to this diagonal element. The support matrix S x is obviously idempotent, that is, it satisfies S S = x x 2 and the image vector x satisfies x = S x x. Furthermore, when spectral data is only available within a subdomain of Fourier space, we introduce a support matrix S k for this subdomain in k-space as well, where the diagonal elements of S k indicate for which points in k-space spectral data is available. The support matrix S k is idempotent and available k-space data is given by the vector S k d.
Having introduced the support matrices S x and S k , we can now formulate our reconstruction problem, which consists of retrieving the image vector x = S x x from the equation where we have introduced the space-and band-limited discrete Fourier transform as A = S k FS x . As pointed out in e.g. [2], in a continuous setting such a reconstruction problem can be formulated in terms of an integral equation of the second kind for the image function restricted to its support. A similar approach can be followed in the general discrete case considered here. In particular, we start with the full DFT of the image and write Applying the inverse Fourier transform F −1 = F H to this equation, we get and by restricting this equation to the object domain by multiplying the above equation on the left by the support matrix S x , we arrive at the equation (8) is the discrete counterpart of the integral equation presented in [2] for one-dimensional problems. If this equation is solved using the Neumann series starting with a vanishing initial guess one arrives at the Papoulis-Gerchberg algorithm or alternating orthogonal projection method [1].
However, as pointed out in [2] for the one-dimensional continuous case, the one-dimensional continuous counterpart of (8) can also be written as a normal equation involving the truncated Fourier transform operator. This approach can be extended to the general case considered here as well. Specifically, given the definition of the space-and band-limited DFT matrix A, it is easily verified that showing that solving (8) is equivalent to solving the normal equation (10). As is well known, any solution that minimises the least-squares objective function satisfies the normal equation (10) as well. From Slepian-Pollak theory [10] we know that in the one-dimensional continuous case, the left and right singular functions of the Fourier transform that is band-and space-limited to the intervals [-k 0 , k 0 ] and [-a 0 , a 0 ] in k-and x-space, respectively, have these intervals as their support. Moreover, the left and right singular functions are complete on L 2 [-k 0 , k 0 ] and L 2 [-a 0 , a 0 ], respectively, and all singular values of the truncated Fourier transform belong to the interval [0, 1] with clustering occurring around zero and one. Some of these one-dimensional spectral results have been extended to higher dimensions [3] (in two dimensions to functions with circular support in ordinary and k-space, for example), but no such theory exists for the truncated two-or three-dimensional discrete Fourier transforms considered here with general x-or k-space indicator functions.
Therefore, let us first consider computing the singular value decomposition (SVD) of matrix A given by A = UΣV H , where the columns u i of U and the columns v i of V are the left and right singular vectors of matrix A, respectively, and Σ is a diagonal matrix with the nonnegative singular values σ i of A on its diagonal arranged in decreasing order. With the SVD of matrix A at our disposal, let v v v be the space spanned by the first k right singular vectors that correspond to the first k largest singular values σ 1 σ 2 ... σ k > 0. We can then take the kth truncated SVD solution as an approximate solution to (10). Care must be exercised when selecting k, however, since a poorly chosen k may lead to an overly smooth reconstruction or the approximate solution x k is heavily affected by noise that is present in the data. For reconstruction problems encountered in practice, however, computing the SVD of matrix A comes at prohibitively high computational costs. Computing matrix-vector products with matrix A, on the other hand, can be carried out at 'FFT speed' and we therefore resort to iterative methods that solve the reconstruction problem. Specifically, instead of the approximate solutions of (12), we take is the kth Krylov subspace generated by A H A and vector A H d. The CGLS algorithm [12] that starts with a vanishing initial guess produces approximations x k that satisfy (13) and we therefore use this algorithm to solve the reconstruction problem. The reason for taking the x k of (13) as approximate solutions comes from the observation that the first Krylov subspace vector (right-hand side vector) predominantly contains contributions from right singular vectors that correspond to the largest singular values. Furthermore, matrix A has many zero singular values and since in CGLS the largest singular values are typically approximated first, we expect that x k mainly contains information about the right singular vectors that correspond to the first k singular values. In fact, in case of noisy data, the regularising effects of CGLS may be used to obtain stable reconstructions by terminating the iterative process after a certain number of iterations (semiconvergence, see e.g. [13]). For completeness, we mention that the approximations x k of (13) may also be generated by the LSQR algorithm [12].
Finally, we note that if we add a sparsifying regulariser to the objective function of (11), we end up with the prototype objective function of CS (1). If the image has a sparse representation in some basis and the imaging artefacts are incoherent then CS may accurately reconstruct the image without any support information. However, when this information is available then CS techniques may not be necessary and solving the normal equation (10) with CGLS and a vanishing initial guess may already give satisfactory reconstruction results. In the following section we present reconstruction results for k-space undersampling patterns typically used in MRI, which illustrate that solving the normal equation with CGLS indeed may be sufficient in these cases and no CS techniques are necessary.

Image Reconstruction
To illustrate the performance of the CGLS imaging algorithm in case support information is known, we consider reconstruction problems for noise-free and noisy 2D and 3D data sets. In our 2D experiments we use a Shepp-Logan phantom with 64 × 64 pixels [14] as the model solution and a brain image from the Kirby 21 dataset [15] with 256 pixels in each spatial direction. These phantoms and their supports are shown in figure 1. The support of the 2D Shepp-Logan phantom consists of 1988 out of a total of 4096 pixels (≈49%), and the support of the 2D brain image has a size of 23 892 out of a total of 65 536 pixels (≈36%).
For the model solutions shown in figure 1, we consider six different undersampling patterns with Fourier or k-space support functions illustrated in figure 2. The first support function ( figure 2(a)) is a square undersampling pattern that leaves out the edges of k-space. The second support function (Figure 2(b)) is a pattern of random lines with all center lines being sampled, while the third pattern (figure 2(c)) consists of completely random lines (so not all center lines are sampled). An undersampling pattern consisting of random points ( figure 2(d)) is also considered, along with a radial pattern and a spiral pattern (figures 2(e) and (f), respectively). In all cases, the undersampling factor is approximately equal to two. We note that in MRI, some of these patterns are more difficult to implement in practice than others due to hardware limitations. Additionally, we remark that in practice, when using a radial or spiral pattern, the sampled points do not necessarily fall on a Cartesian grid. To make the application of an inverse Fourier transform to k-space possible, an iterative process called gridding [16,17] may be used to transform the measured data into a Cartesian format. In this paper, however, we do not take these limitations into consideration and simply assume that we have access to undersampled k-space data in Cartesian format. Figure 3 shows the reconstruction results for the Shepp-Logan phantom and the six different undersampling patterns of figure 2. Specifically, the model solution is shown in the left column, while the reconstructions obtained by simply applying an inverse FFT (IFFT) to incomplete spectral data are shown in the second column. The reconstructions that are obtained using CGLS starting with a vanishing initial guess are shown in the fourth column of figure 3. These results were obtained after 1000 CGLS iterations or when the normalized residual dropped below a tolerance of 10 −10 . We observe that for all cases, CGLS with known support information yields images of improved quality compared to images obtained by simply applying an IFFT. This is also reflected in the peak signal-to-noise ratio (PSNR) values presented in table 1.
For the brain image, the reconstruction results are shown in figure 4 and detailed views of the reconstructions within the domain indicated by the white square in figure 1(c) are shown in figure 5. Comparing the reconstruction results obtained with an IFFT with the reconstructions obtained with CGLS, we again observe that including support information significantly improves the quality of the reconstructions. To quantify this statement, we compute the PSNR values of the various reconstructions shown in figure 4 and the results are presented in table 2. We observe that for all undersampling patterns the PSNR values of the reconstructions obtained with CGLS and support information included are higher than the PSNR values of the reconstructions obtained with an IFFT without support information. We note that the CGLS algorithm performs particularly well for a random points undersampling pattern (The norm of the error of the reconstruction result is in the order of 10 −8 for the brain image). Unfortunately, such a sampling pattern is very difficult or impossible to realise during practical MR scans. Furthermore, in the case of a spiral or radial undersampling pattern, our algorithm using support information significantly increases the PSNR of the reconstruction.

Compressed sensing
The image reconstruction problem may also be solved with CS techniques, of course. No support information is required in this case, but CS imaging algorithms are generally more complex than the CGLS approach considered in this paper.
As an illustration, let us consider the CS problem that consists of finding the image function x that minimises the objective function for the undersampling patterns of figure 2. Here, Ψ is the Daubechies wavelet operator [18] and λ is a regularisation parameter, which is determined in a heuristic manner through numerical experimentation. No  points undersampling patterns. In these cases, CGLS outperforms CS and, as mentioned above, especially the CGLS reconstruction for a random points undersampling pattern is of very high quality.
From figure 4 and 5 and table 2 we observe that for the Kirby data set, CGLS in general again produces reconstructions of a similar quality as the reconstructions obtained with CS, except that here the PSNRs of all CGLS reconstructions are larger than the PSNRs of the corresponding CS reconstructions. A random points undersampling pattern, in particular, produces a CGLS reconstruction that is clearly superior to the CS reconstruction. Moreover, for a random lines undersampling pattern, CGLS produces a reconstruction that is significantly better than the reconstruction obtained with CS as well with a PSNR that is almost twice as large as the PSNR of the CS reconstruction.

Noisy measurements
To demonstrate the performance of the CGLS reconstruction method in case of noisy data, we again consider reconstructing the Shepp-Logan and Kirby phantoms from their corresponding incomplete k-space data, but now the data is contaminated by noise with an SNR of 50. For both phantoms, we restrict ourselves to a spiral undersampling pattern, since the effects of noise are similar for the other undersampling patterns and noisy data is also considered in section 3.3, where we apply the proposed CGLS algorithm to measured data.
In figure 6, the model solutions (first column), the reconstructions obtained with an inverse FFT (second column), and the reconstructions obtained with the CGLS algorithm with support information included (third column) are shown for the Shepp-Logan phantom (first row), the Kirby head model (second row), and for the region indicated by the white square in figure 1(c) (third row). These results clearly show that the CGLS

Reconstructions based on experimental low-field MRI data
In this section, we apply the proposed reconstruction method to a 3D data set, which was acquired using a spinecho sequence on the low-field MRI scanner described in [6]. The scanning parameters of the measurement are reported in table 4. An apple was placed inside the scanner and all imaging was carried out using a single receive coil. We note that the signals that were obtained have an SNR that is much lower than the SNR of signals measured in typical commercial MR scanners (in this case, the SNR is about 6). Applying a 3D inverse Fourier transform to fully sampled k-space data results in a 3D reconstruction of the apple and some slices of this reconstruction are shown in the left column of figure 7. These reconstructions serve as model solutions, since we obviously do not have a perfect or high SNR model solution available in this case.  To illustrate the performance of CGLS, we again consider a k-space undersampling pattern of random lines with the center of k-space completely sampled. The resulting undersampling factor is equal to two and the reconstruction results obtained by simply applying an inverse FFT to the incomplete data set and after 100 iterations of the CGLS algorithm are shown in the second and third column of figure 7, respectively. We observe that an image is obtained which visually resembles the original full k-space solution. The PSNR values in table 5 also show that using support information increases the PSNR from 24.02 (IFFT reconstruction) to 61.92 (CGLS reconstruction).

Conclusions
In this paper we discussed an imaging method that can be applied if the support of an object within a certain FoV is known and its spatial Fourier transform is only known on a certain k-space undersampling pattern. We demonstrated that the CGLS algorithm applied to the corresponding truncated Fourier transform equation produces reconstructions that are essentially of a similar quality as reconstructions obtained by solving a standard CS problem in which support information is not taken into account. In particular, for the Shepp-Logan phantom and a range of two-dimensional k-space undersampling patterns, the CGLS algorithm produces reconstructions with PSNR values that are generally slightly lower than the PSNR values of the corresponding CS reconstructions. However, for a realistic head model the PSNR values of the CGLS reconstructions are typically larger than the PSNR values of the CS reconstructions. Specifically, in the case of a random points undersampling pattern, the CGLS reconstructions of the Shepp-Logan phantom and the head model have significantly larger PSNR values than their CS counterparts. For both phantom models, the PSNR values of the CGLS reconstructions are also larger than the PSNR values of the reconstructions obtained by simply applying an IFFT to the available incomplete data sets.
In 3D, an improvement in the reconstructions was also observed when support information is taken into account. The CGLS algorithm was applied to a measured data set obtained with a low-field MR scanner. For a random lines undersampling pattern the method was able to provide reconstructions with sufficient detail and PSNR values that are significantly larger than the PSNR values of the reconstructions obtained via a standard IFFT. In conclusion, when support information about the object is available, a straightforward application of the CGLS algorithm to a truncated Fourier transform equation definitely improves simple IFFT-based reconstructions and, at least for the undersampling patterns considered here, generally provides reconstructions with a quality similar to reconstructions obtained via standard CS techniques.
Carrying out a Fourier reconstruction using a reduced number of data points invariably leads to imaging artefacts. In this paper, we considered basic CGLS and CS reconstruction algorithms to address this problem. More advanced CS techniques can be used as well, of course, (see [19][20][21][22][23][24][25][26][27][28][29][30], and [31], for example), but the same holds true for the CGLS algorithm. For example, in an MR setting, the use of multiple receive coils can be included in the reconstruction process and parallel imaging techniques analogous to SENSE [32] and GRAPPA [33] may be developed. Furthermore, total variation regularisation may be included in the reconstruction process as well either in an additive or multiplicative manner [34]. Specifically, suppose that we have N receive coils available with coil sensitivity maps represented by the diagonal matrices C i , i = 1, 2,K,N. The data Table 5. PSNR of the 3D apple reconstructions obtained using a simple inverse Fourier transform and using CGLS with support information included.

PSNR (IFFT)
PSNR (Support) Apple image 24.02 61.92 collected by the ith coil is given by S k FS x C i x = d i , i = 1, 2,K,N, and all data can be combined into a single data equation as Ax = d with A = (I N ⊗ S k FS x )C, where I N denotes the identity matrix of order N, With multiplicative total variation regularisation included, the imaging problem is posed as an optimisation problem with an objective function given by

TV
where P · P 2 is the 2-norm and F TV (x) is a multiplicative total variation objective functional (see [34] for details). Regularisation can also be included in the usual additive manner, of course, but multiplicative regularisation has the advantage that no effective regularisation parameter needs to be determined and computed. Future work will focus on the full development of the above image reconstruction technique.