U.S. patent application number 10/682637 was filed with the patent office on 2005-04-14 for mri gradient waveform design using convex optimization.
This patent application is currently assigned to THE BOARD OF TRUSTEES OF THE LELAND STANDFORD JUNIOR UNIVERSITY. Invention is credited to Conolly, Steven M., Hargreaves, Brian A..
Application Number | 20050077895 10/682637 |
Document ID | / |
Family ID | 34422572 |
Filed Date | 2005-04-14 |
United States Patent
Application |
20050077895 |
Kind Code |
A1 |
Hargreaves, Brian A. ; et
al. |
April 14, 2005 |
MRI gradient waveform design using convex optimization
Abstract
A time-optimal MRI gradient design method utilizes constrained
optimization to design minimum-time gradient waveforms that satisfy
gradient amplitude and slew-rate limitations. Constraints are
expressed as linear equations which are solved using linear
programming, L1-norm formulation, or second-order cone programming
(SOCP).
Inventors: |
Hargreaves, Brian A.; (Menlo
Park, CA) ; Conolly, Steven M.; (Palo Alto,
CA) |
Correspondence
Address: |
BEYER WEAVER & THOMAS LLP
P.O. BOX 70250
OAKLAND
CA
94612-0250
US
|
Assignee: |
THE BOARD OF TRUSTEES OF THE LELAND
STANDFORD JUNIOR UNIVERSITY
Palo Alto
CA
|
Family ID: |
34422572 |
Appl. No.: |
10/682637 |
Filed: |
October 8, 2003 |
Current U.S.
Class: |
324/307 ;
324/309 |
Current CPC
Class: |
G01R 33/561 20130101;
G01R 33/4824 20130101 |
Class at
Publication: |
324/307 ;
324/309 |
International
Class: |
G01V 003/00 |
Goverment Interests
[0001] The U.S. Government has rights in the disclosed invention
pursuant to NIH Grant Nos. NIH-HL39297, NIH-HL56394, NIH-AR46904
and NIH-CA50948 to Stanford University.
Claims
What is claimed is:
1. A method of magnetic gradient waveform design for rapid magnetic
resonance imaging using constrained optimization comprising the
steps of: a) selecting a number N of discrete-time waveforms with a
sampling period, .tau., and time T=N.tau., b) expressing
constraints as equations, c) identifying if a solution exists for
T, and if not, increasing T until a solution exists, d) decreasing
T to find a shortest solution, and e) solving the equations to
identify the solution.
2. The method as defined by claim 1 wherein step e) utilizes linear
programming to identify the solution.
3. The method as defined by claim 2 wherein the linear programming
finds the vector x that minimizes a cost function f.sup.Tx subject
to the constraint Ax=<b, where matrix A and vector b are formed
by combining all of the linear constraint equations for the
amplifier and pulse sequence constraints.
4. The method as defined by claim 1 wherein steps c), d), and e)
utilize L1-norm formulation.
5. The method as defined by claim 4 wherein the gradient is one
dimensional.
6. The method as defined by claim 4 wherein the gradient is two
dimensional.
7. The method as defined by claim 4 wherein the gradient is three
dimensional.
8. The method as defined by claim 4 wherein in step a) the gradient
amplifier has current and voltage limits (I.sub.max and V.sub.max)
that result in the following limits on G(p): 22 G ( t ) I max and L
t G ( t ) + RG ( t ) V max where L, R, and .eta. are the gradient
coil inductance, resistance, and efficiency.
9. The method as defined by claim 1 wherein steps c), d), and e)
include adding slack variables to the optimization.
10. The method as defined by claim 9 wherein step e) solves for
gradient Gx(n) and a set of variables H.sub.x(n) that converge to
.vertline.Gx(n).vertline. and a set of variables S.sub.x(n) that
converge to .vertline..alpha.Gx(n)+.beta.Gx(n+1).vertline..
11. The method as defined by claim 10 wherein the slack variables
are forced to converge by adding linear gradient constraints:
-H.sub.x[n]+G.sub.x[n].ltoreq.0 -H.sub.x[n]-G.sub.x[n].ltoreq.0 and
slew-rates constraints:
-S.sub.x[n]+.alpha.G.sub.x[n]+.beta.G.sub.x[n+1].- ltoreq.0
-S.sub.x[n]-.alpha.G.sub.x[n]-.beta.G.sub.x[n+1].ltoreq.0 whereby
the constraints force the Hx(n) and Sx(n) variables to approach an
appropriate absolute value when combined with minimization of the
following pulse function: 23 J ( S x [ n ] , H x [ n ] ) = n = 1 N
( H x [ n ] + S x [ n ] ) .
12. The method as defined by claim 1 wherein step e) utilizes
second-order cone programming (SOCP) to find a solution x that
minimizes a linear pulse function f.sup.Tx subject to a second
order cone constraint .parallel.Ax+b.parallel..sub.2.ltoreq.Cx+d
which is a superset of the linear constraints in linear programming
where matrix A and vector b are formed by combining all linear
constraint equations for amplifier and pulse sequence
constraints.
13. The method as defined by claim 1 wherein steps c), d) and e)
use a binary-search to minimize time in identifying the minimum
value of n which provides a solution.
14. The method as defined by claim 13 wherein the gradient is one
dimensional.
15. The method as defined by claim 13 wherein the gradient is two
dimensional.
16. The method as defined by claim 13 wherein the gradient is three
dimensional.
17. The method as defined by claim 13 wherein in step a) the
gradient amplifier has current and voltage limits (I.sub.max and
V.sub.max) that result in the following limits on G(p): 24 G ( t )
I max and L t G ( t ) + RG ( t ) V max where L, R, and .eta. are
the gradient coil inductance, resistance, and efficiency.
18. The method as defined by claim 1 wherein the gradient is one
dimensional.
19. The method as defined by claim 1 wherein the gradient is two
dimensional.
20. The method as defined by claim 1 wherein the gradient is three
dimensional.
21. The method as defined by claim 1 wherein in step a) the
gradient amplifier has current and voltage limits (I.sub.max and
V.sub.max) that result in the following limits on G(p): 25 G ( t )
I max and L t G ( t ) + RG ( t ) V max where L, R, and .eta. are
the gradient coil inductance, resistance, and efficiency.
22. The method as defined by claim 1 and further including: f)
establishing physical constraints of a gradient amplifier and a
gradient coil.
23. The method as defined by claim 22 and further including: g)
establishing beginning and end gradient boundary constraints and
fixed area for one or more gradient moments.
24. The method as defined by claim 1 and further including: f)
establishing beginning and end gradient boundary constraints and
fixed area for one or more gradient moments.
25. The method as defined by claim 1 wherein step b) includes
gradient heating constraints.
26. The method as defined by claim 25 wherein step b) includes
gradient magnetostimulation constraints.
27. The method as defined by claim 1 wherein step b) includes
gradient magnetostimulation constraints.
Description
BACKGROUND OF THE INVENTION
[0002] This invention relates generally to magnetic resonance
imaging (MRI), and more particularly, the invention relates to the
design of magnetic gradients for use in MRI.
[0003] Magnetic resonance imaging (MRI) requires placing an object
to be imaged in a static magnetic field, exciting nuclear spins in
the object within the magnetic field, and then detecting signals
emitted by the excited spins as they precess within the magnetic
field. Through use of magnetic gradient and phase encoding of the
excited magnetization, detected signals can be spatially localized
in three dimensions.
[0004] Advances in MR system hardware enable the use of new rapid
pulse sequences. Typical rapid sequences include rapid
gradient-echo (FLASH, GRASS) and multi-echo spin-echo (TSE, FSE,
RARE) sequences. Refocused SSFP (True-FISP, FIESTA, balanced-FFE)
sequences, which produce high signal and contrast, are becoming
common as improved gradients allow imaging with minimal artifacts
from off-resonance. All of these sequences demand efficient
gradient waveform design. Efficient acquisition methods include
echo-planar and spiral imaging. Aside from imaging trajectories,
gradient waveform design includes phase-encoding, prewinder and
rewinder gradients. In rapid sequences with short repetition times,
the design of these latter gradients is an important consideration
in improving imaging efficiency, because their duration reduces the
image acquisition duty cycle.
[0005] In particular, the design problem is to minimize gradient
waveform durations subject to both hardware constraints and
sequence constraints such as desired gradient area. Numerous
previous works have presented different methods to optimize
gradients in different situations. Many of these methods are
limited to the design of trapezoidal pulses, and most have been
demonstrated for 1D gradient design. Simonetti et al. presented a
general technique to minimize 1D-gradient properties such as
moments, diffusion-sensitivity 2 or RMS-current of 1D gradients
using numerical optimization. See "An Optimal Design Method for
Magnetic Resonance Imaging Gradient Waveforms," IEEE Trans Med
Imaging 1993; 12:350-360; and "MRI Gradient Waveform Design by
Numerical Optimization," Magn Reson Med 1993; 29: 498-504.
BRIEF SUMMARY OF THE INVENTION
[0006] The design of rapid imaging sequences requires time-optimal
magnetic gradient waveforms to provide maximum spatial resolution
while keeping both repetition times and scan times low. In
accordance with the invention, constrained optimization is used in
designing gradients that satisfy gradient-amplitude and slew-rate
limitations and additional constraints such as beginning and end
amplitude, gradient area and higher-order gradient moments.
[0007] The design is formulated as linear equations for a number,
N, of discrete-time waveforms with a sampling period .tau.. A
minimum time gradient requires the identification of the smallest
value for N for which a design solution exists. A binary search
technique is employed to minimize the time in identifying the
minimum value of N.tau. which provides a solution.
[0008] The linear equations are then solved using known linear
programming, L1-norm formulation, or second-order cone programming
(SOCP).
[0009] The design can also include constraints on the gradient
waveform to avoid overheating and to limit maximum stimulation
thresholds in the waveform algorithm. It is known that the gradient
amplifier, wires, and coils have several heating constraints. Also,
the time varying magnetic field of a gradient amplifier can induce
electric fields of sufficient magnitude to excite nerves. These
constraints can be addressed in the gradient waveform design.
[0010] The invention and objects and features thereof will be more
readily apparent from the following detailed description and
appended claims when taken with the drawings.
BRIEF DESCRIPTION OF THE DRAWINGS
[0011] FIG. 1 illustrates a voltage model and a
constant-slew-rate-limit model.
[0012] FIG. 2 illustrates discrete approximation of a gradient
waveform.
[0013] FIG. 3 illustrates piece-wise-linear approximation to
quadratic constraints.
[0014] FIG. 4 illustrates piece-wise-linear approximation to a
spherical surface.
[0015] FIGS. 5a-5c illustrate comparison of speed of different
formulations for 1D, 2D and 3D, respectively.
[0016] FIGS. 6a-6d illustrate first-moment-nulled gradient
design.
[0017] FIGS. 7a-7b illustrate echo-planar transition waveform
design for gradient waveforms and slew-rates.
[0018] FIGS. 8a-8d illustrate spiral rewinder gradient design.
[0019] FIGS. 9a-9d illustrate moment-nulled spiral rewinder
gradient design.
[0020] FIGS. 10a-10c illustrate "stack-of-spiral" rewinder gradient
design.
[0021] FIG. 11 illustrates gradient power for different oblique
gradient angles.
[0022] FIGS. 12a-12b illustrate on-axis echo-planar transition
waveform design.
[0023] FIGS. 13a-13b illustrate off-axis echo-planar transition
waveform design.
[0024] FIG. 14 is a graph illustrating voltage and current
constraints and magnetic stimulation constraints for safe operation
of a gradient amplifier.
DETAILED DESCRIPTION OF ILLUSTRATIVE EMBODIMENTS
[0025] The design method can be used in accordance with the
invention to design minimum-time gradient waveforms that satisfy
gradient-amplitude and slew-rate limitations. Additional
constraints such as gradient start and end amplitudes, gradient
area or higher-order gradient moments can be included. Importantly,
application of the method can be extended to multi-dimensional
gradient design which is necessary for procedures that use oblique
scan-plane orientations.
[0026] Consider first the general constraints in gradient waveform
design including gradient amplifier voltage and current limits as
well as pulse-sequence requirements. Next, methods of formulating
the problem are described so that it can be solved using three
different types of existing linear and quadratic programming
techniques. All of these techniques result in comparable
solutions.
[0027] A gradient amplifier supplies current to a gradient coil,
producing a continuously-varying gradient, G(t). The amplifier has
current and voltage limits (I.sub.max and V.sub.max) that result in
the following limits on G(t) as described in King et al.,
"Optimized Waveforms for Spiral Scanning," Magn Reson Med 1995;
34:156-160. 1 G ( t ) I max ( 1 ) and L t G ( t ) + RG ( t ) V max
. ( 2 )
[0028] Here L, R and .eta. are the gradient coil inductance,
resistance and efficiency (i.e., in mT/m/A). In MRI systems, the
gradient can be considered a vector, {right arrow over (G)}(t) that
is composed of three components, G.sub.x(t), G.sub.y(t), and
G.sub.z(t). The primary purpose of the invention is to specify
gradient waveforms that can be freely rotated in three dimensions.
This constrains the maximum gradient or change of gradient in any
direction, and the constraints of Eqs. 1 and 2 become (with
G.sub.max=.eta.I.sub.max) 2 ; G ( T ) r; 2 G max ( 3 ) and ; L t G
( t ) + R G ( t ) r; 2 V max . ( 4 )
[0029] Frequently, a simpler "constant-slew-rate-limit" model is
used, where .eta.RI.sub.max is subtracted from the right side of
Eqs. 2 and 4 and the RG(t) or R{right arrow over (G)}(t) term is
omitted. For a single gradient axis, the constraints of the voltage
model and the constant-slew-rate-limit model are shown in FIG. 1
which illustrates limits on gradient amplitude and slew-rate for a
single-axis gradient due to maximum current (I.sub.max) and maximum
voltage (V.sub.max). Limits are shown using the common
constant-slew-rate-limit model (white area) and the more accurate
voltage-limit model (white plus light gray areas). As the ratio of
coil resistance (R) to coil inductance (L) increases, the
constant-slew-rate-model becomes increasingly constraining, and the
voltage-limit model should be used. .eta. is the ratio of gradient
strength to gradient current. As the ratio of coil resistance to
coil inductance increases, the constant-slew-rate-limit model
becomes more constraining.
[0030] Imaging requirements of the pulse sequence constrain the
k-space change (.DELTA.{right arrow over (k)}) that the gradients
will produce in terms of the gyromagnetic ratio, .gamma. as 3 t = 0
T G ( t ) t = 2 k . ( 5 )
[0031] Additionally, there are cases where higher-order moments of
the gradients can be specified, as 4 t = 0 T t q G ( t ) t = 2 m q
( 6 )
[0032] where .DELTA.{right arrow over (m.sub.q)} is a specified
design constraint on the change in the q.sub.th moment of the
gradient.
[0033] Boundary constraints on the gradient waveforms include
having specified initial and final values for the gradient, i.e.,
{right arrow over (G)}(O)={right arrow over (g.sub.0)} and {right
arrow over (G)}(T)={right arrow over (g.sub.f)}. This is required
to design gradient waveforms at the start and end of the sequence,
and to link different waveforms together.
[0034] Practical gradients are generated as discrete-time waveforms
with a sampling period .tau.. A gradient waveform can be
represented as a discrete sequence {right arrow over (G)}[n]={right
arrow over (G)}(n.tau.) with n=1 . . . N. FIG. 2 shows a single
component of the discrete-time vector {right arrow over (G)}[n].
The constraints in Eqs. 1-6 can easily be converted to
discrete-time equivalents, where a derivative is approximated by
the slope between two consecutive samples, and an integral is
approximated by a discrete sum multiplied by .tau.. Specifically,
it is useful to express the voltage constraint (Eq. 4) as 5 ; G [ n
] + G [ n + 1 ] r; 2 ( 7 ) where = R - L , = L , and = V max
[0035] for the voltage model. In the constant-slew-rate-limit
model, 6 = - L , = L , and = V max - RI max .
[0036] Typically the gradients, as well as the desired moments,
consist of two or three separate gradient axes. In basic cases,
these axes can be treated separately and all constraints expressed
as described in the previous section. However, there are numerous
cases where gradients must be designed in a 2D or 3D space that
must be free to rotate.
[0037] First we now define the solution space consisting of three
separate discrete-time waveforms along each axis, G.sub.x[n],
G.sub.y[n] and G.sub.z[n], which must collectively satisfy the
quadratic constraints of Eqs. 3 and 4. The boundary constraints and
moment constraints (Eqs. 5 and 6) are applied separately for each
of the dimensions. Overall, the problem is now quadratically
constrained.
[0038] Before formulating solutions to the gradient design problem,
it is useful to summarize the constraints for discrete time
waveform design. Table 1 lists the four types of constraint for
gradient design problems for D-dimensional, N-point waveform with
Q-moments constrained, and the number of instances of each in
gradient design. The quadratic inequality constraints are the prime
source of complexity in gradient design.
1TABLE 1 Constraint Type Number Current (Gradient Amplitude)
Quadratic Inequality N Voltage (Gradient Slew-Rate) Quadratic
Inequality N Gradient Moments Linear Equality D .times. Q Boundary
Values Linear Equality 2D
[0039] A certain class of constrained optimization problems,
convex-optimization problems, have been studied heretofore. For
many types of convex optimization problems, very efficient and
robust solution methods have been developed. These methods are
guaranteed to find the globally-optimal solution if any solution
exists. Following are described methods of expressing the gradient
design constraints as standard convex optimization problems so that
they may be solved efficiently and reliably.
[0040] Both absolute-value and equality constraints can
equivalently be expressed as two inequality constraints. A single
absolute-value constraint a.vertline.x.vertline..ltoreq.b can be
replaced by two linear constraints, ax.ltoreq.b and -ax.ltoreq.b.
An equality constraint of the form ax=b can be replaced by the two
inequalities ax.ltoreq.b+.epsilon. and -ax.ltoreq.-b+.epsilon..
This is desirable in practical optimization, as it allows the
tolerance .epsilon. to be different for different constraints.
[0041] Using this formulation, the constraints of Eqs. 1, 2, 5 and
6 as well as the boundary constraints can all be expressed as
linear inequality constraints on the discrete-time gradient
waveform.
[0042] Each quadratic constraint may be approximated by multiple
linear constraints. This is done in 2D, for example, by
approximating the circle
G.sub.x[n].sup.2+G.sub.y[n].sup.2=G.sub.max.sup.2 as a series of P
lines defined by 7 cos ( p ) G x [ n ] + sin ( p ) G y [ n ] G max
cos ( P ) = 1 P ( 8 ) where p = 2 p - 0.5 P .
[0043] FIG. 3 shows this piece-wise-linear approximation to a 2D
quadratic constraint with P=8. The white area represents allowable
values of the 2D amplitude, voltage, or slew-rate. For all of the
2D examples herein, P=16 is used. A piece-wise-linear approximation
to a sphere is also possible for 3D quadratically constrained
problems. A simple way to do this is to consider the surface of the
sphere .vertline.{right arrow over (r)}.vertline.=c in the first
octant where the components of {right arrow over (r)} are all
positive (FIG. 4). By placing points on the surface of the sphere,
triangles with roughly equal area are constructed, and the
constraints can be expressed in the form
{right arrow over (n)}.multidot.{right arrow over
(r)}.ltoreq.{right arrow over (n)}.multidot.{right arrow over (V)}
(9)
[0044] where {right arrow over (n)} is a vector normal to the
triangle, i.e., {right arrow over (n)}=({right arrow over
(B)}-{right arrow over (A)})-.times.({right arrow over (C)}-{right
arrow over (A)}) and {right arrow over (V)} is the position vector
of any of the triangle corners, ({right arrow over (A)}, {right
arrow over (B)}, or {right arrow over (C)}).
[0045] For the 3D examples herein, the first quadrant is divided
into 9 planes using the points shown in FIG. 4. Here
piece-wise-linear approximation to spherical surface represents 3D
quadratic constraints. The points shown are permutations of
<1,0,0>, 8 3 2 , 1 2 , 0 , and 1 3 , 1 3 , 1 3 .
[0046] The entire sphere can then be represented by the 72 planes
generated with all permutations of positive and negative signs in
the individual {right arrow over (n)} vectors.
[0047] The problem with only linear inequality constraints can be
solved using linear-programming (LP). The "standard LP" form of a
problem is to find the vector x that minimizes a cost function
f.sup.Tx subject to the constraint Ax.ltoreq.b. The matrix A and
vector b are formed by combining all of the linear constraint
equations for amplifier and pulse sequence constraints. This
formulation of the problem is referred to as "simple-LP"
formulation.
[0048] The cost function in simple-LP formulation is linear in the
gradient values. Minimization of such a cost function is not
particularly useful. As such the actual optimization consists of a
series of feasibility tests as will be described later. For the
simple-LP formulation, f is chosen to be a vector of ones.
[0049] The LP formulation for an N-point waveform uses ND variables
with roughly 2.sup.D+1NP constraints, where D is the number of
dimensions and P is the number of piece-wise-linear constraints on
gradient amplitude and voltage amplitude for which coefficients are
all positive (i.e., the number of constraints in the first quadrant
or octant). The Appendix shows the complete simple-LP formulation
for a 2D freely-rotatable gradient design problem.
[0050] An alternative to simple-LP formulation is the L1-norm
formulation, which alters the number of constraints and can produce
different solutions to the problem, see Xu et al., "Homogeneous
Magnet Design Using Linear Programming," IEEE Trans Med Imaging
2000; 36:476-483.
[0051] The L1-norm formulation uses standard linear programming to
minimize the L1-norm (absolute value) of the gradient and slew
voltage by adding "slack variables" to the optimization problem. In
addition to solving for G.sub.x[n], we solve a set of variables
H.sub.x[n], that converge to .vertline.G.sub.x[n].vertline. and a
set of variables S.sub.x[n] that converge to
.vertline..alpha.G.sub.x[n]+.beta.G.sub.x[n+1- ].vertline..
[0052] The slack variables are forced to converge by first adding
the following linear gradient constraints to the problem:
-H.sub.x[n]+G.sub.x[n].ltoreq.0
-H.sub.x[n]-G.sub.x[n].ltoreq.0 (10).
[0053] Additionally we add the following slew-rate constraints:
-S.sub.x[n]+.alpha.G.sub.x[n]+.beta.G.sub.x[n+1].ltoreq.0
-S.sub.x[n]-.alpha.G.sub.x[n]-.beta.G.sub.x[n+1].ltoreq.0 (11).
[0054] These constraints will force the H.sub.x[n] and S.sub.x[n]
variables to approach the appropriate absolute value when combined
with minimization of the following cost function: 9 J ( S x [ n ] ,
H x [ n ] ) = n = 1 N ( H x [ n ] + S x [ n ] ) . ( 12 )
[0055] For two and three dimensional problems, an appropriate set
of variables H.sub.y[n], S.sub.y[n], H.sub.z[n] and S.sub.z[n]
would be added. Constraints for these variables are expressed as in
Eqs. 10 and 11, and the variables would be added to the cost
function in Eq. 12.
[0056] The simple-LP cost function of all ones simply tries to make
individual gradient axes closer to their negative amplitude limits.
However, the L1-norm cost function is preferable, as it tends to
minimize the magnitude of both gradient current and voltage, which
will improve heating characteristics of the waveform.
[0057] When L1-norm formulation is used, the piece-wise-linear
constraints on gradient amplitude may be reduced by replacing the
constraint on G.sub.x[n], G.sub.y[n], and G.sub.z[n] with
equivalent constraints on H.sub.x[n], H.sub.y[n], and H.sub.z[n].
Similarly, the voltage constraints may be replaced with constraints
on S.sub.x[n], S.sub.y[n], and S.sub.z[n]. Since all of the slack
variables are constrained to be non-negative, only the gradient
amplitude and voltage-limit constraints with positive coefficients
are needed. This reduces the number of gradient amplitude
constraints by a factor of 2.sup.D where D is the number of
gradient axes, and is the primary advantage of the L1-norm
formulation compared to simple-LP formulation.
[0058] The L1-norm formulation uses 3ND variables, and
approximately 2N(P+2D) constraints. For 2D or 3D problems, the
number of constraints is often much smaller than that of the
simple-LP formulation. Additionally, as described by Xu et al.,
L1-norm programming tends to find solutions where some or many of
the variables are zero. This may be advantageous in gradient
design, as solutions with lower gradient or slew duty cycle result
in lower gradient heating.
[0059] Second-order cone programming (SOCP) (see Lobo et al.,
"Applications of second-order cone programming," Linear Algebra
Appl. 1998; 284: 193-228) is a method that finds the solution x
that minimizes a linear cost function f.sup.Tx subject to the
"second-order cone" constraint.
.parallel.Ax+b.mu..sub.2.ltoreq.Cx+d (13).
[0060] First this constraint is a superset of the linear
constraints in linear programming. Second, the cost function is
linear, as with the simple-LP formulation described above, and is
chosen in the same manner, as all ones.
[0061] The SOCP formulation of the gradient design problem begins
with Eqs. 5-6. Although these constraints could be reduced by a
factor of 2.sup.D compared with the simple LP or L1-norm
formulations, this has little effect on convergence time. The main
advantage of the SOCP formulation is that the quadratic limits of
Eqs. 3 and 4 result in a total of 2N-1 constraints, regardless of
the number of gradient dimensions. Additionally, the SOCP
formulation does not approximate the quadratic gradient limits,
which sometimes results in slightly shorter duration solutions.
[0062] Because the time to find a solution increases with the
number of constraints, it is useful to try to further reduce the
number of constraints. Assuming that the slew-rate constraints are
met, there will be gradient amplitude constraints near the
endpoints that are redundant. Specifically, if the (minimum)
distance from the start or end gradient to a linear or quadratic
gradient constraint is .DELTA.G, then that constraint is redundant
for at least a time 10 L G V max + RG max .
[0063] The minimum-time gradient design problem is to find
N.sub.min, the minimum value of N for which a solution exists.
Linear programming functions can be called sequentially with
different values of N. We first assume that either {right arrow
over (G)}(0)=0 or {right arrow over (G)}(n.tau.)=0. If a solution
exists for N, then N.sub.min.ltoreq.N. Otherwise,
N.sub.min.gtoreq.N.
[0064] This suggests using a binary-search technique. This
technique divides the unknown interval containing N.sub.min on each
call. Begin with estimates of the upper and lower bounds N.sub.u
and N.sub.l on N. First these are tested to ensure that they are in
fact upper and lower bounds. If not, both bounds are raised or
lowered accordingly and retested. Once N.sub.u and N.sub.l have
been verified, N is set as close as possible to
N.sub.l+v(N.sub.u-N.sub.l). Depending on whether or not a solution
exists, either N.sub.u or N.sub.l is set to N, and the process is
repeated. Use v=0.8 because an existing solution is found
considerably more quickly than the lack of a solution using
SOCP.
[0065] The above formulation has been implemented using Matlab 6.0
with functions from the Mathworks Optimization toolbox (Mathworks,
Natick, Mass.). Specifically, use the linprog( ) function, which
replaces the old lp( ) function in Matlab, as well as the socp( )
function described above. Matlab functions are available for
general use.
[0066] First, compare the time required to solve the minimum time
gradient-area-nulling problem using the simple-LP, L1-norm and SOCP
methods. All tests used Matlab 6.0 on 1.8 GHz Athlon PCs running
Linux. The binary search formulation was identical for all three
methods, and the comparison was performed for 1D, 2D and 3D
problems.
[0067] Generated were 1000 random gradient waveforms of random
duration and the given dimension with endpoints within the
magnitude constraint for the given solution method. A sample period
of .tau.=20 .mu.s was used, which assumes a gradient amplifier
bandwidth of less than 25 kHz. The gradient and zeroth moment of
the gradient were then rewound to zero using the bisection method
described above. Each step of the binary search method was solved
using simple-LP, L1-norm or SOCP formulation. The times required to
solve the rewinder problem for each waveform and each formulation
method were recorded. In all cases a maximum gradient amplitude of
40 mT/m, and a constant-slew-rate model with 150 T/m/s maximum
slew-rate were used. For 2D and 3D, the rewinders are limited by
the quadratic constraints described above.
[0068] FIGS. 5a-5c compare convergence times for the three methods
for 1D, 2D and 3D design problems. Shown are convergence time for
1D (a), 2D (b), and 3D (c) time-optimal gradient design using
simple-LP (dotted line), L1-norm (dashed line) and SOCP (solid
line). Although simple-LP formulation is the most efficient for 1D
design problems, the L1-norm is significantly better for 2D and 3D
problems. SOCP outperforms L1-norm in 2D and 3D problems, except
when the solution turns out to be very long. Not surprisingly,
convergence time increases with the solution length. SOCP
formulation is fastest regardless of dimension for most practical
waveforms (i.e., duration less than 1 ms), and is significantly
faster than the other methods as the number of gradient dimensions
increases to 2 or 3.
[0069] Presented now are several examples showing applications of
linear programming for time optimal gradient design. We begin with
1D examples that show that linear programming produces the
time-optimal waveforms. Next we show 2D examples, particularly for
use in spiral imaging. Finally we show 3D examples for
double-oblique or oblique spiral scans.
[0070] For a constant amplitude limit, constant-slew-rate limit
gradient, the minimum time waveforms to produce a given zeroth,
first and/or second moment are generally triangular or trapezoidal
waveforms. However, calculation of the segments sometimes must be
solved numerically.
[0071] A first example is the generation of a simple 1D
moment-nulled phase-encoding gradient. The boundary constraints on
the pulse waveform are G.sub.x[0]=G.sub.x[N]=0 and the desired
k-space area and first-moments are .DELTA.{right arrow over
(k)}=k.sub.des and m.sub.1=0. This example sets k.sub.des=400
m.sup.-1, G.sub.max=40 mT/m and maximum slew-rate 11 L = 150 T / m
/ s .
[0072] The resulting gradient waveform, k-space trajectory,
slew-rate and first moment are shown in FIG. 6 which shows (a)
first-moment-nulled gradient waveform, (b) corresponding k-space
trajectory, (c) gradient slew-rate, and (d) gradient first-moment,
all plotted as functions of time for Example A. The gradient
waveform has a duration of 1.048 ms. The gradient is a bipolar
combination of a trapezoidal lobe with a triangular lobe, as can
reasonably be expected given amplitude- and slew-limited
constraints in one dimension. Note that either the gradient
amplitude limit or the slew-rate limit is always met, which should
be the case for a minimum-time gradient in one dimension.
[0073] A second example is to design a minimum-time waveform that
transitions between successive lines of an echo-planar imaging
trajectory. For consistency, the same gradient and slew-limits, 20
mT/m and 200 T/m/s respectively, are used with the assumption that
the gradients can be freely rotated to any oblique in-plane
position. The EPI trajectory uses a 20 mT/m readout gradient, with
0.8 cm.sup.-1 phase-encode area between lines. In our formulation,
we specify G.sub.x[1]=-G.sub.x[N]=20 mT/m, G.sub.y[1]=G.sub.y[N]=0
mT/m, .DELTA.k.sub.x=0 cm.sup.-1 and .DELTA.k.sub.y=0.8 cm.sup.-1,
.tau.-1 .mu.s.
[0074] The result, shown in FIG. 7, has a total duration of 245
.mu.s. In FIG. 7, gradient waveforms (a) and slew-rates (b) for an
echo-planar readout/phase-encode transition are shown. Gradients
and slew-rates are shown for the x-axis (solid line) y-axis (dashed
line) and magnitude (dotted line). This method performs as well as
the time-optimal method described in ref. (9).
[0075] The solutions with both methods look very similar, and the
minor difference between them is most likely simply the way that
the end-point constraint is formulated. Interestingly, this
gradient is not zero at either endpoint. However, because of the
symmetry of the problem, the bisection method of finding the
minimum-time gradient still works.
[0076] Minimum gradient waveform design on one dimension can
generally be solved exactly and analytically. However,
two-dimensional design is not always simple. Consider spiral
imaging, where interleaved spiral waveforms are played by rotating
the waveform in the x-y gradient plane. For the trajectory to be
played at arbitrary rotation angles, the gradient and slew-rate
vector magnitudes must be within slew-rate constraints.
[0077] Spiral waveform design methods have been addressed
previously. The spiral waveform can be rewound by separately
rewinding each axis with both amplitude (current) and voltage
limits set to 1/{square root}{square root over (2)} times the
individual axis limits. This is adequate for applications where the
spiral readout duration is long, but is not time-optimal. In cases
such as spiral steady-state free precession imaging, it is
desirable to rewind the first moment of the gradient as well. In
this example, a time-optimal rewinder for spiral gradient waveforms
as well as a rewinder for a first-moment-nulled spiral waveform are
used.
[0078] For this example, a spiral waveform can be used similar to
that used for SSFP imaging. The design assumes 60 spiral
interleaves to achieve a resolution of 1.25 mm and a field-of-view
of 20 cm with a readout duration of 2.7 ms.
[0079] The optimization goal is to rewind the gradient waveforms
and zeroth moment vector simultaneously to zero as quickly as
possible. The rewinder gradient vector initially equals the final
gradient of the readout, and ends at zero in both dimensions.
Combined with the constraint of rewinding the zeroth moment vector,
this gives six equality constraints. FIG. 8 shows the spiral
gradient waveform and rewinder characteristics. FIGS. 8a-8d
illustrate gradient waveforms (a), k-space trajectory (b),
slew-rates (c) and first moment (d) for a simple rewinder waveform
(gray region) for a 2D-spiral imaging gradient. Each plot shows x
(solid line) and y (dashed line) components, as well as the
magnitude (dotted line). The k-space-only rewinder lasts for 0.66
ms, and has a significant residual first moment. All quantities
reach zero at the end of the waveform. As with the 1D example,
either the gradient magnitude or slew-rate magnitude limit is
reached for the entire duration of the rewinder.
[0080] The spiral rewind problem can be repeated for the case where
the first moments are also rewound at the end of the waveform. This
adds a constraint in each axis for a total of eight equality
constraints. The L1-norm solution is shown in FIGS. 9a-9d. Here are
shown gradient waveforms (a), k-space trajectory (b), slew-rates
(c) and first moment (d) for an m.sub.1-nulled rewinder waveform
(gray region) for a 2D-spiral imaging gradient. Each plot shows x
(solid line) and y (dashed line) components, as well as the
magnitude (dotted line). The m.sub.1-nulled rewinder lasts for 1.23
ms, and rewinds the gradient, k-space location and the first
moment. For nulling of the first moments, the minimum-time waveform
is significantly increased, from 0.78 ms to 1.38 ms. If the two
axes were rewound separately (as is commonly performed), each using
70% of the maximum power, the duration would be closer to 1.8 ms.
In rapid imaging sequences such as SSFP, 0.4 ms is a significant
time reduction, i.e., 10% of the repetition time or 25% of the data
acquisition time.
[0081] There are cases where the gradient vector must be rotated
freely in three dimensions, such as double-oblique scans or
oblique-plane spiral scans. The 3D "stack of spirals" example in
FIGS. 10a-10c pertains to the latter. Shown are gradient waveforms
(a), k-space trajectory (b) and slew-rates (c) for a 3D rewinder
waveform (gray region) for a "stack-of-spirals" imaging gradient.
Each plot shows x (solid line), y (dashed line) and z (dash-dot
line) components, as well as the magnitude (dotted line). The
rewinder gradient lasts for 0.74 ms. For gradient-spoiled
acquisition, it is desirable to play a gradient spoiler on the
z-axis following the spiral readout. To minimize the overall time,
this spoiler should be played out simultaneously with the spiral
rewinder.
[0082] The problem is formulated the same way as that of the simple
spiral-rewinder in Example B with the additional constraints that
in the z-axis, we specify G.sub.x[1]=G.sub.x[N]=0 and a desired
k-space area.
[0083] The LP formulation is useful for optimizing oblique gradient
design. Gradient design is simplest when the three gradient axes
may be designed independently of each other, that is when each
"logical" gradient axis corresponds exactly to one "physical"
gradient axis. This non-oblique case corresponds to the 3D LP
formulation where the gradient constraints are simply
.vertline.G.sub.x[n].vertline..ltoreq.G.sub.max,
.vertline.G.sub.y[n].vertline..ltoreq.G.sub.max, and
.vertline.G.sub.z[n].vertline..ltoreq.G.sub.max.
[0084] Gradient limits can be reset for any rotation angle by
scaling the logical gradient limit by (cos .phi.+sin .phi.).sup.-1,
where .phi. is the angle between major logical and physical axes.
As shown in FIG. 11, this ensures that the gradient vector always
lies within the large gray square. However, this is obviously
inefficient in that a significant amount of gradient power is not
used. FIG. 11 illustrates gradient amplitude usage for different
oblique-scan strategies. The dark gray square represents all usable
gradient amplitude, and the circle represents the
quadratically-constrained amplitude if the gradient is designed to
be freely-rotatable. If gradients are designed completely
independently, then the amplitude limits on each axis must be
reduced so that the worst-case gradient still lies within the dark
gray square. This is shown for a rotation angle of
.theta.=15.degree. (dashed lines). With constraints tailored to the
specific oblique angle, all gradient power is available for any
rotation angle.
[0085] The LP formulation is very adaptable, as the physical
constraints can simply be rotated to the logical coordinate system.
This means that for oblique gradient design all available gradient
power is always used.
[0086] Consider the echo-planar transition of example B where the
logical and physical gradients are rotated by .phi.=0.degree. and
.phi.=30.degree.. Here the design is gradient for each specific
oblique position. The gradient amplitude constraints are
cos .phi.G.sub.x[n]-sin .phi.G.sub.y[n].ltoreq.G.sub.max
-cos .phi.G.sub.x[n]+sin .phi.G.sub.y[n].ltoreq.G.sub.max
sin .phi.G.sub.x[n]+cos .phi.G.sub.y[n].ltoreq.G.sub.max
-sin .phi.G.sub.x[n]-cos .phi.G.sub.y[n].ltoreq.G.sub.max (14)
[0087] where .phi. is the angle of rotation between logical and
physical coordinate systems. The slew-rate constraints are
cos .phi.(.alpha.G.sub.x[n]+.beta.G.sub.x[n+1])-sin
.phi.(.alpha.G.sub.y[n]+.beta.G.sub.y[n+1]).ltoreq..psi.
-cos .phi.(.alpha.G.sub.x[n]+.beta.G.sub.x[n+1])+sin
.phi.(.alpha.G.sub.y[n]+.beta.G.sub.y[n+1]).ltoreq..psi.
sin .phi.(.alpha.G.sub.x[n]+.beta.G.sub.x[n+1])+cos
.phi.(.alpha.G.sub.y[n]+.beta.G.sub.y[n+1]).ltoreq..psi.
-sin .phi.(.alpha.G.sub.x[n]+.beta.G.sub.x[n+1])-cos
.phi.(.alpha.G.sub.y[n]+.beta.G.sub.y[n+1]).ltoreq..psi. (15).
[0088] FIGS. 12 and 13 show the resulting waveforms. Both waveforms
are faster than that shown in FIG. 7, as more gradient power is
available, both for amplitude and slew-rate. The off-axis case
results in a longer waveform, because slewing in two logical axes
must be performed by one physical axis.
[0089] The gradient waveform design can also include constraints to
avoid overheating and to limit maximum thresholds. The gradient
amplifier, wires, and coils have several heating constraints. These
can be expressed as functions of the gradient waveforms and the
derivatives of the waveforms. Often these constraints are binding
in the operation of the scanner. Hence it would be important to be
able to constrain the heating in the waveform design algorithm.
Fortunately it is possible to include constraints on the gradient
waveform and its derivative in the methods disclosed herein.
[0090] The heat dissipated in the coil itself is 12 E coil = 0 T I
2 R coil t ( 16 )
[0091] where R.sub.coil is a constant, and I is proportional to the
gradient waveform, G(t).
[0092] Using the SOCP formulation, the waveforms G(t) for each axis
are approximated by a discrete-time waveform G.sub.i, where i=1 . .
. N. Using this, the quadratic constraint can be expressed directly
and added to the other constraints in the problem, as 13 i = 1 N TG
i 2 < E max ( 17 )
[0093] where E.sub.max is the maximum heat that can be dissipated
in the coil for a given time duration, .eta. is the coil efficiency
(G/cm/A), and T is the gradient sampling step. Alternatively,
instead of constraining the heat, the heat can be minimized by
first adding a slack variable, E, in the SOCP formulation as
follows: 14 i = 1 N TG i 2 < E . ( 18 )
[0094] Then E can be minimized in the SOCP cost function
directly.
[0095] Amplifier heating models vary from system to system, but
generally are of the following form: 15 E amp = 0 T I 2 R sat t + 0
T I V sat t + V pwm I PWM ( I ) . ( 19 )
[0096] The first two terms of this equation can be constrained or
minimized by adding slack variables in the SOCP formulation. The
first term uses one slack variable, E.sub.1, as 16 i = 1 N TR sat G
i 2 2 < E 1 . ( 20 )
[0097] The second term first defines N variables, P.sub.i, as 17 V
sat 2 G i 2 2 < P i . ( 21 )
[0098] Then, neglecting the third term, E.sub.amp can be
constrained or minimized by constraining or minimizing the sum: 18
J = i = 1 N TP i + E 1 . ( 22 )
[0099] The time varying magnetic field of a gradient amplifier has
been known to induce electric fields of sufficient magnitude to
excite nerves. This is known as magnetostimulation or peripheral
nerve stimulation (PNS). Stimulation is a function of gradient
amplitude and the derivative of the gradient amplitude. The FDA
does not advise MRI manufacturers to operate the gradient
amplifiers above stimulation thresholds. Hence, it would be a
contribution to the state of the art of gradient waveform design to
be able to include maximum stimulation thresholds in the waveform
algorithm. Average population studies have determined stimulation
thresholds that are a linear combination of G(t) and dG/dt(t).
[0100] Mathematically, this stimulation limit can be written as 19
G ( t ) + G t < A stim ( 23 )
[0101] where G(t) is the gradient and dG/dt is the first derivative
of the gradient. The stimulation threshold, A.sub.stim, is
determined experimentally, and may vary from patient to patient.
Given A.sub.stim, the stimulation limits can be imposed on the
individual or combined gradient axes in the exact same way as the
amplifier voltage limits. FIG. 14 illustrates the voltage and
current constraints and magnetic stimulation constraints for safe
operation of a gradient amplifier.
[0102] The design of many minimum-time gradient waveforms can be
expressed as a constrained optimization problem. With some
manipulation, the problem can be posed as a standard linear
programming (LP) problem, for which many efficient solving methods
have been developed. The formulation described solves minimum-time
gradient problems where the initial and final gradient values and
desired moments are specified, and where the gradient either starts
or ends at zero.
[0103] Proof of time-optimality is an interesting question. For 1D
cases, this can be done, and results indicate this algorithm
produces the minimum-time solution. For more complex cases where
there is not an analytic solution, it is not possible to prove
time-optimality of the solution. Based on the principles of convex
optimization techniques, it is assumed that the solutions are
time-optimal. Additionally, the simple LP formulation and L1-norm
formulation always give the same duration gradient.
[0104] The minimum-time gradient problem, as expressed above, is a
series of feasibility problems using different length solution
vectors. Shown are three different problem formulations. First, the
problem can be solved simply as a standard LP problem, using
G.sub.x[n], G.sub.y[n], and G.sub.z[n] variables representing the
gradient on each axis. Second, the L1-norm formulation adds
additional slack variables, H.sub.x[n], H.sub.y[n], and H.sub.z[n],
representing the gradient magnitudes, and S.sub.x[n], S.sub.y[n],
and S.sub.z[n] representing the voltage magnitude to further reduce
the number of constraints. The L1-norm formulation also allows a
more physically-meaningful cost function to be expressed. Third,
the second-order cone programming (SOCP) formulation uses
G.sub.x[n], G.sub.y[n], and G.sub.z[n] variables and quadratic
constraints on both gradient amplitude and gradient voltage.
[0105] The simple LP formulation is advantageous compared to
L1-norm only for 1D problems. For 2D or 3D problems for which the
solutions consist of 75 or fewer time samples (1.5 ms with
.DELTA.T=20 .mu.s), the SOCP formulation is the fastest
formulation. Given that the SOCP formulation also uses the
available gradient limits slightly more efficiently than either
simple LP or L1-norm formulation, it is clearly the best option.
However, if SOCP software is not available, the L1-norm formulation
also performs very well for 2D or 3D problems.
[0106] Minimum time, multi-dimensional gradient design is an
increasingly important part of designing MR imaging sequences. The
invention provides a method that can deliver the time-optimal
gradients given a variety of constraints such as amplitude,
slew-rate, end points or gradient moments. The invention uses
standard linear programming tools, which are robust and
efficient.
[0107] While the invention has been described with reference to
specific embodiments, the description is illustrative of the
invention and is not to be construed as limiting the invention.
Various modifications and applications may occur to those skilled
in the art without departing from the true spirit and scope of the
invention as defined by the appended claims.
[0108] While the invention has been described with reference to
specific embodiments, the description is illustrative of the
invention and is not to be construed as limiting the invention.
Various modifications and applications may occur to those skilled
in the art without departing from the true spirit and scope of the
invention as defined by the appended claims.
Appendix A: Simple-LP Gradient Design Formulation
[0109] To actually solve the gradient design problem using linear
programming or SOCP, the matrix A and vector b are constructed
using all of the constraints discussed in the text. For simple-LP
formulation, the constraint matrix form is as follows: 20 ( 1 0 0 0
0 0 0 0 - 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 - 1 0 0 0 0 0 0 1
0 0 0 0 0 0 0 - 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 - 1 C 1 0 0
0 S 1 0 0 0 C P 0 0 0 S P 0 0 0 0 C 1 0 0 0 S 1 0 0 0 C P 0 0 0 S P
0 0 0 0 0 C 1 0 0 0 S 1 0 0 0 C P 0 0 S P C 1 C 1 0 0 S 1 S 1 0 0 C
P C P 0 0 S P S P 0 0 0 0 C 1 C 1 0 0 S 1 S 1 0 0 C P C P 0 0 S P S
P / 2 / 2 / 2 / 2 0 0 0 0 - / 2 - / 2 - / 2 - / 2 0 0 0 0 0 0 0 0 /
2 / 2 / 2 / 2 0 0 0 0 - / 2 - / 2 - / 2 - / 2 ) ( G x [ 1 ] G x [ 2
] G x [ N - 1 ] G x [ N ] G y [ 1 ] G y [ 2 ] G y [ N - 1 ] G y [ N
] ) ( G max ' G max ' G max ' G max ' G max ' G max ' ' ' ' ' k x +
- k x + k y + - k y + )
[0110] Horizontal lines separate four different types of
constraints: boundary constraints, current constraints, voltage
constraints and area constraints. Here C.sub.p=cos(.theta..sub.p)
and S.sub.p=sin(.theta..sub.- p) with .theta..sub.p defined as in
Eq. 8. Also, 21 G max ' = G max cos ( P ) and ' = cos ( P ) .
[0111] In the above expression, E can have different values for
different constraints. (As described in the text, .gamma. is the
gyromagnetic ratio, .tau. is the sampling period, and .psi. is
defined in Eq. 7.
* * * * *