Research Article  Open Access
Andrew Thomas Curtis, Christopher Kumar Anand, "Random Volumetric MRI Trajectories via Genetic Algorithms", International Journal of Biomedical Imaging, vol. 2008, Article ID 297089, 6 pages, 2008. https://doi.org/10.1155/2008/297089
Random Volumetric MRI Trajectories via Genetic Algorithms
Abstract
A pseudorandom, velocityinsensitive, volumetric space sampling trajectory is designed for use with balanced steadystate magnetic resonance imaging. Individual arcs are designed independently and do not fit together in the way that multishot spiral, radial or echoplanar trajectories do. Previously, it was shown that secondorder cone optimization problems can be defined for each arc independent of the others, that nulling of zeroth and higher moments can be encoded as constraints, and that individual arcs can be optimized in seconds. For use in steadystate imaging, sampling duty cycles are predicted to exceed 95 percent. Using such pseudorandom trajectories, aliasing caused by undersampling manifests itself as incoherent noise. In this paper, a genetic algorithm (GA) is formulated and numerically evaluated. A large set of arcs is designed using previous methods, and the GA choses particular fit subsets of a given size, corresponding to a desired acquisition time. Numerical simulations of 1 second acquisitions show good detail and acceptable noise for largevolume imaging with 32 coils.
1. Introduction
Reconstruction of magnetic resonance imaging (MRI) from data sampled using noncartesian sampling has recently received increasingly mathematically sophisticated treatment, for example [1], with notable improvements in reconstruction speed and accuracy. The application of novel design techniques for noncartesian sampling trajectories, however, has received less attention.
In [2], a novel pseudorandom volumetric space trajectory design method was presented. This methodology, henceforth referred to as Durga, combines a number of ideas in trajectory design and general sampling design for the first time, including randomness, [3, 4] constrained optimization [5] to balance trajectories for steadystate imaging [6–9], genetic algorithms [10], undersampling to trade acquisition time for (structured) noise [11–13], and targetoriented design rather than patterns of symmetric interleaving [14, 15]. By combining these ideas, Durga achieves significantly better efficiency as measured by sampling dutycycle for a balanced steadystate pulse sequence.
Flowinsensitive space trajectories are inherently more efficient than trajectories requiring a rewinder to balance first and possibly higher moments. Trajectories which sample three dimensions in space, like Durga, can further increase efficiency by not rewinding slice select gradients, and simply starting and stopping sampling offset from the center of space immediately after and before the excitation pulse. This is summarized in Table 1 by a comparison of spiral with velocity compensating rewinder [5], Teardrop [7] which incorporates inplane velocity compensation into the readout, and Durga [2]. These numbers are relative to gradient peak/slew limits of 40 mTm^{−1} and 150 Tm^{−1}s^{−1} for Spiral [5] and Durga [2] and 27 mTm^{−1} and 72 Tm^{−1}s^{−1} for Teardrop [7].

Volumetric imaging has several inherent advantages over thinslice imaging, including isotropic resolution, reduced effect of inflow, and the ability to completely correct for distortions in gradients, mainfield inhomogeneity, and eddy currents. It is also easier to completely compensate for velocity effects by nulling higher moments. Momentnulled planar trajectories are often only effectively nulled in one or two dimensions, because they are used with slice select and possibly phase encoding gradients which are difficult to null.
In the case of rapid imaging, however, the most important property of sampling in more dimensions is the ability to under sample without introducing aliasing in the form of ghosts [11, 16]. Using Durga, one can essentially trade off under sampling for unstructured noise on a continuous basis, breaking the dependence on Nyquist sampling. The trajectories used in the illustrations to reconstruct volumes are under sampled in time relative to perfect Cartesian sampling (ignoring gradient limits, and only considering sampling bandwidth).
This paper presents two innovations over [2]. (1)A genetic algorithm for choosing subsets of trajectory arcs (Figure 1) corresponding to a , designed using one of the methods in [2]. This significantly increases the simplicity and flexibility of implementation and improves the quality of the solution, as measured by the point spread function.(2)Multicoil numerical simulations, including the effect of noise, showing low levels of aliasing artifacts with extreme undersampling. Although this sampling strategy is compatible with and benefits from multicoil imaging, the reduction in aliasing possible by using an iterative SENSE reconstruction [17] is much smaller than what could be expected based on experience with other trajectories.
2. Methods
2.1. Individual Trajectory Design
Sampling trajectories in D space are generated by solving a secondorder cone optimization problem (SOCP), originally formulated in [2]. An SOCP is a minimization problem in which the variables are constrained to lie inside a set defined by quadratic functions of the variables. Putting a bound on the length of a vector is a common special case, and the one exploited in this trajectory design problem: where the variables, are the position in space at time ; controls the deviation from the constraints (1); is the time step; and are peak (3) and slew (4) constraints on the gradients; is the maximum resolution used in the reconstruction (5); the sum in (6) models the phase error caused by constant flow—errors from nonconstant flow can be handled by adding norms of higher moments; are the targets near which the trajectory should pass (7) at time , see Figure 2; and the s determine the relative priorities placed on meeting the different goals.
In the trajectories designed for this paper, the goals are random points on the boundary sphere or at . Multiple traversals of are not optimal for sampling, but included to facilitate calibration of gradient distortion and mainfield inhomogeneity. This problem nulls the first moment at the center of the rf pulse. Higher moments, nulling at different points in the trajectory, or bounding the size of the moment with a tolerance could be encoded in the same way. The trajectories in this paper are designed to begin and end at , for steadystate imaging, but this is not a limitation of the method. See [2] for a detailed description of the software used to solve this optimization problem, and a discussion of variations in the design, including iterative designs in which goals are placed at lowdensity points for previous trajectories.
Note that the constraint on the first moment will cause an optimal trajectory to traverse parts of space not near a shortest path between the goals, , which increases the coverage of space, but there are no explicit objectives for coverage, intertrajectory spacing, or trajectory length. In [2], coverage is handled by pseudorandom assignment of goals from the vertices of a regular triangulation of the limit sphere in space.
2.2. Combining Sets of Trajectories
In this paper, coverage is assured by optimizing its effect on the point spread function (psf), which captures sampling artifacts (exactly for uniform coils, and approximately for surface coils). The ratio of the norm of the central point in the psf and the norm of the rest of the psf is maximized over subsets of a large fixed set of trajectories, each optimized as described above. Since the set of subsets of a fixed size grows factorially with the number of trajectories, only a small subset can be tested. A genetic algorithm uses information about previous subsets to test only subsets likely to have better psfs.
This approach differs from formulations, where the variables (referred to as chromosomes in genetic algorithms) are local or nonlocal control points on individual trajectories [10]. Such approaches lead to a much larger search space, and increased complexity of each search step, since the solution of the SOCP problem is the most expensive part of the optimization. In the present method, one SOCP problem is solved per trajectory (and may be solved in parallel) before application of the GA search.
Individual trajectories are 5.6 milliseconds, including 0.1 millisecond for rf excitation, and designed with gradient limits typical of a whole body clinical imager (peak gradient of 40 mTm^{−1} and a gradient slew rate 150 Tm^{−1}s^{−1}). For ease of comparison, we will choose enough trajectories for 1 second of sampling (178 with the above parameters).
As a first step, many trajectories are designed by randomly selecting goals on the boundary sphere (which are interleaved with fixed goals at ). Of these, many are clearly less suitable than others. To filter out the least suitable, a small number are designed to determine an approximate threshold for the longest %, and the % closest to their own goals. Only trajectories above the thresholds are used in the initial set.
This initial large set is then partitioned into 400 sets of the desired size (168) to seed the GA. These seed sets are then evolved through 45 generations, using chromosomes consisting of integer indices into the initial set. An additional 10 trajectories were added using a “density threading” [2], whose aim is to increase sampling of voids along the plane by specifically creating trajectories to pass through areas of low sampling density. Computation time refers to a PowerMac G5 2.5, running a C program calling SOCP, http://www.stanford.edu/~boyd/socp, and the Genetic Algorithm Utility Library, http://gaul.sourceforge.net.
The objective of the optimization is to improve the quality of the psf, as measured by Fourier transforming the (weighted) sampling density, and dividing the central value in the psf by the total energy. The method of determining the density must be the same as used in reconstructing the images, since different methods can potentially produce different artifacts even with the same trajectories [18].
To calculate the fractional sampling, both in space and in time, nearest neighbour resampling was used to calculate a sampling density, and the number of nonzero elements in the space array was computed. For example, 1 second trajectory set, 1.54% of the voxels in a space had nonzero values, giving one measure of undersampling. This corresponds to 51.7% of the number of the samples that would be collected in 1 second with a receiver bandwidth of 500 kHz, giving a measure of the efficiency relative to the (unrealizable) maximum sampling rate.
2.3. Image Reconstruction
Numerical simulations were reconstructed by simulating irregularly sampled data (500 kHz sampling) along the trajectories. Solid phantoms were approximated by cubes. Receiver coils were approximated by magnetic dipoles, with sensitivities evaluated at the center of each phantom cube. In this approximation, the signal in space corresponding to each cube is a product of sinc functions, scaled by the complex sensitivity value.
The irregular samples were then divided by the sampling density in space and resampled to a regular grid using convolution.
For multiple coils with sensitivities , the transformed images , were combined with two methods: weightingbysensitivitiesand conjugate gradient reconstructions [17].
A lot of analysis has been done on resampling irregular data with density variations, see [18–20] and references therein. Convolution with a fixed kernel [21] is the simplest resampling method to implement and to analyze. The optimal piecewiselinear kernel described in [22] was used. Estimating the density to use for correction was done by resampling a constant set of data points using a widththree triangular kernel. The key property of this triangular function is that it is a partition of unity for samples spaced 1, 2, or 3 grid points apart.
3. Results
When the described algorithm was used to generate the 168shot volume representing 1 second of imaging time, 12 hours were taken to create a working base of 4000 trajectories using randomly generated goals. An additional 13000 trajectories were rejected by the threshold test. Another 72 hours were required for the GA to select a fit subset.
Figure 3 compares the results of the GA to a histogram of 3000 randomly formed subsets. The GA set is 7.8 standard deviations better than the set, and is shown as a banded bar. Also shown is the result of adding 10 densitythreaded trajectories [2].
The estimated sampling density used for density correction is shown in Figure 4. For display purposes, the corrected sampling density is shown, in which the uniformity of the bright pixels reflects the quality of the correction. Recall that only 1.54% of grid points in space contain a sample point at this reconstructed resolution.
(a)
(b)
A better measure of the quality of the set of trajectories, together with the density correction, is given by the psf (Figure 4), which predicts very little blurring and a low level of noiselike aliasing.
Because it is difficult to inspect more than a planar crosssection of the psf at a time, it is hard to appreciate the extent to which the noiselike aliasing accumulates in three dimensions. This is demonstrated using a numerical phantom consisting of four rings of varying sizes, Figure 5.
The phantom was chosen to be readily identifiable in crosssection, and in volume/surface rendering. The apparent noise will increase as a function of the total signaling volume; the phantom represents a midway point between contrastenhanced MR angiography (with a small blood volume producing signal) and anatomical imaging. Four different reconstructions of the central  crosssection are shown in Figure 6.
The uniform reconstruction shows significantly more alias noise in the periphery of the image than the 32 surfacecoil reconstruction. For measurement noise, this is as expected, since the coils are more sensitive near the periphery. For the noiselike aliasing artifact, a similar phenomenon is in operation. The true image is multiplied by the coil sensitivity, which is then corrected when performing multicoil combination. The artifacts are also modulated by the sensitivity, but since they are delocalized (by definition), they are not corrected during combination, and hence partially cancel each other out. This is more apparent where the sensitivities vary (the periphery) than where they are relatively uniform (the center).
The noisy multiple and single surface coil images (lower row) show typical loss of sensitivity by the center of the phantom for a single channel, which is reduced in the combined image. The 32coil noisy image shows visually equal noise levels to the first two images, but with different (higher) frequency composition.
Iterative methods have been successfully used in planar imaging to reduce spiral artifacts by using a priori information about multiple coil sensitivities, notably using the conjugate gradient method [17]. This does not produce visible reductions in aliasing in this case. What it can do is reducing the effect of blurring, as is visible in the reconstruction of a voxel uniform box. Figure 7 shows the result of the first step common to steepestdescent and the conjugate gradient method. Figure 8 shows the residual decreasing for the first step, and starting to diverge thereafter, with two coil configurations, and simple and complex phantoms. This is not very surprising. Adapted methods have been proposed for illposedness [24] and roundoff errors [25], but are beyond the scope of this paper.
A common crosssection is displayed of the original solution, the result after one step, and a scaled version of the gradient (using, in the electronic edition, different colors for positive and negative pixel values). What aliasing artifact is present is either unchanged or enhanced by the gradient step, while softening of the edges is clearly reduced.
4. Discussion
The psf and numerical phantoms presented here show that Durga is a very promising approach to designing trajectories for volumetric imaging. Pseudorandom trajectory design visibly eliminates coherent aliasing artifacts in numerical simulations.
Using the GA increases the flexibility of this method and increases the quality of the solutions. There is a marked difference in quality between randomly generating subsets and using the GA to improve a population. After filtering out inferior individual trajectories, randomly selecting subsets produced surprisingly little variation in the quality of the psf. Initial attempts in using genetic algorithms based on randomly chosen initial points and unfiltered trajectories were not able to beat the heuristic in [2], so success of the GA depends very much on a careful formulation, and even then, it is unlikely that the global optimum will be reached. Unfortunately, it is also considerably more expensive. Additional work should focus on improving individual trajectories before starting the GA.
The modest improvement produced by steepest descent means that the sampling patterns are so efficient that little additional information can be extracted by a parallel reconstruction taking advantage of geometric coil information, at least with the present coils. We conjecture that the randomness optimized to make the psf flat produces a large cluster of small eigenvalues in the operator being iterated in the CG step, causing the CG to begin diverging after one iteration.
We are planning modifications of the basic trajectory design to quantify the first effect, and working with partners to collect data to evaluate the second. In any case, multiple coil reconstruction using coil sensitivities does reduce the apparent noise, which is important to end users.
We will also investigate compressed sensing reconstructions [26] which require incoherent aliasing artifacts such as those presented in this work, because in such cases “randomness is too important to be left to chance” [26].
References
 T. Knopp, S. Kunis, and D. Potts, “A note on the iterative MRI reconstruction from nonuniform $k$space data,” International Journal of Biomedical Imaging, vol. 2007, Article ID 24727, 9 pages, 2007. View at: Publisher Site  Google Scholar
 C. K. Anand, A. T. Curtis, and R. Kumar, “Durga: a heuristicallyoptimized data collection strategy for volumetric magnetic resonance imaging,” Engineering Optimization, vol. 40, no. 2, pp. 117–136, 2008. View at: Publisher Site  Google Scholar
 K. S. Nayak and D. G. Nishimura, “Randomized trajectories for reduced aliasing artifact,” in Proceedings of the 6th Annual Meeting of the International Society for Magnetic Resonance in Medicine (ISMRM '98), p. 670, Sydney, Australia, April 1998. View at: Google Scholar
 R. Ahmad, Y. Deng, D. S. Vikram et al., “Quasi Monte Carlobased isotropic distribution of gradient directions for improved reconstruction quality of 3D EPR imaging,” Journal of Magnetic Resonance, vol. 184, no. 2, pp. 236–245, 2007. View at: Publisher Site  Google Scholar
 B. A. Hargreaves, D. G. Nishimura, and S. M. Conolly, “Timeoptimal multidimensional gradient waveform design for rapid imaging,” Magnetic Resonance in Medicine, vol. 51, no. 1, pp. 81–92, 2004. View at: Publisher Site  Google Scholar
 A. Oppelt, R. Graumann, H. Barfuss, H. Fischer, W. Hartl, and W. Shajor, “Fisp—a new fast mri sequence,” Electromedica, vol. 54, no. 1, pp. 15–18, 1986. View at: Google Scholar
 C. K. Anand, M. Thompson, D. Wu, and T. Cull, “Teardrop, a novel trajectory for truefisp,” in Proceedings of the 9th Annual Meeting of the International Society for Magnetic Resonance in Medicine (ISMRM '01), p. 1804, Glasgow, UK, April 2001. View at: Google Scholar
 C. K. Anand, T. Ren, and T. Terlaky, “Optimizing Teardrop, an MRI sampling trajectory,” Optimization Methods and Software. In press. View at: Publisher Site  Google Scholar
 K. S. Nayak, B. A. Hargreaves, B. S. Hu, D. G. Nishimura, J. M. Pauly, and C. H. Meyer, “Spiral balanced steadystate free precession cardiac imaging,” Magnetic Resonance in Medicine, vol. 53, no. 6, pp. 1468–1473, 2005. View at: Publisher Site  Google Scholar
 B. M. Dale, J. S. Lewin, and J. L. Duerk, “Optimal design of $k$space trajectories using a multiobjective genetic algorithm,” Magnetic Resonance in Medicine, vol. 52, no. 4, pp. 831–841, 2004. View at: Publisher Site  Google Scholar
 K. K. Vigen, D. C. Peters, T. M. Grist, W. F. Block, and C. A. Mistretta, “Undersampled projectionreconstruction imaging for timeresolved contrastenhanced imaging,” Magnetic Resonance in Medicine, vol. 43, no. 2, pp. 170–176, 2000. View at: Publisher Site  Google Scholar
 J. Du, S. B. Fain, F. R. Korosec, T. M. Grist, and C. A. Mistretta, “Timeresolved contrastenhanced carotid imaging using undersampled projection reconstruction acquisition,” Journal of Magnetic Resonance Imaging, vol. 25, no. 5, pp. 1093–1099, 2007. View at: Publisher Site  Google Scholar
 J. Du, F. J. Thornton, S. B. Fain et al., “Artifact reduction in undersampled projection reconstruction MRI of the peripheral vessels using selective excitation,” Magnetic Resonance in Medicine, vol. 51, no. 5, pp. 1071–1076, 2004. View at: Publisher Site  Google Scholar
 R. Mir, A. Guesalaga, J. Spiniak, M. Guarini, and P. Irarrazaval, “Fast threedimensional $k$space trajector design using missile guidance ideas,” Magnetic Resonance in Medicine, vol. 52, no. 2, pp. 329–336, 2004. View at: Publisher Site  Google Scholar
 J. Spiniak, A. Guesalaga, R. Mir, M. Guarini, and P. Irarrazaval, “Undersampling $k$space using fast progressive 3D trajectories,” Magnetic Resonance in Medicine, vol. 54, no. 4, pp. 886–892, 2005. View at: Publisher Site  Google Scholar
 A. V. Barger, W. F. Block, Y. Toropov, T. M. Grist, and C. A. Mistretta, “Timeresolved contrastenhanced imaging with isotropic resolution and broad coverage using an undersampled 3D projection trajectory,” Magnetic Resonance in Medicine, vol. 48, no. 2, pp. 297–305, 2002. View at: Publisher Site  Google Scholar
 K. P. Pruessmann, M. Weiger, P. Börnert, and P. Boesiger, “Advances in sensitivity encoding with arbitrary $k$space trajectories,” Magnetic Resonance in Medicine, vol. 46, no. 4, pp. 638–651, 2001. View at: Publisher Site  Google Scholar
 J. G. Pipe and P. Menon, “Sampling density compensation in MRI: rationale and an iterative numerical solution,” Magnetic Resonance in Medicine, vol. 41, no. 1, pp. 179–186, 1999. View at: Publisher Site  Google Scholar
 F. Wajer, R. Lethmate, J. van Osch, D. GraveronDemilly, and D. van Ormondt, “Simple formula for the accuracy of gridding,” in Proceedings of the 9th Annual Meeting of the International Society for Magnetic Resonance in Medicine (ISMRM '01), Glasgow, UK, April 2001. View at: Google Scholar
 V. Rasche, R. Proksa, R. Sinkus, P. Börnert, and H. Eggers, “Resampling of data between arbitrary grids using convolution interpolation,” IEEE Transactions on Medical Imaging, vol. 18, no. 5, pp. 385–392, 1999. View at: Publisher Site  Google Scholar
 J. I. Jackson, C. H. Meyer, D. G. Nishimura, and A. Macovski, “Selection of a convolution function for Fourier inversion using gridding,” IEEE Transactions on Medical Imaging, vol. 10, no. 3, pp. 473–478, 1991. View at: Publisher Site  Google Scholar
 C. K. Anand, T. Terlaky, and B. Wang, “Rapid, embeddable design method for spiral magnetic resonance image reconstruction resampling kernels,” Optimization and Engineering, vol. 5, no. 4, pp. 485–502, 2004. View at: Publisher Site  Google Scholar
 A. Rosset, L. Spadola, and O. Ratib, “OsiriX: an opensource software for navigating in multidimensional DICOM images,” Journal of Digital Imaging, vol. 17, no. 3, pp. 205–216, 2004. View at: Publisher Site  Google Scholar
 B. Eicke, A. K. Louis, and R. Plato, “The instability of some gradient methods for illposed problems,” Numerische Mathematik, vol. 58, no. 1, pp. 129–134, 1990. View at: Publisher Site  Google Scholar
 Z. Strakoš and P. Tichý, “On error estimation in the conjugate gradient method and why it works in finite precision computations,” Electronic Transactions on Numerical Analysis, vol. 13, pp. 56–80, 2002. View at: Google Scholar
 M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: the application of compressed sensing for rapid MR imaging,” Magnetic Resonance in Medicine, vol. 58, no. 6, pp. 1182–1195, 2007. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2008 Andrew Thomas Curtis and Christopher Kumar Anand. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.