Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
IEEE Trans Med Imaging. Author manuscript; available in PMC 2010 August 26.
Published in final edited form as:
PMCID: PMC2928662

A Fast Method for Designing Time-Optimal Gradient Waveforms for Arbitrary k-Space Trajectories

Michael Lustig, Student Member, IEEE, Seung-Jean Kim, Member, IEEE, and John M. Pauly, Member, IEEE


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 [1]–[6].

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 [7], [8] 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 [1], [9]–[11]. 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 [12] or spirals [1], [9]–[11] 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 [13], [14], or provide waveforms that traverse k-space from one point to another but not on a specific path [15], [16].

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.


A. Curve Parametrization

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 h 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.

Fig. 1
Properties of the arc-length parametrization of a curve in R3.

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


B. The Time-Optimal Gradient Waveform Design Problem

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 C(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 s(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 [16]. 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.

C. Formulation in the Arc-Length Parametrization

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


It follows from (5) and (8) that


This means that it is sufficient to design the gradient magnitude along the path.

We start with the formulation of the hardware constraints in the arc-length parametrization. It follows from (8), (9), (10), and (11) that





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.

Finally, problem (12) can be equivalently formulated in the arc-length parametrization as


Here the variable is the time function s(t). Once we know the optimal solution s[open star](t) of problem (16), we can find the solution to the original problem (12) using


D. The Time-Optimal Solution in the Phase-Plane

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 (s vs. s). This plane is often referred to as the phase plane in optimal-control theory literature [17]–[19]. 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 [17]–[19].)

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, s) are given in (15).)

The optimal velocity υ (s) is simply given by


A proof of optimality is given in the Appendix.

E. The Time-Optimal Gradient Waveform Solution

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. 35 has little practical value, it was chosen for its educational value as it shows clearly and simply the stages of the design.

Fig. 2
Outline of the time-optimal gradient design algorithm.
Fig. 3
Example of stage I of the design algorithm: Conversion from arbitrary parametrization to arc-length parametrization.
Fig. 4
Example of stage II of the design algorithm: Calculation of the time-optimal “velocity” in the phase-plane.
Fig. 5
Example of stage III of the design algorithm: Conversion from the arc-length parametrization to the time-optimal parametrization and the calculation of the gradient waveforms. Note that either the gradient magnitude or the slew-rate are maximized at every ...

F. Additional Constraints

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 general, intermediate gradient and slew-rate magnitude constraint can be applied by replacing the constraints in Eq. (10) and Eq. (11) to be a function of the arc length, s:





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 [20]. 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

A. A Line in k-space

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.

Fig. 6
Examples showing several k-space trajectories with their corresponding time-optimal gradient waveforms and gradient magnitude and slew-rate. (a) A line. (b) Dual density spiral. (c) Rosette trajectory. (d) Randomly perturbed variable density spiral.

B. Dual Density Spiral Design

Spiral trajectories are used in real-time and rapid imaging applications. In some applications, for example MR fluoroscopy [21] and autocalibrated parallel imaging [22], 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 [set membership] [0, 1.2]cm−1 and 10cm for r [set membership] [1.8, 6]cm−1. The transition region of r [set membership] [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.

C. Time-Optimal Rosette Trajectory Design

The rosette trajectory was first introduced in [5] 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 [23], where scan time is crucial.

In [5] 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.

D. Randomized Trajectories for Compressed Sensing Application

It was shown in [4] 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.

A. Numerical issues

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.

B. Computational Complexity

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 [18] that


Therefore, TυTυ. It can be shown through some arguments similar to those used to prove Theorem 3 in [18] that υ satisfies the constraints of (18). Therefore, υ is optimal for (18).

Contributor Information

Michael Lustig, Magnetic Resonance Systems Research Laboratory, Stanford University, Stanford, CA 94305 ; ude.drofnats.lrsrm@gitsulm..

Seung-Jean Kim, Information Systems Laboratory, Stanford University, Stanford, CA 94305 (ude.drofnats@mikjs).

John M. Pauly, Magnetic Resonance Systems Research Laboratory, Stanford University, Stanford, CA 94305 ; ude.drofnats@yluap.


1. Meyer CH, Hu BS, Nishimura DG, Macovski A. Fast spiral coronary artery imaging. Magn Reson Med. 1992;vol. 28(no. 2):202–213. [PubMed]
2. Pipe JG. An optimized center-out k-space trajectory for multishot MRI: Comparison with spiral and projection reconstruction. Magn Reson Med. 1999;vol. 42(no. 4):714–720. [PubMed]
3. Scheffler K, Hennig J. Frequency resolved single-shot MR imaging using stochastic k-space trajectories. Magn Reson Med. 1996 Apr;vol. 35(no. 4):569–576. [PubMed]
4. Lustig M, Lee JH, Donoho DL, Pauly JM. Faster imaging with randomly perturbed, under-sampled spirals and l1 reconstruction. Proceedings of the 13th Annual Meeting of ISMRM; Miami Beach. 2005. p. 685.
5. Noll DC. Multishot rosette trajectories for spectrally selective MR imaging. IEEE Trans Med Imaging. 1997 Aug;vol. 16(no. 4):372–377. [PubMed]
6. Gurney PT, Hargreaves BA, Nishimura DG. Design and analysis of a practical 3D cones trajectory. Magn Reson Med. 2006 Mar;vol. 55(no. 3):575–582. [PubMed]
7. Dale BM, Lewin JS, Duerk JL. Optimal design of k-space trajectories using a multi-objective genetic algorithm. Magn Reson Med. 2004 Oct;vol. 52(no. 4):831–841. [PubMed]
8. Mir R, Guesalaga A, Spiniak J, Guarini M, Irarrazaval P. Fast three-dimensional k-space trajectory design using missile guidance ideas. Magn Reson Med. 2004 Aug;vol. 52(no. 2):329–336. [PubMed]
9. King KF, Foo TK, Crawford CR. Optimized gradient waveforms for spiral scanning. Magn Reson Med. 1995;vol. 34(no. 2):156–160. [PubMed]
10. Glover GH. Simple analytic spiral k-space algorithm. Magn Reson Med. 1999;vol. 42(no. 2):412–415. [PubMed]
11. Kim D-h, Adalsteinsson E, Spielman DM. Simple analytic variable density spiral design. Magnetic resonance in medicine. 2003;vol. 50(no. 1):214–219. [PubMed]
12. Heid O. The fastest circular k-space trajectories. Proceedings of the 10th Annual Meeting of ISMRM; Honolulu. 2002. p. 2364.
13. Simonetti OP, Duerk JL, Chankong V. An optimal design method for magnetic resonance imaging gradient waveforms. IEEE Trans Med Imaging. 1993;vol. 12(no. 2):350–360. [PubMed]
14. Simonetti OP, Duerk JL, Chankong V. MRI gradient waveform design by numerical optimization. Magn Reson Med. 1993;vol. 29(no. 4):498–504. [PubMed]
15. Dale BM, Duerk JL. Time-optimal control of gradients. Proceedings of the 10th Annual Meeting of ISMRM; Honolulu. 2002. p. 2361.
16. Hargreaves BA, Nishimura DG, Conolly SM. Time-optimal multi-dimensional gradient waveform design for rapid imaging. Magn Reson Med. 2004;vol. 51(no. 1):81–92. [PubMed]
17. Bobrow SDJ, Gibson J. Time-optimal control of robotic manipulators along specified paths. The International Journal of Robotics Research. 1985;vol. 4(no. 3):3–17.
18. Kim S-J, Choi D, Ha I. A comparison principle for state-constrained differential inequalities and its application to time-optimal control. IEEE Transactions on Automatic Control. 2005;vol. 50(no. 7):967–983.
19. Shin K, Mckay N. Minimum-time control of robotic manipulators with geometric path constraints. IEEE Transactions on Automatic Control. 1985;vol. 30(no. 6):531–541.
20. Boyce WE, DiPrima RC. Elementary Differential Equations and Boundary Value Problems. 6th ed. John Wiley&Sons, Inc.; 1997.
21. Spielman DM, Pauly JM, Meyer CH. Magnetic resonance fluoroscopy using spirals with variable sampling densities. Magnetic resonance in medicine. 1995;vol. 34(no. 3):388–394. [PubMed]
22. Heberlein K, Hu X. Auto-calibrated parallel spiral imaging. Magnetic resonance in medicine. 2006;vol. 55(no. 3):619–625. [PubMed]
23. Buc EK, Kohler S, Hancu I. 13C multi-spectral 2D rosette imaging. Proceedings of the 15th Annual Meeting of ISMRM; Berlin. 2007. p. 1379.