Group Sparse Reconstruction of Multi-Dimensional Spectroscopic Imaging in Human Brain in vivo

Four-dimensional (4D) Magnetic Resonance Spectroscopic Imaging (MRSI) data combining 2 spatial and 2 spectral dimensions provides valuable biochemical information in vivo; however, its 20–40 min acquisition time is too long to be used for a clinical protocol. Data acquisition can be accelerated by non-uniformly under-sampling (NUS) the ky − t1 plane, but this causes artifacts in the spatial-spectral domain that must be removed by non-linear, iterative reconstruction. Previous work has demonstrated the feasibility of accelerating 4D MRSI data acquisition through NUS and iterative reconstruction using Compressed Sensing (CS), Total Variation (TV), and Maximum Entropy (MaxEnt) reconstruction. Group Sparse (GS) reconstruction is a variant of CS that exploits the structural sparsity of transform coefficients to achieve higher acceleration factors than traditional CS. In this article, we derive a solution to the GS reconstruction problem within the Split Bregman iterative framework that uses arbitrary transform grouping patterns of overlapping or non-overlapping groups. The 4D Echo-Planar Correlated Spectroscopic Imaging (EP-COSI) gray matter brain phantom and in vivo brain data are retrospectively under-sampled 2×, 4×, 6×, 8×, and 10× and reconstructed using CS, TV, MaxEnt, and GS with overlapping or non-overlapping groups. Results show that GS reconstruction with overlapping groups outperformed the other reconstruction methods at each NUS rate for both phantom and in vivo data. These results can potentially reduce the scan time of a 4D EP-COSI brain scan from 40 min to under 5 min in vivo. Algorithms 2014, 7 277


Introduction
Magnetic Resonance Imaging (MRI) exploits the resonant frequency of 1 H protons within water in vivo to generate anatomical images of the human body.Magnetic Resonance Spectroscopic Imaging (MRSI) is a similar imaging technique to MRI; however in lieu of the resonant frequency of 1 H protons in water, the resonant frequencies of 1 H protons in metabolites including, but not limited to, creatine (Cr), choline (Cho), glutamate (Glu), glutamine (Gln), lactate (Lac), aspartate (Asp), and N-acetyl-aspartate (NAA) are used to generate a metabolic image instead [1,2].Each of these metabolites has a unique resonance spectrum in vivo caused by their chemical environment and covalent bonding structure that can be used to identify and quantify their concentrations within each voxel [3].Using this concentration information, the biochemical compositions of healthy and diseased tissues can be determined without the need for invasive biopsies and the altered metabolism of cancers can be detected [4].
The four-dimensional (4D) Echo-Planar Correlated Spectroscopic Imaging (EP-COSI) is an MRSI pulse sequence that allows for the simultaneous acquisition of two spatial and two spectral dimensions, (k y ,k x ,t 2 ,t 1 ), in one scan in vivo [5].This 4D sequence provides a 2D spectrum for each voxel in a slice, which differs from other MRSI sequences that only provide a single spectral dimension per spatial voxel [6,7].The overlap of resonance peaks within a single spectral dimension is a major impediment to identifying individual metabolites and the increased spectral dispersion offered by a second spectral dimension can disentangle complex overlapping spectral peaks [8].Two spectral dimensions have been shown to increase the specificity and sensitivity of tumor grade classification when used with dynamic contrast enhanced (DCE)-MRI in the breast [9].However, 4D EP-COSI acquisitions are slow compared with most MRI sequences and can require up to 40 min for a typical scan, which is too long to be used on a routine clinical basis.
The 4D EP-COSI data acquisition is a rasterized scan that acquires a 2D spatio-temporal plane, k x − t 2 , from the 4D volume during each TR.The second spatial (k y ) and spectral (t 1 ) dimensions are incrementally acquired between rasters until the entire 4D volume is sampled.k y is acquired using standard phase encoding techniques from MRI and t 1 is acquired as a series of phase shifted 1D spectra [5].The 4D EP-COSI data acquisition can be accelerated in vivo by non-uniformly under-sampling (NUS) the incrementally acquired dimensions, k y and t 1 .However, NUS produces artifacts in the spatio-spectral mixed domain that can be removed by reconstructing the missing samples in the k y − t 1 plane through non-linear, iterative reconstruction [10,11].
Previous work has accelerated MRI acquisition by under-sampling k y in k-space and reconstructing the missing samples with Compressed Sensing (CS) reconstruction [12,13].CS has been applied to under-sampled, spatio-temporal and spatio-spectral mixed domains in dynamic MRI [14] and 3D spatial, 1D spectral MRSI [15].The feasibility of under-sampling the mixed-domain k y − t 1 plane of a 4D MRSI dataset and reconstructing the missing samples with Total Variation (TV) denoising or Maximum Entropy (MaxEnt) reconstruction has been shown previously [16,17].NUS rates as low as 5× were demonstrated in vivo and showed that it is possible to accelerate the acquisition of 4D MRSI data down to a clinically acceptable 5-10 min scan time.
The current work derives a solution to the Group Sparse (GS) reconstruction problem for 4D MRSI data within the Split Bregman iterative framework and compares the results with CS, TV, and MaxEnt reconstructions of 4D EP-COSI phantom and in vivo brain data under-sampled at 2×, 4×, 6×, 8×, and 10× [18,19].The Split Bregman algorithm has been modified to allow for arbitrary coefficient grouping patterns of overlapping or non-overlapping groups with no a priori assumptions on the grouping patterns used.The Split Bregman algorithm has been previously used for TV-based reconstruction of NUS MRSI data [16] and for multi-channel reconstruction of NUS MRI data where the multiple measurement vector (MMV) problem was solved by extending the algorithm to accommodate row-wise grouping of jointly sparse samples [20].GS reconstruction has previously been applied to the under-sampled k y − t planes of 2D cine and perfusion cardiac MRI and provided superior results to CS [21,22].
In this article, overviews of the Split Bregman algorithm and GS reconstruction problem are provided in Sections 2.2 and 2.3, respectively.The derivation of the solution to the GS reconstruction problem within the Split Bregman framework is provided in Section 2.4.A brief overview of the CS, TV, and MaxEnt reconstruction problems is provided in Section 2.1, and a description of the non-linear algorithms and parameter settings used to solve them is given in Section 3.2.The experimental setup for the 4D EP-COSI gray matter brain phantom and in vivo brain scans is provided in Section 3.1.Results from retrospectively under-sampled phantom and in vivo datasets then reconstructing the NUS k y − t 1 plane with CS, TV, MaxEnt, and GS with overlapping and non-overlapping groups, are shown in Sections 4.1 and 4.2.These experiments show that GS reconstruction with overlapping coefficient groups provides reconstructions with lower metabolite Root Mean Square Errors (RMSE) than the other methods used, as well as qualitatively superior spectra down to 8× NUS in both phantom and in vivo datasets.This acceleration translates to a clinically acceptable 5 min 4D EP-COSI brain scan that previously took 40 min to complete.

CS, TV, and MaxEnt Based MRSI Reconstruction
The under-sampling of 4D MRSI datasets and subsequent non-linear reconstruction by CS, TV, or MaxEnt have been described by previous work and will only be covered here briefly [16,17].Each of these methods can be used to iteratively reconstruct missing samples from an NUS MRSI dataset by enforcing a fidelity constraint in the sample domain and either a maximum sparsity, minimum differences, or maximum entropy signal model in some transform domain [10,23,24]: where u ∈ C N is the reconstructed spatial, spectral-domain (Y, X, F 2 , F 1 ) data, F is the 4D Fourier operator, R ∈ R M xN is the under-sampling mask that determines which samples were acquired in the data, and σ is the standard deviation of noise in f .The choice of φ(u) depends on the signal model used during reconstruction.The transform sparsity signal model used in CS uses φ(u) = Ψu 1 , which is the l 1 norm of the reconstructed data in some transform domain, Ψ, and attempts to remove NUS artifacts by reducing the amplitude of each transform coefficient and increasing the sparsity of Ψu [25].The minimum TV signal model uses φ(u) = ∇ y • Ψu + ∇ F 1 • Ψu 1 , where ∇ y and ∇ F 1 are anisotropic differential operators applied to Ψu along the y and F 1 dimensions, respectively, and removes NUS artifacts by minimizing the differences between adjacent transform coefficients.Because the mixed-domain k y − t 1 plane of a 4D EP-COSI dataset is self-sparse, no sparsifying transforms are applied to u for the CS and TV reconstruction so Ψ = I [16].The MaxEnt signal model uses φ(u) = −S1 /2 , which is the spin-1 /2 entropy derived from the statistical distributions defined by the density matrices of spin-1 /2 nuclei, such as 1 H [26].It removes NUS artifacts by modeling each transform coefficient in the spatial-spectral domain as the byproduct of an ensemble of 1 H spins whose phase dispersion determines the signal amplitude; in-phase spin ensembles have high signal amplitude and low phase entropy, while high entropy ensembles have low phase coherence and low signal amplitude.Because Shannon entropy is derived for discrete processes that result from a Poisson-based distribution, it cannot be used for reconstructing NMR data, which do not follow a Poisson distribution [26,27].

Split Bregman Algorithm
The Split Bregman algorithm is from the class of Alternating Direction Method of Multipliers (ADMM) that split a constrained problem into a sequence of simpler unconstrained sub-problems [19,28].It differs from continuation-based methods by keeping the values of any Lagrange multipliers fixed between iterations and modifies the data instead.This has the benefits of increased numerical stability and a lower dependence on the initial choice of Lagrange multiplier values.If we wish to use the Split Bregman algorithm to solve a problem of the form: where E(u, z) is a convex, non-differentiable function and H(u, z) can be assumed to be of the form Au = z, we first convert it to an unconstrained problem: where conventional continuation-based methods would increase λ → ∞ as H(u, z) → 0 to find a solution; however, we apply the Bregman distance relaxation to Equation (3) and split it into sub-problems instead.The Bregman distance for the function E(u) at the point u k is: where the Bregman parameter p k u is the sub-gradient of E(u) at u k .The Bregman distance and iteration scheme are then applied to Equation (3) by defining the u, z, and Bregman parameter updates at each iteration over k: λ is never increased and under fairly weak assumptions on E(u, z) and H(u, z), ∇ u H and ∇ z H → 0 as k → ∞ so that p k+1 → p k and the Bregman parameters converge [19].However, because H(u, z) is of the form Au = z, Equation ( 5) can be simplified to the equivalent problem [29]: where b k+1 is a Bregman parameter that ensures Au → z as b k+1 → b k converges without increasing λ and sacrificing stability.Because E(u, z) is convex and non-differentiable, Equation ( 6) is split into its u and z sub-problems, which are solved separately at each iteration, decoupling u from z [28]:

Group Sparse Reconstruction
GS reconstruction is an extension of CS that exploits the correlations among adjacent transform coefficients caused by their structured sparsity [18,30,31].Structured sparsity, or block sparsity as it is also known, is the tendency of large transform coefficients to be adjacent to each other and form clusters rather than be evenly distributed across the transform domain.In GS reconstruction, adjacent transform coefficients are reconstructed together in groups rather than individually, as is done in CS.By reconstructing groups of coefficients, the GS signal model correlates individual transform coefficients with their neighbors, allowing them to influence each other during reconstruction.Because of the structured sparsity relationships among transform coefficients, it has been shown the GS signal model provides results equal to or better than CS in an RMSE sense, which does not exploit these relationships [30,32].
GS reconstruction can be formed as a constrained convex optimization problem that uses the l 2,1 -norm as the objective function instead of the l 1 -norm used in CS.The l 2,1 -norm over a vector u ∈ C N is defined as: where u gi ∈ C P = {u j , u k ...u n } and j, k, n ∈ {S = 1...N, P ≤ N } is one of the L groups of u.Because the u gi groups may overlap, u 2,1 is non-convex and variable substitution must be used on u to ensure that GS reconstruction is a convex optimization problem: where z is defined as z = Gu and G ∈ R LP ×N is the group matrix of 1s and 0s that determines which transform coefficients from u belong to each group in z [33].G is defined as having a single 1 per row so no z gi overlap but each z gi may contain any number of coefficients from u.The remaining terms are unchanged from Equation (1).Because each transform coefficient from u may be within more than one z gi group and separate reconstructions must be created for each version of that coefficient, the set of transform coefficient groups, z, may contain more points than u and is an over-complete reconstruction of u.
The GS problem reduces to the CS problem when each group contains one coefficient, thus making G = I and z = u.As implemented for these experiments, the grouping pattern must offer complete coverage of the entire 4D MRSI volume so that each spatial, spectral-domain coefficient is contained in at least one group.

Split Bregman Based Group Sparse Reconstruction
The GS reconstruction problem is generally considered difficult to solve due to its mixed norm structure and the non-smooth nature of the l 2,1 objective function.However, by applying Split Bregman iterative relaxation, we can solve Equation ( 9) for the optimal MRSI reconstruction to within the standard deviation of noise with ADMM convergence guarantees [19].In order to apply the Split Bregman algorithm, Equation ( 9) is formed as an unconstrained problem and both constraints converted to equality: where Relaxing Equation ( 9) into Equation (10), forms E(u, z) as a combination of the objective function and fidelity constraint through the Lagrange multiplier, µ, which would normally cause the reconstruction result to depend heavily on the initial value for µ; however, as will be shown, the Bregman iterative framework reduces the dependence on µ by treating the sampled data, f , as a Bregman parameter and updating it accordingly throughout the reconstruction [29].By following the process defined in Equations ( 5)- (7) and splitting E(u, z) into its l 2,1 and l 2 components, we can derive the u and z sub-problems and a Bregman parameter update that solves the unconstrained problem in Equation (10): The u k+1 sub-problem is now a function of u with z k treated as a constant, and the z k+1 sub-problem treats u k+1 similarly.The z k+1 sub-problem is non-differentiable; however, its equivalent problem can be solved as: where gshrink is group-wise shrinkage over each group, instead of individual coefficients, and has the closed form solution: where x gi = (Gu k+1 ) gi + (b k z ) gi .The u k+1 sub-problem is differentiable; therefore, optimality conditions can be derived for u and simplified to: where G G is a diagonal matrix with each (G G) ii being the number of groups that contain transform coefficient u k+1 i [33].If each transform coefficient is in the same number of groups, hence G G is a multiple of I, the left-hand side of Equation ( 15) is circulant and can be easily inverted by the Fourier transform [19]: where If each transform coefficient is not in the same number of groups and G G is not a multiple of I, u k+1 in Equation ( 15) can be solved for by the Gauss-Seidel method or other iterative solvers for linear systems.
Iterating the u and z sub-problems and Bregman update over k in Equation ( 12) solves the unconstrained problem from Equation (10), where the Lagrange multiplier, µ, was introduced; however, to solve the constrained problem in Equation ( 9) and reduce the dependence of the reconstruction result on the initial value of µ, an outer iteration is completed over i that updates f based on changes to the reconstruction, u k+1 : where f is the initial value of the sampled data.By replacing f with f i in Equation ( 12) then applying Equation ( 17) during each outer iteration, f i+1 → f i and RFu k+1 → f to within the standard deviation of the noise, as the algorithm converges to the solution of the constrained problem [29].
The pseudocode for the Split Bregman GS MRSI reconstruction algorithm is shown in Algorithm 1 and incorporates Equations ( 12)- (17).The inner loop iterates over k and updates u and z, and the outer loop iterates over i and updates the sampled data f i as a Bregman parameter.
The number of iterations over k is application-dependent, but for the 4D MRSI GS problem, we use M = 15 iterations over k for each iteration over i and iterate over i until the normalized residual error between u k+1 and f i is less than 10 −6 .
In addition to the phantom data comparisons, each reconstruction method was also evaluated on retrospectively under-sampled in vivo 4D EP-COSI brain data acquired with the same sequence parameters, except for TR/TE = 1.7 s/23 ms and outer volume suppression (OVS) prior to each TR to suppress the scull marrow fat resonances.The brains of two healthy human volunteers aged 23 and 36 years were scanned using the Siemens 3T Trio scanner with a 12-channel head coil, using the same post-processing and retrospective NUS sample masks as the phantom scans prior to reconstructing with CS, TV, MaxEnt, and GS.

MRSI Scan Reconstruction
Separate reconstructions of each NUS phantom and in vivo scan were performed using CS, TV, MaxEnt, GS with non-overlapping groups (GS 1 ), and GS with 50% overlapping groups (GS 2 ) [34].The grouping patterns were within the spectral domain and comprised of equal-sized groups of (F 2 , F 1 ) = (8, 4) coefficients as illustrated by Figure 1, which shows the coefficient grouping support of a single data point among nine adjacent groups for the non-overlapping and 50% overlapping cases.The CS and TV reconstruction problems were solved using the Split Bregman algorithm [19] and MaxEnt was solved using the Cambridge algorithm [17,35].Because GS simplifies to the CS problem with group sizes of one coefficient, (F 1 , F 2 ) = (1, 1), the CS and GS problems were solved using the same code base, but with different groupings.In order to simplify any necessary parameter tuning, the optimal Lagrange parameters used by the CS reconstruction were used by the GS reconstructions, but scaled by 1  G , where G is the number of coefficients per group.The optimal TV, CS, and GS Lagrange parameters were chosen by empiric determination and because the Split Bregman algorithm reduces the dependence of the reconstruction results on the initial value of the Lagrange parameters used, the same values were used for the phantom and in vivo scans at each of the NUS rates tested.The Lagrange parameters used by each reconstruction type are shown in Table 1.
The MaxEnt algorithm converged when the normalized difference between the gradients of the entropy (S) and residual (C) were parallel: ∇S / S 2 − ∇C / C 2 ≤ 0.001 [17].The CS/GS and TV implementations of the Split Bregman algorithm used 15 inner loops per outer loop and converged when the normalized residual error between u k+1 and f i was less than 10 −6 .The maximum number of outer loops was limited to 25 but was never reached.

Gray Matter Brain Phantom
Figure 2 shows a select 2D COSY spectrum from the central voxel of a 4D EP-COSI gray matter phantom scan that was retrospectively under-sampled 8× using the phantom mask in Figure 3 and reconstructed by CS, TV, MaxEnt, GS 1 , and GS 2 .Because of the homogeneity of the phantom spectra, a single representative 2D spectrum was chosen to illustrate the qualitative effectiveness of each reconstruction method.The 1D projections of the F 1 and F 2 dimensions are shown to the right and above each 2D spectrum, respectively.The spectrum from the fully sampled phantom scan is shown in the top left of the figure and the 8× NUS spectrum with zeros in place of missing samples is shown to its right.The contour levels of each spectrum are the same as those used for the fully sampled scan to allow noise levels and peak heights to be compared between reconstructions.As can be seen in the 8× NUS spectrum, the broad point spread function (PSF) of the under-sampling mask caused the large diagonal peaks of NAA, Cr, Cho, and mI to alias along the F 1 dimension and obscure the much smaller NAA, Glx, mI, and Asp cross-peaks, which have also been spread along F 1 by convolution with the PSF.Because the t 2 dimension was fully sampled, the F 2 projections of the fully sampled, 8× NUS, and reconstructed spectra have similar peak amplitudes and line widths; however, the t 1 under-sampling caused the spectral peaks shown in the fully sampled F 1 projection to collapse in the 8× NUS projection.The F 1 projections of reconstruction method show reasonable agreement with the fully sampled dataset, indicating the peak energy of major diagonals, such as NAA, Cr, Cho, and mI, has been deconvolved from the broad PSF by each of the reconstructions.
There are noticeable differences among the reconstructions that should be highlighted.The noise floors of the CS and TV reconstructions are higher than all of the other spectra, which makes identifying cross peaks difficult.The CS and TV reconstructed diagonal peaks have lost their Lorentzian line shape and the cross peaks are present but "choppy" in appearance, further indicating issues with proper line shape reconstruction.The MaxEnt and GS reconstructions have noise floors equal to or lower than the fully sampled dataset.Cross peaks and diagonals are not "choppy" in appearance, and most of the cross peaks shown in the fully sampled dataset can be identified.The MaxEnt reconstruction has lower peak amplitudes than the GS reconstructions and possibly TV, but its low noise floor means that low amplitude peaks are not obscured, which was confirmed by adjusting the contour levels down.The peak amplitudes and line shapes of the GS reconstructions are the most similar to the fully sampled dataset from all of the reconstructions, as illustrated by the low amplitude NAA and Glx lower and upper cross peaks (LCP, UCP), whose complex multiplet structures have been properly reconstructed with reasonable line shapes and widths.
4D EP-COSI datasets are self-sparse along the Y −F 1 plane [16].There are large regions between the metabolite peaks that are dominated by noise, therefore Root Mean Square Errors (RMSE) calculated over the entire dataset can be heavily biased by changes to the noise distribution from non-linear reconstruction methods like CS, TV, MaxEnt, and GS.Previous work has demonstrated that constraining the RMSE calculations to the metabolite peaks and ignoring the noise regions is a suitable quantitative metric of reconstruction accuracy that coincides with qualitative observations because it measures amplitude, line width, and line shape simultaneously with minimal influence from the noise distribution [17].Therefore, the mean metabolite peak reconstruction RMSEs with respect to the fully sampled data of six EP-COSI gray matter brain phantom scans versus NUS rate are shown in Figure 3 along with the 2×-10×NUS masks used to under-sample the k y − t 1 planes of the phantom and in vivo datasets used in these experiments.The RMSEs were calculated for each metabolite peak using the ppm ranges shown in Table 2.The means were calculated from the central 4 × 4 VOI of the six phantoms scans used, yielding 96 voxel measurements per metabolite peak.The standard deviations are not shown because they were 2-3 orders of magnitude lower in value than the means.
As can be seen in Figure 3, the mean reconstruction RMSEs for the large diagonal peaks, NAA 2.0 ppm, Cr 3.0 ppm, and Cr 3.9 ppm, are shown at the top of the figure, while smaller diagonal peaks and cross peaks are shown in the middle and bottom, respectively.For the large diagonals, all of the reconstruction RMSEs follow similar trends at each NUS rate, which is in agreement with the qualitative results from Figure 2 that indicated each of the reconstruction methods was able to accurately reconstruct the diagonal peaks.However, for the smaller diagonals and cross peaks, the general trend at lower NUS rates is for CS to have the highest RMSE, followed by TV, MaxEnt, GS 1 , and GS 2 , but at higher NUS rates, MaxEnt has RMSEs equal to or lower than GS 1 .The RMSEs of GS 1 for Cho and mI do not follow the general trends of the large diagonal peak RMSEs, nor those of the smaller diagonals and cross peaks; the Cho RMSEs for GS 1 are higher than all other methods for each NUS rate and one of the highest for mI.However, the RMSEs for GS 1 of the other metabolites are some of the lowest of the methods tested.
Table 2. F 1 ,F 2 ppm ranges used to calculate the metabolite RMSEs shown in Figure 3.These results demonstrate that GS reconstruction has the lowest RMSEs for each metabolite peak and NUS rate, with either GS 2 or GS 1 providing the lowest mean RMSE.The mean RMSEs for GS 2 are consistently lower than GS 1 , indicating that overlapping groups may be a more reliable method of metabolite peak reconstruction than non-overlapping groups.

In Vivo Brain
The reconstruction results for an in vivo 4D EP-COSI brain scan of a healthy 23 year old male are shown in Figure 4.A representative 2D COSY spectrum was chosen from a central 2 × 2 × 2 cm 3 voxel of the brain that shows each of the metabolites of interest with low fat and water contamination, as illustrated by the green square on the T1-weighted image in the bottom left of Figure 4.The fully sampled spectrum in the top left corner of the figure was under-sampled 8× using the sample mask from Figure 3, then reconstructed by CS, TV, MaxEnt, GS 1 , and GS 2 as shown.
As is seen in Figure 2, the 8× NUS spectrum shows significant aliasing of the diagonal peaks along F 1 and a collapse of the peaks within the F 1 projection.The CS and TV reconstructions have higher noise floors than the other reconstructions and narrower line widths, as can be seen at the base of each diagonal peak.The GS reconstructions produced the most accurate reconstructions of the peak line widths and were able to maintain the noise level.The NAA and Glx cross peaks can be partially seen in each of the reconstructions; however, they are obscured by noise in the CS and TV reconstructions and have low amplitudes in the MaxEnt and GS 1 reconstructions.GS 2 was able to reconstruct the cross peaks to a reasonable degree in spite of the large water tail along F 1 at F 2 = 4.72 ppm that nearly obscured the lower NAA cross peak.The acquisition VOI used during the scan is shown as the white bolded line in the fully sampled image on the left and the black bolded line in the images on the right.The contour levels used for the fully sampled spectra were used for the 8× NUS and reconstructed spectra.As can be seen, the NAA/NAAG and Glx peaks are well resolved within the VOI of the fully sampled image with minimal spatial aliasing along Y and noticeably more aliasing along X, which was caused by the PSF of the EP-COSI pulse sequence.The peak amplitudes of the NAA/NAAG and Glx peaks are roughly equivalent over the VOI, with the exception of reduced peak amplitudes in the anterior-left and posterior-right voxels.The 8× NUS metabolite peak amplitudes are lower than in the fully sampled image, a trend also seen in Figures 2  and 4, and the increased spatial aliasing along Y is evident from the presence of NAAG and Glx peaks outside of the VOI.The CS and TV reconstructions show increased noise both inside and outside of the VOI, with CS failing to properly reconstruct many of the peaks in the posterior row of the VOI with the "choppy" line shapes that were seen in Figures 2 and 4   Table 3 shows the mean metabolite peak RMSEs calculated over the 4 × 4 voxel VOIs of the 4×, 6×, and 8× NUS reconstructions of the in vivo 4D EP-COSI brain scan shown in Figures 4 and 5.The RMSEs were calculated using the ppm locations shown in Table 2 and the lowest values calculated for each metabolite peak and NUS rate are highlighted.As can be seen, the majority of highlighted RMSEs are from the GS 2 reconstructions with a handful highlighted for the GS 1 and MaxEnt reconstructions.These results are similar to what was shown in the metabolite RMSE plots in Figure 2, which shows that GS 2 and GS 1 have most of the lowest values, followed by MaxEnt.Additionally, the values of the RMSEs for the large diagonal peaks of NAA and Cr are more tightly grouped than the smaller diagonal and cross peak RMSEs, and the high GS 1 RMSEs for Cho 3.2 ppm indicate it had trouble reconstructing the metabolite peak, as in Figure 2.
Table 3. Mean metabolite peak RMSEs (dB) for the 4×, 6×, and 8× NUS reconstructions calculated over the 4 × 4 VOI of the in vivo 4D EP-COSI brain scan shown in Figures 4  and 5.The lowest RMSE for each metabolite peak and NUS rate is highlighted.

Discussion
In this article, an algorithm to solve the GS reconstruction problem using the Split Bregman framework was derived that allows for arbitrary grouping patterns of transform coefficients and was applied to reconstructing retrospectively under-sampled 4D EP-COSI data.The results for GS reconstruction using overlapping and non-overlapping, equal sized groups were compared against the results from CS, TV, and MaxEnt reconstructions at NUS rates of 2×, 4×, 6×, 8×, and 10× on gray matter phantom and in vivo brain data.Figures 2, 4, and 5 show that spectra reconstructed by GS reconstruction using overlapping coefficient groups (GS 2 ) are the most qualitatively similar to the fully sampled spectra with respect to the noise floor, metabolite peak amplitudes, line widths, and spatial localization.The metabolite peak RMSEs calculated over the VOIs of the phantom and in vivo datasets in Figure 3 and Table 3 show that the GS 2 reconstructed spectra also have the lowest values for the majority of peaks evaluated, indicating that they are more accurately reconstructed than the other tested methods.The agreement between these results indicate that GS 2 is the most effective means of reconstructing NUS 4D EP-COSI phantom and in vivo data tested.The higher noise floor and dynamic range from water and fat contamination in the in vivo scans did not prevent GS 2 from successfully reconstructing the highly under-sampled data.
Compared with the other non-GS reconstruction methods, GS 1 produced high quality spectra and was able to reconstruct the metabolite peaks to a higher degree of fidelity; however, as indicated by the mean metabolite peak RMSEs in Figure 3 and Table 3, the Cho 3.2 ppm diagonal peak was not reconstructed by GS 1 as accurately as it was by the other reconstruction methods.This discrepancy was caused by the nature of the grouping method used for GS 1 , which used equally sized non-overlapping groups of (F 2 , F 1 ) = (8, 4) coefficients.The Cho 3.2 ppm diagonal peak was located at the intersection of multiple coefficient groups; therefore, the structural block of coefficients that comprised the peak were divided into multiple groups and were not correlated with one another.Instead of increasing the sparsity of the coefficients around the peak and reducing it at the peak, this division had the opposite effect and prevented the Cho 3.2 ppm peak from being properly reconstructed.In contrast, GS 2 was able to reconstruct the Cho 3.2 ppm diagonal peak to a higher degree of fidelity with the same sized

Figure 1 .
Figure 1.Illustration of the coefficient grouping support of a single data point for equal sized non-overlapping and 50% overlapping groups.

Figure 5
Figure5shows the spatial distribution of a 1.5 × 1.5 ppm region of a 2D COSY spectrum centered at F 2 = 2.5 ppm taken from each voxel of the fully sampled, 8× NUS, and reconstructed in vivo 4D EP-COSI brain scan shown in Figure4.The fully sampled data is superimposed onto a T1 weighted image from the middle of the EP-COSI slice that was acquired during the same study.Each spectral region contains the metabolite diagonals and cross peaks of NAA/NAAG and Glx in that region.The acquisition VOI used during the scan is shown as the white bolded line in the fully sampled image on the left and the black bolded line in the images on the right.The contour levels used for the fully sampled spectra were used for the 8× NUS and reconstructed spectra.As can be seen, the NAA/NAAG and Glx peaks are well resolved within the VOI of the fully sampled image with minimal spatial aliasing along Y and noticeably more aliasing along X, which was caused by the PSF of the EP-COSI pulse sequence.The peak amplitudes of the NAA/NAAG and Glx peaks are roughly equivalent over the VOI, with the exception of reduced peak amplitudes in the anterior-left and posterior-right voxels.The 8× NUS metabolite peak amplitudes are lower than in the fully sampled image, a trend also seen in Figures2 and 4, and the increased spatial aliasing along Y is evident from the presence of NAAG and Glx peaks outside of the VOI.The CS and TV reconstructions show increased noise both inside and outside of the VOI, with CS failing to properly reconstruct many of the peaks in the posterior row of the VOI with the "choppy" line shapes that were seen in Figures2 and 4. The MaxEnt reconstruction shows better denoising properties than the CS and TV reconstructions with most of the NAA/NAAG and Glx peak amplitudes and line widths restored within each voxel of the VOI.The NAA/NAAG and Glx peaks are properly localized within the VOIs of the MaxEnt and GS reconstructions to within the same levels of spatial aliasing seen in the fully sample image.The line widths and amplitudes of the GS 1 and GS 2 reconstructed peaks are the most qualitatively similar to the fully sampled scan.
Figure5shows the spatial distribution of a 1.5 × 1.5 ppm region of a 2D COSY spectrum centered at F 2 = 2.5 ppm taken from each voxel of the fully sampled, 8× NUS, and reconstructed in vivo 4D EP-COSI brain scan shown in Figure4.The fully sampled data is superimposed onto a T1 weighted image from the middle of the EP-COSI slice that was acquired during the same study.Each spectral region contains the metabolite diagonals and cross peaks of NAA/NAAG and Glx in that region.The acquisition VOI used during the scan is shown as the white bolded line in the fully sampled image on the left and the black bolded line in the images on the right.The contour levels used for the fully sampled spectra were used for the 8× NUS and reconstructed spectra.As can be seen, the NAA/NAAG and Glx peaks are well resolved within the VOI of the fully sampled image with minimal spatial aliasing along Y and noticeably more aliasing along X, which was caused by the PSF of the EP-COSI pulse sequence.The peak amplitudes of the NAA/NAAG and Glx peaks are roughly equivalent over the VOI, with the exception of reduced peak amplitudes in the anterior-left and posterior-right voxels.The 8× NUS metabolite peak amplitudes are lower than in the fully sampled image, a trend also seen in Figures2 and 4, and the increased spatial aliasing along Y is evident from the presence of NAAG and Glx peaks outside of the VOI.The CS and TV reconstructions show increased noise both inside and outside of the VOI, with CS failing to properly reconstruct many of the peaks in the posterior row of the VOI with the "choppy" line shapes that were seen in Figures2 and 4. The MaxEnt reconstruction shows better denoising properties than the CS and TV reconstructions with most of the NAA/NAAG and Glx peak amplitudes and line widths restored within each voxel of the VOI.The NAA/NAAG and Glx peaks are properly localized within the VOIs of the MaxEnt and GS reconstructions to within the same levels of spatial aliasing seen in the fully sample image.The line widths and amplitudes of the GS 1 and GS 2 reconstructed peaks are the most qualitatively similar to the fully sampled scan.

Table 1 .
Lagrange multiplier values used during reconstruction.The GS reconstructions use the CS parameter values scaled by the number of coefficients per group.The MaxEnt reconstruction does not use Lagrange multipliers.