Electric potential energy optimized 3D radial sampling trajectories for MRI

Comparing optimization strategies

Four ELECTRO trajectories were produced to demonstrate the effects of using the multi-stage optimization strategy and the effects of reducing the number of sphere sizes to be optimized. Computation time for the non-multi-stage optimization with all sphere sizes limited this comparison to trajectories with only 2500 readouts. Each optimization was run for 30,000 iterations. First, these four trajectories were analyzed by evaluating the full objective function \(\:G\left(S\right)\) for each iteration of the optimization, even if the full \(\:G\left(S\right)\) was not used for that stage or optimization strategy (Fig. 4A). The multi-stage optimization strategies initially ignore most of the terms in \(\:G\left(S\right)\), resulting in a delayed drop in \(\:G\left(S\right)\) compared to the non-multi-stage optimization strategies. However, the multi-stage optimization strategies both converged to lower values of \(\:G\left(S\right)\).

In the multi-stage optimization, each point experiences forces from at most \(\left( {2m – 2} \right)\) other points, where \(\:m\) is the largest sphere size in the current stage. Without multi-stage optimization, each point experiences forces from every other point. This makes it difficult for a point to reach its optimal position because it would need to push through many other points and the forces pushing a point to its optimal position get weaker as the point gets closer. The fewer force interactions in the earlier stages of the multi-stage approach allows points to move more easily to their approximate optimal positions. As the optimization proceeds to the later stages these points stay approximately in the same position, while obtaining low energy configurations for both the small and large spheres. Figure S1 in the Supplementary Material contains videos showing the point movements for each of the four optimization strategies. The multi-stage optimization strategy is also more computationally efficient because it has fewer forces to compute in all stages except for the last, where the number of forces to be calculated is the same as a non-multi-stage optimization.

For the final solution of each type of optimization, \(\:{H}_{m}\left(S\right)\) was evaluated and compared in Fig. 4B and C. Figure 4B shows the range of sphere sizes from 2 to 64, while Fig. 4C shows the range of sphere sizes from 2437 to 2500. To enable visual comparison, \(\:{H}_{m}\left(S\right)\) was normalized by dividing by the number of spheres \(\:\left(N-m+1\right)\) and by dividing by the number of interactions in a size \(\:m\) sphere, \(\:\left({m}^{2}-m\right)/2\). Figure 4B shows that the energies for the smaller sphere sizes were indeed lower for the multi-stage optimizations. Figure 4B and C show that optimizing for only the sphere sizes from Narayana’s cows sequence resulted in lower energies for the smaller sphere sizes, but higher energies for the larger sphere sizes. This trade-off may be favorable for cine and real-time imaging where it is more important for smaller bundles of readouts to be well spread out. For the non-multi-stage optimizations, optimizing for only Narayana’s cows sphere sizes resulted in a lower \(\:G\left(S\right)\) (Fig. 4A), likely from the lower energies of the smaller sphere sizes. On the other hand, for the multi-stage optimizations, optimizing for all sphere sizes resulted in convergence to the lowest value of \(\:G\left(S\right)\) out of the four methods. However, the multi-stage optimization with Narayana’s cows sphere sizes requires much less computation than that with all sphere sizes. The former reached the final stage in 7 s, at iteration 1824, while the latter reached the final stage in 330 s, at iteration 3956. To keep computation times feasible for the larger ELECTRO trajectories, we proceed with the rest of our analysis using the multi-stage optimization with Narayana’s cows sphere sizes.

ELECTRO trajectories with increasing N

ELECTRO trajectories using the multi-stage optimization with Narayana’s cows sphere sizes were produced for 10,000, 20,000, 40,000, and 80,000 readouts. The optimizations were set to run for 10,000 iterations, except for the ELECTRO trajectory with 80,000 readouts, which was set to run for 20,000 iterations because it had more stages and it converged more slowly. The global NMNA of the first \(\:N\) readouts and the last \(\:N\) readouts were evaluated for each trajectory and the values are shown in Fig. 5. All ELECTRO trajectories appear to have a drop in NMNA towards the highest sphere sizes, suggesting that under the current scheme the weights \(\:{\alpha\:}_{m}\) for the highest \(\:m\) are too small or there should be more sphere sizes included in \(\:G\left(S\right)\) in the highest \(\:m\) range. This drop in NMNA correlates with the higher energies for the larger sphere sizes in Fig. 4C. To avoid this potential sampling bias, these trajectories can simply be truncated given the NMNA is fairly consistent for most \(\:m\), e.g., take only the first 60,000 readouts from the ELECTRO trajectory optimized for 80,000 readouts. The 80,000 readout ELECTRO trajectory took 7.8 h to reach the final stage (at iteration 10,026) and 76.3 h to complete all 20,000 iterations.

ELECTRO trajectories were produced for 10,000, 20,000, 40,000, and 80,000 readouts. The global NMNA was evaluated for the first \(\:N\) and last \(\:N\) readouts.

Regional point distribution analysis

The regional NMNA was evaluated for five spherical caps on the unit sphere. Since the spatial distortion from the inverse cylindrical equal-area projection varies as a function of the polar angle \(\:\theta\:\) but not the azimuthal angle \(\:\phi\:\) (Fig. 2A–C), four spherical caps were positioned at \(\:\theta\:=0^\circ\:,30^\circ\:,60^\circ\:,90^\circ\:\), each with \(\:\phi\:=0^\circ\:\) and cap angle \(\:\beta\:=15^\circ\:\). The fifth spherical cap was the whole sphere. The NMNA was evaluated as a function of the first \(\:N\) readouts of the ELECTRO, supergolden, plastic, 2D bit-reversed, and Halton trajectories. This was evaluated up to 80,000 spokes for every trajectory except for the 2D bit-reversed trajectory, which only had 40,000 spokes so that the same trajectory can be used in the next analysis (i.e., Fig. 7). The two parameters used for the 2D bit-reversed trajectory were 20 and 2000. Figure 6A– E compare the regional NMNA between the different spherical caps for each trajectory. Figure 6F compares the global NMNA between the different trajectories. The supergolden and plastic trajectories have the largest variation in regional NMNA between the four regions of the sphere, and large changes in regional NMNA as a function of \(\:N\). The NMNA of the 2D bit-reversed trajectory generally increases as a function of \(\:N\), until near uniformity when all readouts are included. The 2D bit-reversed, Halton, and ELECTRO trajectories have relatively little variation in regional NMNA between the four regions. The Halton and ELECTRO trajectories have relatively consistent values of NMNA for all \(\:N\) evaluated, though the values for the ELECTRO trajectory are higher. Like the global NMNA, the regional NMNA has drops at large sphere sizes for the ELECTRO trajectory. The supergolden trajectory at 40,000 readouts has greater clustering in the region at \(\:\theta\:=0^\circ\:\) than the Halton trajectory. This can be seen both visually (Fig. 1A) and in their regional NMNA values (1.28 for supergolden and 1.33 for Halton), but the global NMNA of the supergolden trajectory remains higher than that of the Halton trajectory (1.37 for supergolden and 1.24 for Halton). This highlights the importance of regional analysis and how local information can be lost due to global averaging.

Regional point distribution analysis of the trajectories using NMNA. (A–E ) Regional NMNA of five center-out radial trajectories. Four regions at different polar angles were evaluated. (F) The NMNA averaged over the whole sphere for each of the five trajectories.

Uniformity of the smaller spheres

For each of the sphere sizes from 2 to 1000, the global NMNA was calculated for each of the \(\:\left(N-m+1\right)\) spheres. The average and standard deviation of the values were evaluated for each sphere size, and the results are shown in Fig. 7A and B respectively. Computation time limits this analysis to the lower range of sphere sizes, though this range includes realistic temporal resolutions for cine and real-time MRI applications. Trajectories with 40,000 readouts were included in this comparison. The ELECTRO trajectory with 80,000 readouts was truncated by taking the first 40,000 readouts. To measure the flatness of the curves, \(\:\sigma\:\), in Fig. 7A the standard deviation of the values in each curve was calculated. In the range of sphere sizes from 2 to 1000, the ELECTRO trajectory has stable values for the average NMNA (\(\:\sigma\:\)=0.005) that are consistently higher than those from the Halton and 2D bit-reversed trajectories (except for the sphere size of 2, for which the Halton trajectory has a higher average NMNA). The supergolden and plastic trajectories have average NMNA values that dramatically oscillate above and below the values of the ELECTRO trajectory (\(\:\sigma\:\)=0.090 and \(\:\sigma\:\)=0.070 respectively). The 2D bit-reversed and Halton curves have moderate flatness values at \(\:\sigma\:\)=0.014 and \(\:\sigma\:\)=0.019 respectively. The standard deviation of the NMNA (Fig. 7B) is generally lowest for the supergolden and plastic trajectories, followed by the ELECTRO trajectory and then the Halton and 2D bit-reversed trajectories. Compared to the supergolden trajectory, the ELECTRO trajectory had significantly higher (Welch’s t-test, \(\:p) average NMNA values for 771 sphere sizes, while the supergolden trajectory had significantly higher values for 227 sphere sizes. Compared to the plastic trajectory, the ELECTRO trajectory had significantly higher average NMNA values for 336 sphere sizes, while the plastic trajectory had significantly higher values for 662 sphere sizes. Between the supergolden and plastic trajectories, the supergolden trajectory had significantly higher values for 222 sphere sizes while the plastic trajectory had significantly higher values for 777 sphere sizes.

Comparison of NMNA values for the lower range of sphere sizes containing realistic temporal resolutions for cine and real-time MRI applications. Five center-out radial trajectories with 40,000 readouts are compared. (A) Average NMNA and (B) standard deviation of NMNA.

Point spread function analysis

PSFs were calculated for the 40,000 spoke 3D center-out radial trajectories. Radial spokes with instantaneous gradient ramping were used (i.e., no uneven sampling on the spoke). The density compensation function was derived geometrically, assuming all points in a spherical shell occupied the same volume (no compensation for uneven angular spacing). The same density compensation function was used for each trajectory. The non-uniform fast Fourier transform⁴³ from the SigPy library ( was used.

The three orthogonal center planes of the PSFs (Fig. 8) show the different aliasing artifacts from each trajectory and their directional dependence. Most noticeably, the PSF of the supergolden trajectory has disk-like artifacts in the mid z-plane, likely caused by the spiral sample clustering around the z-axis (Fig. 1A). The supergolden, plastic, 2D bit-reversed, and Halton trajectories all have one type of aliasing pattern in their mid x- and y-planes and a different type of aliasing pattern in their mid z-planes. This is likely a result of using the spherical polar and azimuthal angles in their derivation, which results in a symmetry about the z-axis. Of these four trajectories, the differences in aliasing pattern can be seen within the 1 FOV radius for the supergolden, plastic, and Halton trajectories, while the differences for the 2D bit-reversed trajectory are mainly outside of the 1 FOV radius. The ELECTRO trajectory is the only trajectory that has PSF mid planes that are indistinguishable from each other. The PSFs of the Halton and ELECTRO trajectories have mostly streak-like aliasing, correlating with the random-like pattern of their readout directions, while the PSFs of the supergolden and plastic trajectories have blade-like aliasing, which may correspond to their readouts having patterns that are regularly repeated azimuthally. Having incoherent aliasing artifact is favorable for L1-regularized image reconstructions such as compressed sensing⁶ and total-variation minimization⁴⁴, while not having any directional dependence simplifies the assumptions to be made in these reconstructions. Any exponentially increasing sequence could be tried instead (e.g. Fibonacci or Padovan), where the base of the exponent would be chosen to control the trade-off between computation time and spacing between the numbers in the sequence. Narayana’s cows sequence was chosen because its base sits between those of the Fibonacci and Padovan sequences.

Point spread function comparison between five different 3D center-out radial trajectories. A narrow display window of [0, 0.004] was used. Images show the three orthogonal center planes at 2x FOV. The center planes of the ELECTRO PSF are indistinguishable from each other.

Evaluation of experimental factors on image quality

As described in the Methods section, numerical simulations of image acquisition and reconstruction were performed using the different trajectories to evaluate the effect of realistic experimental factors on image quality. Image quality was quantified based on the root mean squared error (RMSE) between the ground truth image (digital thoracic phantom) and the reconstructed images.

Of the factors tested, an uncorrected gradient delay produced the bigger change in RMSE across all trajectories. For example, for a sphere size of 40,000, the RMSE for all trajectories in the absence of gradient delay was less than 1% of the maximum signal of the ground truth image. Under conditions of experimentally realistic and uncorrected gradient delays, this error rose over six to eight-fold; however, following retrospective gradient delay correction, RMSE values returned to their baseline values. Similar behaviour was found for a sphere size of 20,000, with RMSE values rising from 1 to 2% to 5–6%, but again returning to baseline after correction. Finally, for a sphere size of 10,000, baseline RMSE values ranged between 5 and 8%, rising to 7.5–9.5%, and returned to baseline after correction. These findings demonstrate that, although RMSE values varied across sphere sizes and between trajectories, gradient delay errors after retrospective correction had a small influence on image quality.