|Home | About | Journals | Submit | Contact Us | Français|
A fast and simple algorithm for designing time-optimal waveforms is presented. The algorithm accepts a given arbitrary multi-dimensional k-space trajectory as the input and outputs the time-optimal gradient waveform that traverses k-space along that path in minimum time. The algorithm is non-iterative, and its run time is independent of the complexity of the curve, i.e. the number of switches between slew-rate limited acceleration, slew-rate limited deceleration and gradient amplitude limited regions. The key in the method is that the gradient amplitude is designed as a function of arc length along the k-space trajectory, rather than as a function of time. Several trajectory design examples are presented.
Advances in the field of Magnetic Resonance Imaging (MRI) such as gradient hardware, high field systems, optimized receiver coil arrays, fast sequences and sophisticated reconstruction methods provide the ability to image faster than ever. New acquisition methods are being explored in which k-space is scanned in nontraditional trajectories –.
One of the design challenges in rapid imaging is to minimize the gradient waveform duration, subject to both hardware and sequence constraints. The preferred solution would be to design both the k-space trajectory and the gradient waveforms at the same time. However, because of the complexity and the large number of control variables, often only approximate heuristic methods ,  are used as the optimal solution is generally not tractable.
A simpler, more common approach is to first choose a sampling trajectory and then design the gradient waveforms for it. For example, spirals are often designed this way , –. In this approach, the problem is to find the gradient waveform that will traverse k-space from one point to another along a specific path and in minimum-time, while meeting the hardware gradient magnitude and slew-rate constraints. For trajectories such as linear, circular  or spirals , – the solution is quite simple; first operate in the slew-rate limited regime till the maximum gradient is reached, and then operate in the maximum gradient regime. For the general arbitrary trajectory, the solution is not simple anymore because there can be numerous switching points between slew-limited acceleration, slew-limited deceleration and gradient-limited regimes. In the current literature, there is no general methodology to design gradient waveforms for arbitrary trajectories that guarantees time-optimality and is computationally inexpensive.
It is important to mention that some of the optimal designs that exist in the literature either solve for 1D waveforms , , or provide waveforms that traverse k-space from one point to another but not on a specific path , .
In this paper we develop a fast and simple algorithm based on optimal control theory that provides the complete solution to the time-optimal gradient waveform for arbitrary k-space paths. Using this method a user need only prescribe a parametric curve in k-space (arbitrary parametrization) and the algorithm will output the gradient waveform that traverses the k-space path in minimum time.
We start by describing some properties of planar and volumetric curves that are essential to the derivation of the time-optimal gradient waveform design.
Suppose we are given a specified path, curve, or trajectory C from C0 to C1 in k-space. Suppose the curve C is described as a function of some parameter p:
Here, p = 0 corresponds to the initial point and p = pmax corresponds to the end point:
The first derivative of the curve with respect to its parametrization is the velocity or tangent vector of the parametrization, which is denoted as
The second derivative of the curve with respect to its parametrization is the acceleration vector of the parametrization
From here onwards we denote h′(p) as the derivative of the function h with respect to a general parameter p. We make an exception when using the notation to specifically indicate that it is a time derivative.
A very useful parametrization is the Euclidean arc-length s parametrization:
where L is the length of the path. This parametrization describes the coordinates as a function of the Euclidean distance along the curve. An important property of this parametrization is that it has unit velocity, i.e.,
for all s. Another important property is that the magnitude of the acceleration is the curvature κ(s) of the curve, i.e.
The curvature of the curve at a given point is the reciprocal of the radius of an osculating circle that has the same first and second derivatives with the curve at that point. Figure 1 illustrates the properties of the arc-length parametrization.
When a curve is given in an arbitrary parametrization C(p) = (x(p), y(p), z(p)), it is always possible to convert into the arc-length parametrization by using the relation
In this paper we aim to design a gradient waveform as a function of time. This is equivalent to designing a time parametrization of the curve that describes the k-space coordinates as a function of time. Specifically, we aim to design a time function p = s(t) in the arc-length parametrization such that
where T is the traversal time. The time trajectory in k-space is given by the composite function (t) = C(s(t)).
First, we derive the relation between the gradient waveform and the curve parametrization. The gradient waveform is proportional to the velocity in the time parametrization and is given by
where γ is the gyro-magnetic ratio. Here, we use (t) to indicate the time derivative of s. The gradient slew-rate is proportional to the acceleration in the time parametrization. Using the chain rule we obtain
The design variable in the MRI system is the gradient waveform. The gradients are subject to hardware as well as sequence constraints. For the hardware constraints we assume the frequently used slew-limited model as described in . In this model, the gradient amplitude is subject to the maximum amplitude of the system
It is also subject to the maximum slew-rate of the system
For the sequence constraints, the gradient waveform is constrained to follow a specific trajectory in k-space such that
It is also constrained to have an initial value. For simplicity, we assume an initial value of
Other optional sequence constraints such as final or intermediate values are possible, but are not assumed here. We discuss some of these optional constraints later in this section.
Now that we have derived the hardware as well as the sequence constraints, we consider the problem of finding the time-optimal gradient waveform that satisfies them. The time-optimal problem can be formulated as
Here the variable is the gradient g(t) defined over the time interval [0, T]. The objective is to minimize the traversal time T along the trajectory.
The time-optimal solution is always either gradient or slew-rate limited. Solving for the time-optimal waveform in this formulation requires one to find the optimal switching times between the maximum gradient and maximum slew-rate. This procedure is difficult for complex curves. Instead, we look at an alternative equivalent formulation in which the solution becomes simpler.
Let s denote the arc length from the initial point. Because the k-space path is given as a constraint, we need only to design the time function s(t). Note that s is always increasing, so
This means that it is sufficient to design the gradient magnitude along the path.
A complete derivation of (14) and (15) is deferred to the Appendix. Intuitively, the constraint in (14) accounts for the geometry of the trajectory, and is related to the maximum velocity at which a curve can be approached without violating the acceleration constraint, and is independent of past or future velocities. Equation (15) is a dynamic constraint: It is a differential inequality that describes the allowed change in the velocity at a specific point on the path given the velocities in its proximity.
Up until now, we have only provided the formulation of the problem. Here, we provide a complete solution to (16). The solution is obtained in the velocity vs. arc-length plane ( vs. s). This plane is often referred to as the phase plane in optimal-control theory literature –. The outline is the following: We first find the optimal velocity as a function of arc length (υ⭑(s)) and then find the optimal time function s⭑(t) which can be used to derive the optimal gradient waveform. (This reformulation method has been used in solving time-optimal path planning problems for robotic manipulators along specified paths –.)
In the phase-plane, we represent the velocity as a function of arc length:
Note that the traversal time T is a function of the velocity given by
Also, note that the acceleration is also a function of the velocity given by
Then, the time-optimal control problem (16) amounts to solving the following optimization problem in the phase-plane:
where the optimization variable is the function υ(s) defined over [0,L].
The optimal solution υ⭑ to this problem describes the relation between the optimal time function s⭑ and its derivative in the phase plane:
Using this relation, we can readily recover s⭑ from υ⭑.
We now describe a complete solution to the optimization problem (18). In order to find the optimal velocity, we need to integrate two ordinary differential equations (ODEs). The first ODE is given by
which we integrate forward with the initial condition υ+(0) = 0. The second ODE is
which we integrate backwards with the final condition υ−(L) = υ+(L). (Recall that α(s) and β(s, ) are given in (15).)
The optimal velocity υ⭑ (s) is simply given by
A proof of optimality is given in the Appendix.
Here we describe the procedure to recover the time-optimal gradient waveform g⭑(t) from the optimal velocity υ⭑ (s). We can obtain s⭑(t) by computing the inverse function of t⭑ (s), using the relation
In particular, the traversal time, which is the optimal value of (16), is given by
It follows from (8) that the time-optimal gradient waveform is
A summary of the design algorithm is given in Fig. 2. Figs. 3 through through55 show a simplified example of the stages of the algorithm applied to a simple trajectory. Figure 3 shows the transition from the arbitrary parametrization to the arc-length parametrization. Figure 4 shows the calculation of the gradient “velocity” as a function of arc length in the phase-plane. Figure 5 shows the transition from the arc-length parametrization to the time parametrization and the calculation of the gradient waveforms. Although the trajectory in Figs. 3–5 has little practical value, it was chosen for its educational value as it shows clearly and simply the stages of the design.
To handle nonzero initial gradient g(0), it suffices to integrate the ODE (19) forward with the initial condition υ+(0) = γ|g(0)| instead of υ+(0) = 0. However, if the value of g(0) is infeasible, the outcome of the design will be the maximum feasible one.
In the same way, to handle final gradient value g(L), it suffices to integrate the ODE (20) backward with the final condition υ−(0) = γ|g(L)| instead of υ−(L) = υ+(L). Again, if the value is infeasible, the outcome will be the maximum feasible one.
In this section we present a few examples demonstrating some of the applications of our method to gradient waveform design.
The design algorithm was implemented in Matlab (The MathWorks, Inc., Natick, MA, USA) and in the C programming language. All simulations were performed on a Mandriva Linux workstation with an AMD Athlon 3800+ 64bit processor and 2.5GB memory. Derivative operations were approximated by finite differences. Numerical integrations were approximated by the trapezoid method. The ODEs were solved using a 4th order Runge-Kutte method . We used the cubic-spline interpolation method for interpolating the curve when needed. In all our designs, we assumed gradients capable of 40 mT=m, slew-rate of 150 mT=m=ms and a sampling rate of 4 µs.
A Matlab and a C programming language implementation of the algorithm as well as the function to reproduce the results will be available upon publication at http://www-mrsrl.stanford.edu/~mlustig/software/.
As a sanity test we have applied the algorithm to a horizontal line in k-space. Constraining zero initial and final gradient amplitude we expected the resulting gradient waveform to be a slew and gradient limited trapezoid. The k-space line was chosen to be 10 cm−1 long. The result of the experiment are shown in Fig. 6a. As expected, the resulting gradient waveform is a trapezoid. Note, that the trapezoid is either slew limited or gradient limited, which guarantees time optimality in this simple case.
Spiral trajectories are used in real-time and rapid imaging applications. In some applications, for example MR fluoroscopy  and autocalibrated parallel imaging , spirals with varying densities are desired. Here we demonstrate an alternative design by a simple curve parametrization followed by the time-optimal gradient design.
The spiral is parametrized using the following equation
where Nitlv is the number of spiral interleaves and FOV (r) is a function that describes the supported field of view (FOV) as a function of radius in k-space.
To test our method, we designed a dual density spiral with 16 interleaves, 0.83mm resolution (maximum k-space radius of 6cm−1), a FOV of 55cm for r [0, 1.2]cm−1 and 10cm for r [1.8, 6]cm−1. The transition region of r [1.2, 1.8]cm−1 was linearly interpolated. The results of the design for a single interleaf are presented in Fig. 6b. As expected the gradient waveforms are gradient and slew limited. The transition between the two density regions is clearly seen in the magnitude gradient plot.
The rosette trajectory was first introduced in  as a spectrally selective imaging method. Because of the multiple crossings of the trajectory, off-resonance often causes destructive interference and therefore off-resonance voxels do not show up in the image. Recently, rosettes have been reconsidered as a rapid sampling trajectory for CS-MRI with application to hyperpolarized 13C imaging , where scan time is crucial.
In  Noll describes rosettes as the parametric curve in time
The gradient waveforms according to Noll are
The advantage of this design is that there is an analytic expression for the gradient waveform, and no further calculations are needed. However, it has a couple of limitations. The first is that the waveforms are neither gradient limited nor slew limited and are not time-optimal. The scan efficiency is especially reduced when designing high resolution trajectories, where the waveforms constraints are set mostly by the maximum gradient amplitude. The second is that at t = 0, Eq. (27) does not have a zero value, therefore a correction is needed to ramp the gradients from zero initial value to meet the slew-rate constraints. Instead, we use the parametrization of Eq. (26) and apply our design method to get the desired time-optimal gradient waveforms. To test the design we used ω1 = 1.419, ω2 = 0.8233, kmax = 12 as the parameters in Eq. (26). The results of applying our design are presented in Fig. 6c. The faint lines in the figure show the result for Noll’s design. The figure shows that as expected, the waveforms of our proposed method are gradient and slew limited. It also shows that for this specific trajectory Noll’s design requires a traversal time of 11.32 ms. This is a 20% increase over our time-optimal design which requires only 9.46 ms.
It was shown in  that under-sampled randomly perturbed spiral trajectories can provide faster imaging when used with a special non-linear reconstruction. When under-sampling to save scan time, it is essential that the gradient waveforms be time-optimal. Therefore, we applied our design algorithm to a randomly perturbed variable density spiral where the optimal waveform requires numerous switching between acceleration-deceleration slew-rate limited regions and gradient magnitude limited regions. We designed a 4 interleave variable density spiral, chosen to have a resolution of 1mm and a FOV of 20cm at the k-space origin, that linearly decreases to 5cm on the periphery. The spiral was perturbed in the radial direction by a randomly generated smooth waveform Figure 6d shows the result of the randomly perturbed spiral design. Again, as expected the waveforms are gradient and slew limited.
The results in the previous section show the main advantage of our design method; it can be applied to any type of k-space trajectory: from simple rectilinear to the complicated randomized trajectories. The trajectory and gradient design is simplified to designing a parametric curve, and then designing the gradients for it. The latter, is also a disadvantage of our method, since a poor choice of trajectory that has sharp curvatures may lead to a poor design. Nevertheless, as shown in the examples, the simplicity and generality of the method makes it a powerful design tool.
The differential equations Eq. (19) and Eq. (20) have the gradient magnitude in the denominator. When the gradient magnitude is very small, the right hand side of the equation can become very large and cause inaccuracies in the integration. This can be mitigated by choosing a smaller step size, which increases the computation time and memory consumption. However, since slowing gradients are associated with large curvature locations in the k-space trajectory that are known in advance, one can use a variable step size that is adapted to the curvature.
The proposed method is non-iterative and provides a direct time-optimal solution. The computational complexity is linear with respect to the length of the trajectory curve, and requires a solution to an ODE propagated forward and backwards in time. Our implementation used a fixed step 4th order Runge-Kutte method, which required a relatively small step size to maintain accuracy. Nevertheless, our C programming language implementation required about half a second to design the waveforms in the examples.
We have provided a fast, simple and non-iterative method for designing the time-optimal gradient waveforms for any k-space trajectory. It is the complete solution for the gradient waveform design with a k-space path constraint. We have demonstrated some simple time-optimal waveforms by first designing a parametric curve and using the design method to produce the gradient waveforms.
This work was supported by NIH grants R01 HL074332, R01 EB002992, R01 HL067161 and GE Medical Systems.
We start by observing that
The tangent and the normal vectors are always orthogonal, that is, C″ (s(t))TC′ (s(t)) = 0. Also, recall that |C′ (s(t))| = 1 and that |C″ (s(t))| = κ(s(t)). The slew-rate constraint can therefore be expressed as
which in turn implies that
The right-hand side is non-negative when
Finally, the gradient magnitude is limited |g(t)| ≤ Gmax resulting in the desired second constraint,
Let υ denote any velocity profile along the trajectory that satisfies the constraints of (18). Recall that the traversal time Tυ of υ is given by
It follows from Lemma 3 in  that
Michael Lustig, Magnetic Resonance Systems Research Laboratory, Stanford University, Stanford, CA 94305 ; Email: ude.drofnats.lrsrm@gitsulm..
Seung-Jean Kim, Information Systems Laboratory, Stanford University, Stanford, CA 94305 (Email: ude.drofnats@mikjs).
John M. Pauly, Magnetic Resonance Systems Research Laboratory, Stanford University, Stanford, CA 94305 ; Email: ude.drofnats@yluap.