Perturbed spiral real‐time phase‐contrast MR with compressive sensing reconstruction for assessment of flow in children

we implemented a golden‐angle spiral phase contrast sequence. A commonly used uniform density spiral and a new ‘perturbed’ spiral that produces more incoherent aliases were assessed. The aim was to ascertain whether greater incoherence enabled more accurate Compressive Sensing reconstruction and superior measurement of flow and velocity.


| INTRODUCTION
Phase-contrast magnetic resonance (PCMR) is a proven method of measuring blood flow in the clinical environment. 1,2 Such acquisitions are usually cardiac-gated, enabling collection of high spatio-temporal resolution data. In children, cardiac gating is often combined with signal averaging to allow free-breathing acquisition with minimal respiratory artifact. However, this is a time-consuming approach that significantly prolongs scan time. Accelerated breath-hold PCMR is an alternative, 3 but some children find even short breath-holds difficult. Thus, a rapid free-breathing approach is desirable.
One solution is real-time imaging and several real-time PCMR sequences have been described. 4 Most rely on a combination of parallel imaging (i.e., SENSE 5 or GRAPPA 6 ), efficient k-space filling (i.e., spiral 7 or EPI 8 ) and temporal undersampling (i.e., UNFOLD 9 or k-t BLAST 4 ) to reduce acquisition time. Unfortunately, these methods fail at very high acceleration rates, limiting the achievable spatio-temporal resolution.
It has been shown that compressive sensing (CS) can reconstruct high quality images from heavily undersampled k-space data. 10 However, a prerequisite of CS is that aliasing is incoherent. One way this can be achieved is by combining non-Cartesian trajectories (i.e., radial 11 or spiral 12 ) with golden-angle rotations. Spiral trajectories are of particular interest for PCMR due to short TEs and highly efficient filling of k-space. Golden-angle spiral imaging with CS reconstruction has been shown to be sufficient for real-time cine data. 11,12 However, we have shown in a pilot study 13 that CS reconstruction can cause temporal blurring of PCMR data that could result in underestimation of clinically important metrics. A possible solution may be modification of k-space sampling to produce more incoherent aliases, which should improve data conditioning for the CS reconstruction.
The general aim of this study was to implement a perturbed spiral PCMR acquisition. The specific aims were to (1) find the optimum perturbed spiral trajectories for CS reconstruction of PCMR data using point spread functions (PSFs) and in silico simulations, and (2) validate the developed technique in vivo.

| Perturbed spiral design
All trajectories were designed using a modification of the method described by Pipe. 14 The starting point for the new trajectory design was a uniform density spiral sequence with 36 evenly spaced interleaves required to completely fill k-space. To perform 1-sided velocity encoding, each readout was acquired twice (velocity encoding and compensation acquisitions). The initial reference trajectory parameters were set to FOV = 450 × 450 mm, voxels = 1.76 × 1.76 × 6.0 mm, TR/TE = 8.54/1.93 ms, and velocity encoding = 200 cm/s. These uniform density readouts were continuously rotated by the golden angle (~222º), resulting in a golden-angle spiral acquisition (GAS uniform , Figure 1). To achieve the desired temporal resolution (< 30 ms), the GAS uniform sequence was 18-times undersampled (2 interleaves per frame). Spiral aliases are observed as concentric rings 15,16 in their PSFs (Figure 2), and the position of these aliases depends on the radial undersampling of the spiral trajectory. Thus, it is possible to increase the incoherence of these aliases by modifying the radial undersampling. To achieve this, we developed a trajectory design algorithm that sinusoidally varied the radial acceleration (α r ) as a function of normalized distance from k-space center (r). This can be described formulaically as follows (and graphically in Supporting Information Figure S1): where α 0 is the maximum prescribed radial acceleration, β is the number of oscillations in r ∈ 0.5, 0 − 0.5 between the center and edge of k-space, and ϕ 0 is an additional phase offset. The parameters c 0 and c 1 (c 0 , c 1 [0, 1] ∩ c 0 ≤ c 1 ) divide k-space into 3 sections: (1) a central section (r < c 0 ) with 2-times oversampling (α r = 0.5) and no oscillations; (2) a transition section (r ≥ c 0 ∩ r < c 1 ) with linearly increasing radial acceleration and oscillation amplitude; and (3) an outer section (r ≥ c 1 ) with the maximum acceleration and oscillation amplitude ( r = 0.5 0 for = 0 or r = 0 − 0.5 for > 0). In this study, ϕ 0 was varied between 0º and 360º with 10º increments. This produced a set of 36 perturbed trajectories for any given set of α 0 , c 0 , c 1 , and β with indices (j) between 0 and 35. The exact perturbed trajectory used for any given readout (with an index [i] between 0 and infinity) was chosen using the following series index translation: As with GAS uniform , consecutive readouts were rotated by the golden angle, resulting in a perturbed golden-angle spiral (GAS perturbed ) sequence. It should be noted that even though perturbed trajectories were reused (every 36 readouts), they were always in different k-space position due to the golden-angle rotation. Varying the values of α 0 , c 0 , c 1 , and β allows different perturbations to be designed, and 2 examples of perturbed spirals trajectories are shown in Figure 1 (1)
F I G U R E 1 Trajectory visualization. Shown are the uniform golden-angle spiral acquisition (GAS uniform ), perturbed golden-angle spiral acquisition (GAS perturbed ), and an additional example presenting possible perturbations induced to the spiral trajectory. The k x -k y positions of composite trajectories are presented for 3 consecutive imaging frames. Also plotted are variations in the radial distance (r) of individual samples with their coordinates (k x and k y ) (see Supporting Information Figure S2 for accompanying k x and k y versus time plots).
One issue with perturbed spirals is that the number of readout samples required to reach the edge of k-space may vary. Consequently, we restricted the number of readout samples to the number in the shortest readout for a given set of parameters (α 0 , c 0 , c 1 , and β). This was done to ensure a constant TR throughout the acquisition and resulted in some readouts terminating before reaching k-space edge ( Figure 1). The GAS perturbed can also have different sampling acceleration (SA = N i /N t , in which N i represents the image pixels and N t the total trajectory samples) when compared with GAS uniform . Therefore, the number of readouts combined into an imaging frame had to be adjusted for each set of parameters, to ensure comparable acquisition times and SA. This was incorporated into test scripts using the GAS uniform trajectory as a reference point (N uniform t = 5862, SA uniform = ∼ 11.2). Due to the impact of multiple excitation times and inability to split readout samples, the maximum number of readouts was set to 3 and the minimum sampling acceleration to 93% of SA uniform .

| Point spread function evaluation
The GAS petrubed should produce more incoherent aliases than GAS uniform . However, the level of incoherency will depend on α 0 , c 0 , c 1 , and β, which must be optimized. The PSF is used commonly to ascertain the features of a given sampling pattern and optimize trajectory parameters. In this study, we used the amount of energy leakage to the PSF side lobes 17 as an incoherence metric (× point-bypoint multiplication): The calculations were done for all combinations of α 0 , c 0 , c 1 , and β given in Table 1. The trajectory (GAS perturbed , Figure 1) with the lowest average energy leakage (Figure 3,  indicating higher incoherence of artifacts) was selected for further tests.

| Reconstruction
Compressive sensing solves a set of nonlinear equations (representing the imaging process) through minimization of a cost function. The cost function used in this study is The first term enforces data consistency, where ρ are the image data, E is the encoding matrix (the multicoil nonuniform Fourier transform operator), and y are the acquired k-space data. The additional terms enforce sparse results through L1 norm regularization. In this study, finite difference operators (or total variation) were applied in space and time (TV [spatial], TV t [temporal]) as the sparse transforms.
The optimization was performed using a nonlinear conjugate gradient algorithm. For fast reconstruction, the described CS algorithm was implemented on an external graphics processing unit-equipped computer (Tesla K40c; NVIDIA, Santa Clara, CA) with online communication to the native reconstructor. 18 Acquired data were reconstructed in blocks of 90 frames with coil-sensitivity maps calculated from the time-averaged (flow-compensated) data from each block. These blocks were overlapped by 3 frames on each side to counter potential jump discontinuities. If there was no adjacent block (the start and end of the first and last block), the expansion was achieved by mirroring frames.
Gridding of non-Cartesian samples onto a rectilinear grid 19,20 requires information about the density of the samples. This is not provided by the described trajectory generation algorithm. Consequently, density compensation coefficients were calculated using the method described in Bydder, 21 which was chosen because it required no finetuning. The method uses a linear optimization process to find density distributions for a set of arbitrary trajectory points. The sample's density calculation was implemented as the first step of the graphics processing unit-based MRI reconstruction process.

| In silico model
The GAS uniform and chosen GAS perturbed sampling patterns were first evaluated in an in silico model, enabling comparison with ground-truth data. The in silico model (Supporting Information Figure S3) consisted of a cross section of the body through the ascending aorta. Aortic velocity, distension, signal intensity, and in-plane motion were modeled on data extracted from a high temporal resolution (~20 ms) Cartesian-gated PCMR data set acquired in a healthy volunteer. Respiratory motion was modeled using a function consisting of expansion (inhalation), a brief pause, and contraction (exhalation). The respiratory rate was set to 10 breaths per minute with the maximum body contraction to 98% of its original size. Simulated MR data were created by generating the image data at each of the readout time positions. The readout time was set to the GAS uniform TR of 8.54 ms. The complex phase component was generated through interpolation of the flow curve. The produced models were scaled with synthetically generated coil sensitivity maps (12 coils as in the in vivo experiment). This was then transformed into k-space and gridded onto the tested trajectory.

| In vivo study
The GAS uniform and chosen GAS perturbed sampling patterns were also assessed in a patient population consisting of 20 children referred for cardiac clinical MR (7 females and 13 males; age range: 6-16, median: 12.5 years). The only exclusion criterion was irregular rhythm (i.e., multiple ectopic beats or atrial fibrillation). The National Research Ethics Committee approved the study (06/Q0508/124), and written consent was obtained from all patients or legal guardians of children.
All imaging was performed on a 1.5T MR scanner (Magnetom Avanto; Siemens Medical Solutions, Erlangen, Germany) using the standard 2 spine coils and 1 body matrix coil setup (giving a total of 12 coil elements) used in all children at our institution. A vector electrocardiographic system was used for cardiac gating in the reference Cartesian-gated PCMR acquisitions. The imaging plane for aortic flow assessment was located in the ascending aorta just above the sino-tubular junction. The reference standard flow acquisition was a conventional free-breathing Cartesian retrospectively gated PCMR sequence with the following parameters: FOV = 350 × 262 mm, voxels = 1.82 × 1.82 × 6.0 mm, TR/TE = 4.4/1.9 ms, flip angle = 30º, velocity encoding = 200 cm/s, number of signal averages = 2, GRAPPA = 2, and temporal resolution = 18.5 ms.
Both real-time GAS PCMR acquisitions were set to FOV = 450 × 450 mm, voxels = 1.76 × 1.76 × 6.0 mm, flip angle = 20º, and velocity encoding = 200 cm/s. The bandwidth per pixel was optimized separately to minimize the trajectory errors' impact on image quality. This was done empirically based on a single in vivo case. These adjustments affected the length of a readout. The final TR/TE values were 6.7/1.9 ms for the GAS perturbed and 7.5/1.9 ms for GAS uniform acquisitions, resulting in a temporal resolution of approximately 26.6 ms (about 2.4 seconds for 90 frames) and approximately 29.9 ms (about 2.7 seconds for 90 frames), respectively. A relatively large FOV was chosen to ensure that even in older children there was no signal from outside the FOV.
The regularization parameters (λ 1 , λ 2 ) were selected empirically as a trade-off between image quality and minimization of spatial and temporal blurring.

| Flow quantification
The aorta was segmented on the magnitude images using a semi-automatic method based on the optical flow registration 23 with manual operator correction using in-house plugins for Osirix software (OsiriX Foundation, Switzerland). 24 The resultant regions of interest (ROIs) were transferred to the phase images to produce flow and velocity curves. Maximum velocity was taken as the peak of the velocity curve. Stroke volume was calculated by integrating the resultant flow curve over a single r-r interval. As multiple heartbeats are evaluated with real-time PCMR, the stroke volume and peak velocity are averaged across all r-r intervals.

| Image quality
All quantitative analyses were carried out using in-house plug-ins for OsiriX software, version 9.0. 24 True quantification of SNR and velocity-to-noise ratio (VNR) in images acquired with non-Cartesian trajectories is nontrivial in the clinical environment due to the uneven distribution of noise. 25,26 Therefore, in this study, estimated SNR and VNR were calculated as previously described. 27 In summary, an ROI was drawn in stationary tissue, and estimated noise was calculated as the average SD of the pixel intensity (σ s ) or velocity (σ v ) through time. Final estimates of SNR and VNR were made by dividing the mean signal intensity from an ROI drawn at peak systole by σ s and σ v , respectively.
Quantitative edge sharpness was calculated in peak systole by measuring the average maximum gradient of the normalized pixel intensities across the aortic wall. The image data were resampled onto evenly spaced perpendicular lines crossing the vessel border (marked with the ROIs used to extract the velocity data). Lanczos resampling 28,29 was used with a 0.5-mm step between samples on the lines with a distance of 20 mm. Furthermore, the smooth noise robust differentiation 30 was applied to extract the maximum gradient on the projections.
In real-time data, the SNR, VNR, and edge sharpness measurements were performed in all peak systole frames, and the averaged values were used in comparisons.
Subjective image quality scoring for the GAS uniform and GAS perturbed sequences was done by 2 independent, experienced observers (V.M. and D.K.) who were presented with the magnitude data for each patient in a blinded, randomized manner. The Cartesian data were not included, as they were obvious to the observer, which risked bias. The images were graded on a Likert scale (

| Statistical analysis
All statistical analyses were performed using R software (R Foundation for Statistical Computing, Vienna, Austria) and a p-value of less than .05 indicated a significant difference. All of the results are expressed as mean ± SD. Differences among the 3 imaging techniques were assessed using the 1-way repeated-measures analysis of variance. The imaging techniques were treated as the repeated measure factor. Significant results were further investigated with post hoc pairwise comparison using the Tukey method.
Qualitative image scores were compared using 1-way analysis of variance, as previous work has shown that there is a lower chance of type II errors compared with nonparametric tests for Likert scale data. 31 It is therefore more likely to detect differences among the techniques. The scores provided by the observers were treated as individual factor measures.

| Trajectory optimization
Energy leakage results for the range of GAS perturbed trajectories (Table 1) are presented in graphical form in Figure 3. The optimal trajectory was found for the following parameters: c 0 = 0.2, c 1 = 0.9, 0 = 2.5, and = 0.33. This corresponds to a trajectory that is oversampled at the center with a low frequency oscillation in radial acceleration that slowly increases in amplitude (Figure 1 and Supporting Information Figure S1). Consequently, each trajectory covers a different portion of k-space and contributes unique information. These temporal sampling density distribution changes are visible in the PSFs as changes in distribution of the side lobes between frames (Supporting Information Figure S4 and Supporting Information Video S1). The PSFs of the GAS uniform , the optimized GAS perturbed trajectory, and 1 of the nonoptimized GAS perturbed trajectories are shown in Figure 2 and Supporting Information Video S1. This visually demonstrates the greater incoherence provided by GAS perturbed trajectories compared with GAS uniform and the importance of optimizing the perturbation to increase incoherence.

| In silico model
Magnitude and phase data reconstructed with optimal regularization parameters (the lowest NRMSE) for GAS perturbed (NRMSE = 0.03, 1 = 7.5e −5 , 2 = 7.5e −5 ) and GAS uniform (NRMSE = 0.06, 1 = 4.0e −4 , 2 = 5.0e −6 ) sequences are shown in Figure 4 along with the velocity curves extracted from these data sets. There is significant blurring of the GAS uniform velocity curve compared with the ground truth, resulting in underestimation of the peak velocity. However, there is minimal blurring of the GAS perturbed velocity curve with good agreement of peak velocity.
The effect of regularization parameters on NRMSE for GAS perturbed and GAS uniform is shown in Figure 5. For GAS uniform , increasing temporal regularization (λ 1 ) reduces NRMSE, whereas increasing spatial regularization (λ 2 ) has a small detrimental effect. This pattern can also be appreciated in the extracted velocity curves at different levels of regularization. At low levels of temporal regularization, curves exhibit artifacts due to unresolved coherent aliasing ( Figure 5). Increasing temporal regularization removes these artifacts, but results in temporal blurring. Changing spatial regularization has minimal effect on the shape of the velocity curve for GAS uniform . For GAS perturbed , all NRMSE results are lower compared with GAS uniform , with less variation as a function of regularization ( Figure 5). The GAS perturbed velocity curves with different regularization levels are very similar, with none exhibiting artifacts due to coherent aliasing. The highest NRMSEs were found with high levels of temporal regularization, which caused temporal blurring. Conversely, the lowest NRMSEs were found with high levels of spatial regularization and low levels of temporal regularization.

| Feasibility
The PCMR data were successfully acquired in all 20 children during free breathing. Reconstruction time for each real-time block of 90 frames from GAS perturbed and GAS uniform was about 52 seconds, and all 270 frames were available for viewing on the scanner in about 160 seconds. The regularization levels were optimized separately for GAS uniform and GAS perturbed . The parameters were set to 1 = 5.0e −5 and 2 = 1.0e −5 for both reconstructions based on a visual assessment of results from a single subject. The mean heart rate of the study population was 81 ± 12 (range: 60-108) beats per minute. The reference standard Cartesian free-breathing gated acquisition required 84 heart beats to complete, resulting in 63 ± 10 seconds (range: 47-84 seconds) acquisition time.

| In vivo flow quantification
Examples of velocity and flow curves generated by the Cartesian, GAS uniform , and GAS perturbed acquisitions from the same subject are shown in Figure 6. As in the in silico results, there is substantial blurring of the velocity curves derived from GAS uniform data. This resulted in significantly lower (P < .001) peak velocity measured from the GAS uniform data (68.7 ± 18.4 cm/s) compared with the Cartesian reference (72.4 ± 18.0 cm/s). This bias was also associated with relatively broad limits of agreement (bias: −3.7, limits: −10.4 to 3.0 cm/s; Figure 7). There was much less blurring of the velocity curve derived from GAS perturbed data (72.3 ± 18.6 cm/s) with no significant difference (P = .98) in the peak velocity compared with the Cartesian reference. In addition, there were narrower limits of agreement (bias: −0.1, limits: −4.4 to 4.1 cm/s; Figure 7).
Aortic stroke volume (Figure 7) results showed no statistical difference between Cartesian (73.2 ± 23.7 mL) and both GAS perturbed (71.4 ± 23.4 mL, P = .19) and GAS uniform (74.5 ± 26.0 mL, P = .40) acquisitions. The GAS perturbed F I G U R E 4 In silico results. The imaging results of the in silico reconstruction present a frame at a peak velocity and a temporal cross section through the simulated ascending and descending aorta, as marked. The bottom plots compare the extracted mean flow aortic velocities for the regularization parameters giving the lowest normalized RMS error (NRMSE) | 2085 KOWALIK et AL.

F I G U R E 5
In silico regularization-level optimization. The NRMSE results for GAS perturbed (top left) and GAS uniform (top right) are presented. The effects of regularization are shown with plots of 4 flow curves: the least accurate, the flow curve produced with the combination of regularization parameters that resulted in the worst NRMSE; the lowest spatial regularization ( 2 = 0.25e −5 ), the best result while varying only the temporal regularization (λ 1 ); the lowest temporal regularization ( 1 = 2.5e −5 ), the best result while varying only the spatial regularization (λ 2 ); and the most accurate, the combination of regularization parameters that produced the best NRMSE result. These were plotted against the flow curve extracted from the fully sampled spiral trajectory acquisition produced a small insignificant underestimation of −1.8 mL (limits of agreement: −9.4 to 5.8 mL), whereas the GAS uniform acquisition produced a small, insignificant overestimation of 1.3 mL (limits of agreement: −8.8 to 11.4 mL).

| In vivo image quality
Representative imaging results are shown in Figure 8 and Supporting Information Videos S2-S4. No significant difference (P = .28) was found in the subjective image scoring between GAS perturbed (3.5 ± 0.6) and GAS uniform (3.3 ± 0.7) real-time imaging.

| DISCUSSION
The main findings of the study were as follows: (1) CS reconstruction applied directly to GAS uniform PCMR data results in clinically significant blurring of velocity data; (2) the proposed GAS perturbed trajectory produces better conditioned data for CS reconstruction, as indicated with lower PSF energy leakage; and (3) this resulted in more accurate measurement of peak aortic velocity in silico and in vivo.
The benefit of combining CS with parallel imaging is that it allows much higher acceleration factors compared with parallel imaging alone (i.e., SENSE 5 ) or temporal encoding (i.e., UNFOLD 9 or k-t BLAST 4 ). However, the performance of CS is dependent on how well the algorithm's requirements are met, with the most difficult being incoherent aliasing. This is because most MR systems cannot produce the sharp, nonsmooth changes in gradient moments needed for true random sampling. Several studies have been undertaken to identify MR data sampling patterns that are conducive to CS. 17,[32][33][34][35][36][37][38] It should be noted that gradient waveform design optimization for more complex trajectories has been done. However, these can be time-consuming algorithms with runtimes measured in minutes. 39 A commonly taken approach is to combine non-Cartesian sampling with golden-angle rotations. Examples include both radial 11 and spiral 12 golden-angle acquisitions that use temporal total variation L1 regularization to remove imaging artifacts without introducing clinically important temporal blurring. However, real-time PCMR is a more challenging problem, as it requires higher acceleration factors (18 times in this study) due to the need for velocity-encoded and compensated readouts. In this study, we have shown that phase data are more susceptible to temporal regularization and that application of CS reconstruction to GAS uniform data resulted in significant blurring of velocity curves.
One possible solution is to perturb spiral trajectories to generate more incoherent artifacts, as suggested by Lustig et al. 40,41 However, to our knowledge, this approach has not been applied in clinical studies. This may be because in most applications adequate data quality can be obtained using more conventional methods. This is not the case for real-time PCMR; therefore, we designed a family of perturbed spiral trajectories that could be implemented on a standard clinical scanner. Using PSF energy leakage, we showed that spirals that were oversampled at the center with a low frequency oscillation in radial acceleration that slowly increased in amplitude had the greatest incoherence. When this GAS perturbed sampling pattern was tested in silico, it was demonstrated to be better conditioned for CS reconstruction than GAS uniform . In particular, a high level of temporal regularization was required to remove coherent aliasing artifacts in GAS uniform data, resulting in significant temporal blurring. On the other hand, GAS perturbed data did not exhibit coherent aliasing artefacts, even at low levels of regularization. In addition, the greater spatial incoherence of the GAS perturbed aliases enabled spatial regularization to be used to further improve the image quality.
In the clinical study, both the GAS perturbed and GAS uniform acquisitions produced good-quality magnitude images with little residual artifacts. This is in keeping with previous studies that combined spiral imaging with CS reconstruction 12 and should be expected, as regularization was optimized for image quality. However, velocity and flow curves extracted from GAS uniform data were blurred in time compared with the reference standard Cartesian-gated data. This did not affect the quantification of stroke volume, as temporal blurring has a minimal effect on the integral of the flow curve. Nevertheless, it did result in significant underestimation of the peak mean velocities, which limits the clinical utility of GAS uniform realtime PCMR. In contrast, the GAS perturbed acquisition produced significantly less blurred flow and velocity curves. This resulted in good agreement with the Cartesian reference for both stroke volume and peak velocity quantifications. We believe that this demonstrates that the greater incoherence provided by perturbed spiral trajectories allows more accurate reconstruction of real-time PCMR data. This in turn widens the clinical utility of this technique in children with heart disease.
One well-recognized problem with spiral imaging is reduction in image quality due to trajectory errors. This can be mitigated by keeping readout lengths short and optimizing sampling bandwidth per pixel, as done in this study. There are also spiral deblurring algorithms available, although they were not used in this study as they would increase reconstruction time. A concern with perturbed spirals is that they might result in even greater trajectory errors. However, we saw no difference in qualitative image scores or edge sharpness between GAS perturbed and GAS uniform data. This suggests that trajectory errors were not seriously exacerbated by perturbing the spiral trajectory. Nevertheless, edge sharpness was slightly lower in the spiral acquisitions compared with the Cartesian data. This might be due to trajectory errors, the effect of spatial regularization, and the slightly lower acquired spatial resolution of the spiral acquisitions.
Another problem with spiral real-time approach is the lower SNR compared with the reference standard Cartesian acquisition. This is to be expected due to heavy undersampling of the real-time data and use of 2 signal averages for the Cartesian data. However, this does not appear to affect