Imagine you are trying to draw a complex picture using a paint brush strapped to a remote-controlled car. To finish the picture as fast as possible, you want to floor the accelerator, but if you go too fast, your car tips over and the paint splatters everywhere.
In the world of MRI, we do something remarkably similar. Instead of a car, we have magnetic gradients, and instead of a physical canvas, we have k-space—a digital map where all our raw image data lives. To “draw” our image, we have to steer a single point of focus across this map as quickly as possible.
Following this analogy, if MRI is the car, then k-space is the canvas, and we leave a trail of data points \(k(t)\) as we drive, if you floor the gas pedal (the Gradient \(G(t)\) your position on the canvas changes faster, and you accumulate all that speed over time:
\[k(t) = \gamma \int_0^t G(\tau) d\tau\]
The most frantic parts of this “drive” are the prephasors and spoilers. Think of these as the high-speed sprints needed to get the car from the “Pit Lane” to the starting line, or the emergency braking required to park it after the race. It turns out that we can use some clever math to optimize these sprints and brakes, ensuring we get the best possible image quality without crashing the car. In this post, we’ll explore how linear and quadratic programming can help us navigate k-space more efficiently.
MRI constraints#
In MRI, we have physical constraints on how fast we can change our gradients (acceleration) and how strong they can be (speed). These constraints are like the speed limits and handling capabilities of our car. If we exceed these limits, we risk “crashing” our MRI system, which can lead to poor image quality or even damage the MR System, there also some concern about Peripheral Nerve Stimulation. In practical terms, we have:
- Maximum gradient strength \(G_{\text{max}}\)
- Maximum slew rate \(S_{\text{max}}\)
- The gradient raster time \(\Delta T\) which is the minimum time interval between changes in the gradient.
Then there are the boundary conditions: We need to reach a specific point in k-space, and be at a given speed at the beginning and end of our trajectory, that is \(G(0) = g_{\text{start}}\) and \(G(N \Delta T) = g_{\text{end}}\)
Since the position in k-space is the integral of the gradient, we can express our constraints in terms of the gradient values at discrete time points. For example, if we have \(N\) time points, we can write:
\[\begin{array}{rl} k(N \Delta T) &= \gamma \Delta T \sum_{n=0}^{N-1} G(n \Delta T)\\ \forall n, \,\,|G(n \Delta T)| &\leq G_{\text{max}} \\ \forall n, \,\,|G(n \Delta T) - G((n-1) \Delta T)| &\leq S_{\text{max}} \Delta T \\ G(0) &= g_{\text{start}} \\ G((N-1) \Delta T) &= g_{\text{end}} \end{array}\]
- To be more precise the area of k-space should probably be computed with a trapezoidal rule.
The textbook solution for reaching a target k-space location is the trapezoidal waveform. It consists of a ramp-up at maximum slew rate, a plateau at maximum gradient strength, and a ramp-down. While this is simple to implement and often “good enough,” it is not necessarily optimal when you have strict boundary conditions \(G_{\text{start}}\neq 0, G_{\text{end}}\neq 0\) or when the duration N is fixed. A rigid trapezoid might overshoot the target or fail to utilize the full available gradient area efficiently, especially when we transition between two active segments of a sequence.
Framing as a linear programming Problem#
The formulation above only contains linear combination of the gradients values at times \(n \Delta T\), and linear constraints on these values, so we can use linear programming to find the optimal trajectory.
We can stack our discrete gradient values into a vector \(\bm{g}=[G_0, G_1, \dots, G_{N-1}]^\top\). The k-space constraint becomes a linear equality \(\bm{A}_{eq} \bm{g} = \bm{b}_{eq}\) With
\[A_{eq} = \begin{bmatrix} 1 & 1 & \dots & 1 \\ 1 & 0 & \dots & 0 \\ 0 & 0 & \dots & 1 \\ \end{bmatrix}, \quad b_{eq} = \begin{bmatrix} \Delta k \\ g_{\text{start}} \\ g_{\text{end}}\end{bmatrix}\]
And the inequality constraint \(A_{ub} g \leq b_{ub}\), we split this in too \( G_{n+1} - G_n \leq S_{\text{max}}\) and \(G_{n} - G_{n+1} \leq S_{\text{max}}\) :
\[A_{ub} = \left[\begin{array}{rrrr} -1 & 1 & 0 & \dots \\ 0 & -1 & 1 & \dots \\ \hline 1 & -1 & 0 &\dots \\ 0 & 1 & -1 & \dots \end{array}\right], \quad b_{ub} = \begin{bmatrix} S_{\text{max}}\Delta T \\ \vdots \\ S_{\text{max}}\Delta T\end{bmatrix}\] Combined with the bound \([-G_{\max}, G_{\max}]\) , we pass these to a off-the-shelf solver like scipy.optimize.linprog
Quadratic Problem#
A major limitation of the linear formulation is the non-uniqueness of the solution. Since the LP objective is often just a sum, many different “jagged” paths satisfy the constraints equally. In MRI, we prefer smooth waveforms. To achieve this, we shift to Quadratic Programming (QP), adding a cost on the variation of the gradient. In our implementation, we specifically minimize the second derivative (curvature), which is equivalent to minimizing the norm: \(\|G_{i+1} - 2 G_{i} + G_{i-1}\|^2_2\)
To minimize the variation we fix the endpoints \(g_{\text{start}}\) and \(g_{\text{end}}\)and define a reduced vector of free variables \(u = G_{1\dots N-2}\) . We build a sparse second-difference matrix \(R\) such that \(Rg\) yields the curvature, the quadratic cost then becomes: \[ f(u) = \frac12 u^\top H u + q^\top u + c \] Here \(H=2R^\top R\).
For the solver we will rely on osqp and therefore we adopt the following matrix formulation: \[ \min_{\mathbf{u}} \frac{1}{2} \mathbf{u}^T H \mathbf{u} + q^T \mathbf{u} \quad \text{s.t.} \quad A_{ub} \mathbf{u} \leq b_{ub} \]
Finding the minimum number of point over all axes.#
A strong limitation of the previous formulation is that we have to fix the number of points N beforehand, but in practice we want to find the minimum number of points needed to reach our target in k-space. We want the fastest “sprint” possible.
Luckily, if the solution works for N, it will also work for any other bigger N, and so we can find the “critical” N by doing a binary search over a range of estimation. We can use \(N_{\min} = \lfloor \Delta k / G_{\max} \rfloor\) and \(N_{\max} = N_{\min} + 2 \cdot \lfloor G_{\max} / (S_{\max} \Delta T)\rfloor\), That is the minimal time (considering no ramping gradients), and maximal time (considering maximal ramps up and down).
Finding the minimum N is an expensive endeavor, as it is not possible to check the feasibility of a linear programming problem without actually solving it (or at least attempting to). Each step of our binary search requires invoking the solver, which can become a bottleneck when generating complex trajectories with thousands of segments.
So far the problem is 1D, but in practice we have 3D gradients, and we want to find the optimal trajectory in 3D k-space. While the physics of each axis \(X,Y,Z\) are independent, the timing is not. Some axes may be “harder” to solve than others because they require a larger k-space displacement. However, all axes must arrive at the destination at the same time to keep the sequence coherent. We solve this by finding the critical Nx,Ny,Nz for each axis independently and then taking \(N = \max(N_x, N_y, N_z)\) as the final duration for all three.
Scaling things up#
Typical trajectories consist of a large number of shots (like in a multi-shot spiral or radial sequence), and each one needs a unique prephasor and a spoiler. This means we need to solve this optimization problem hundreds or thousands of times just to prep a single scan.
The most expensive part is finding the number N of points. However, as shots in a trajectory are often similar (e.g., rotating spokes in a radial scan), we can keep track of the critical N for a given set of constraints.
We achieve this by quantizing the input parameters. By defining a “quantum” (typically half the maximum gradient step), we can map the continuous k-space targets and boundary gradients into an integer space: \[\Delta k_q=\text{trunc}(\Delta k/\text{q})\] With \(q = 0.5 S_{\max} \Delta T\), that is half of a maximal step in gradient waveform. This quantization allows us to use memoization effectively, we store the computed \(N\) for a specific triplet of \(\Delta k_q, g_{\text{start}_q}, g_{\text{end}_q}\) in a cache if we encouter a similar constraint later, we skip the expensive binary search and solve the optimization directly with the cached N.
Conclusion#
Efficiently navigating k-space is a balancing act between physical limits and mathematical optimization. By framing the problem as a Quadratic Program, we can generate smooth, hardware-friendly gradients automatically. If you want to use this in your own work, check out the mri-nufft API; the connect_gradient function handles all this heavy lifting under the hood. This is also what drive the export to pulseq sequence
While this is nothing revolutionary per se, proposing similar convex optimization for MRI gradients date back to the late 90s and early 2000s. Having a robust, open-source implementation makes these techniques accessible to the broader community.
Beyond just connecting points, these same principles can be used for gradient projection. If you have an existing trajectory that is “illegal” (exceeds slew rate), you can project it onto the feasible set defined by your constraints, “fixing” the trajectory while staying as close as possible to your original design. But that is a story for another time.