# Random Volumetric MRI Trajectories via Genetic Algorithms

Authors: Andrew Thomas Curtis; Christopher Kumar Anand
Journal: International Journal of Biomedical Imaging
 (2008)

---

## Abstract

A pseudorandom, velocity-insensitive, volumetrick-space sampling trajectory is designed for use with balanced steady-state magnetic resonance imaging. Individual arcs are designed independently and do not fit together in the way that multishot spiral, radial or echo-planar trajectories do. Previously, it was shown that second-order 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 steady-state imaging, sampling duty cycles are predicted to exceed 95 percent. Using such pseudorandom trajectories, aliasing caused by under-sampling 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 large-volume imaging with 32 coils.


---

## Body

## 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 k-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 steady-state imaging [6–9], genetic algorithms [10], under-sampling to trade
acquisition time for (structured) noise [11–13], and target-oriented design rather than patterns of
symmetric interleaving [14, 15]. By combining these ideas, Durga achieves
significantly better efficiency as measured by sampling duty-cycle for a
balanced steady-state pulse sequence.Flow-insensitivek-space trajectories are inherently more
efficient than trajectories requiring a rewinder to balance first and possibly
higher moments. Trajectories which sample three dimensions in k-space, like Durga, can further increase
efficiency by not rewinding slice select gradients, and simply starting and
stopping sampling offset from the center of k-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 in-plane velocity compensation into
the readout, and Durga [2]. These numbers are relative to gradient peak/slew
limits of 40 mTm−1 and 150 Tm−1s−1 for Spiral [5] and Durga [2]
and 27 mTm−1 and 72 Tm−1s−1 for Teardrop [7].Table 1
Duty-cycle
comparison for three balancedk-space sampling patterns.
SpiralTeardropDurgaReadout (ms)2.403.435.50Excitation (ms)1.201.470.10Rewinder (ms)1.40——TR (ms)5.004.935.60Duty cycle0.480.780.98Volumetric imaging has several inherent advantages
over thin-slice imaging, including isotropic resolution, reduced effect of
in-flow, and the ability to completely correct for distortions in gradients,
main-field inhomogeneity, and eddy currents. It is also easier to completely
compensate for velocity effects by nulling higher moments. Moment-nulled 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 2563 volumes are 33× 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 (Figure1) corresponding to a TR, 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
under-sampling.
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.Figure 1
Tubes depicting the178 selected trajectories.


## 2. Methods

### 2.1. Individual Trajectory Design

Sampling trajectories in3D k-space are generated by solving a second-order
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:
(1)minτ,(2)k0=0=kN,(3)∥ki+1−ki∥2≤GmaxδT,(4)∥ki+2−2ki+1+ki∥2≤SmaxδT2,(5)∥ki∥2≤1/R,(6)∥∑i=1Ni(ki+1−ki)∥<λflowτ,(7)∥kti−Gi∥2<{λ0τGi=0,λgoalτotherwise, where the variables, ki∈ℝ3,i=1,…,N−1, are the position in k-space at time i; τ controls the deviation from the constraints
(1); δT is the time step; Gmax and Smax are peak (3) and slew (4) constraints on the
gradients; R 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; Gi are the targets near which the trajectory
should pass (7) at time tj,
see Figure 2; and the λs determine the relative priorities placed on
meeting the different goals.Figure 2
Diagram showing a trajectory going through specified
goals at maximum resolution.δ is a gap in the sampling, chosen to correspond
with the rf excitation.In the trajectories designed for this paper, the goalsGj are random points on the boundary sphere |k|=1/R or at k=0.
Multiple traversals of k=0 are not optimal for sampling, but included to
facilitate calibration of gradient distortion and main-field 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 k=0, for steady-state 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 low-density points for previous trajectories.Note that the
constraint on the first moment will cause an optimal trajectory to traverse
parts ofk-space not near a shortest path between the
goals, Gi,
which increases the coverage of k-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 k-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 theL2 norm of the central point in the psf and the L2 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−1s−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 atk=0). 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 20%, and the 20% 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 kx−ky 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.Theobjective 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 ink-space and in time, nearest neighbour
resampling was used to calculate a sampling density, and the number of nonzero
elements in the k-space array was computed. For example, 1 second trajectory set, 1.54% of the voxels in a 2563k-space had nonzero values, giving one measure
of under-sampling. 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 ink-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 ink-space and resampled to a regular grid using
convolution.For multiple coils with sensitivitiesSi, the transformed images ρi, were combined with two methods: weighting-by-sensitivities(8)∑coilsSi¯⋅ρi∑coilsSi¯⋅Siand 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
piecewise-linear 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
width-three 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.


## 2.1. Individual Trajectory Design

Sampling trajectories in3D k-space are generated by solving a second-order
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:
(1)minτ,(2)k0=0=kN,(3)∥ki+1−ki∥2≤GmaxδT,(4)∥ki+2−2ki+1+ki∥2≤SmaxδT2,(5)∥ki∥2≤1/R,(6)∥∑i=1Ni(ki+1−ki)∥<λflowτ,(7)∥kti−Gi∥2<{λ0τGi=0,λgoalτotherwise, where the variables, ki∈ℝ3,i=1,…,N−1, are the position in k-space at time i; τ controls the deviation from the constraints
(1); δT is the time step; Gmax and Smax are peak (3) and slew (4) constraints on the
gradients; R 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; Gi are the targets near which the trajectory
should pass (7) at time tj,
see Figure 2; and the λs determine the relative priorities placed on
meeting the different goals.Figure 2
Diagram showing a trajectory going through specified
goals at maximum resolution.δ is a gap in the sampling, chosen to correspond
with the rf excitation.In the trajectories designed for this paper, the goalsGj are random points on the boundary sphere |k|=1/R or at k=0.
Multiple traversals of k=0 are not optimal for sampling, but included to
facilitate calibration of gradient distortion and main-field 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 k=0, for steady-state 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 low-density points for previous trajectories.Note that the
constraint on the first moment will cause an optimal trajectory to traverse
parts ofk-space not near a shortest path between the
goals, Gi,
which increases the coverage of k-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 k-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 theL2 norm of the central point in the psf and the L2 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−1s−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 atk=0). 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 20%, and the 20% 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 kx−ky 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.Theobjective 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 ink-space and in time, nearest neighbour
resampling was used to calculate a sampling density, and the number of nonzero
elements in the k-space array was computed. For example, 1 second trajectory set, 1.54% of the voxels in a 2563k-space had nonzero values, giving one measure
of under-sampling. 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 ink-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 ink-space and resampled to a regular grid using
convolution.For multiple coils with sensitivitiesSi, the transformed images ρi, were combined with two methods: weighting-by-sensitivities(8)∑coilsSi¯⋅ρi∑coilsSi¯⋅Siand 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
piecewise-linear 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
width-three 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 168-shot 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.Figure3 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 density-threaded trajectories [2].Figure 3
Comparison of randomly selected trajectory sets versus
solutions from the GA. Note that the banded bars indicate horizontal position
not counts.The estimated sampling density used for density
correction is shown in Figure4. 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 k-space contain a sample point at this
reconstructed resolution.Above: cross-section of the computed sampling
density used in reconstruction. (Below) cross-section of the psf,
expanded using 8X Fourier interpolation.
(a)(b)A better measure of the quality of the set of
trajectories, together with the density correction, is given by the psf (Figure4), which
predicts very little blurring and a low level of noise-like aliasing.Because it is difficult to inspect more than a planar
cross-section of the psf at a time, it is hard to appreciate the extent to
which the noise-like aliasing accumulates in three dimensions. This is
demonstrated using a numerical phantom consisting of four rings of varying
sizes, Figure5.Figure 5
Surface rendering (OsiriX [23]) of solid rings
reconstructed with 178 arcs and 32 noisy coils.The phantom was chosen to be readily identifiable in cross-section, 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
contrast-enhanced MR angiography (with a small blood volume producing signal)
and anatomical imaging. Four different reconstructions of the centralx-y cross-section are shown in Figure 6.Figure 6
Four identical cross-sections showing, reconstruction
with a uniform coil (top left), with 32 dipole coils (top right), with one dipole coil
including noise (bottom right), and 32 dipole coils including noise (bottom left).The uniform reconstruction shows
significantly more alias noise in the periphery of the image than the 32 surface-coil reconstruction. For measurement
noise, this is as expected, since the coils are more sensitive near the
periphery. For the noise-like 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 32-coil 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 163 voxel uniform box. Figure 7 shows the result of
the first step common to steepest-descent 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 ill-posedness [24] and round-off errors
[25], but are beyond
the scope of this paper.Figure 7
Reconstruction of a cubic phantom using 178 trajectories. (Left to right)
Original 32-channel reconstruction, image showing some edge enhancement after
one gradient step, and the gradient (not to scale).Figure 8
L2 difference between true image and reconstruction
after conjugate gradient iterations.A common cross-section 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].



---

*Source: 297089-2008-07-01.xml*