Vds spiral

examples/examples_vds_spiral.py

@@ -0,0 +1,71 @@
+"""Variable density Spiral, based on the MATLAB implementation of Brian Hargreaves."""
+
+# ---
+# jupyter:
+#   jupytext:
+#     formats: py:percent
+#     text_representation:
+#       extension: .py
+#       format_name: percent
+#       format_version: '1.3'
+#       jupytext_version: 1.16.0
+#   kernelspec:
+#     display_name: sim
+#     language: python
+#     name: sim
+# ---
+
+# %%
+from mrinufft.trajectories.trajectory2D import initialize_2D_vds_spiral
+from mrinufft.trajectories.utils import Gammas
+
+# %%
+samples = initialize_2D_vds_spiral(
+    Nc=1,
+    Fcoeff=[100, -24],
+    raster_time=0.000004,
+    res=1,
+    oversamp=4,
+    gmax=4,
+    smax=15000,
+    gamma=4258,
+    in_out=True,
+)
+
+# %%
+samples.shape
+
+# %%
+from mrinufft.trajectories.display import display_2D_trajectory
+
+# %%
+display_2D_trajectory(samples, linewidth=0.01)
+
+# %%
+from mrinufft.trajectories.display import display_gradients
+
+# %%
+display_gradients(samples, show_constraints=True)
+
+# %%
+max(samples.flatten())
+
+# %%
+len(samples[0])
+import matplotlib.pyplot as plt
+
+# %%
+plt.plot(samples[0, :, 0], samples[0, :, 1])
+
+# %%
+len(samples[0])
+
+# %%
+len(samples[0]) * 4e-6 * 3
+
+# %%
+import numpy as np
+
+np.sqrt(-1)
+
+# %%

src/mrinufft/trajectories/maths.py

@@ -2,6 +2,8 @@
 
 import numpy as np
 import numpy.linalg as nl
+from .utils import DEFAULT_GMAX, DEFAULT_SMAX, Gammas
+
 
 CIRCLE_PACKING_DENSITY = np.pi / (2 * np.sqrt(3))
 EIGENVECTOR_2D_FIBONACCI = (0.4656, 0.6823, 1)
@@ -341,3 +343,71 @@ def generate_fibonacci_sphere(nb_points, epsilon=0.25):
     fibonacci_sphere[:, 1] = np.sin(theta) * np.sin(phi)
     fibonacci_sphere[:, 2] = np.cos(phi)
     return fibonacci_sphere
+
+
+####################
+# Variable Density #
+####################
+
+
+def findq2r2(
+    r,
+    r1,
+    T,
+    Ts,
+    N,
+    Fcoeff,
+    rmax,
+    smax=DEFAULT_SMAX,
+    gmax=DEFAULT_GMAX,
+    gamma=Gammas.H,
+):
+    """
+    Calculate the second derivative of r and q (theta) for vds spirals.
+
+    slew-limited or FOV-limited r(t) and q(t) waveforms such that
+    :math:`k(t) = r(t)exp(i*q(t))` Where the FOV is a function of k-space radius (r)
+    FOV = Fcoeff(1) + Fcoeff(2)*r/rmax + Fcoeff(3)*(r/rmax)^2 + ... ; F(1) in cm. F(2)
+    in cm^2. F(3) in cm^3, etc.
+    """
+    F = 0  # FOV function value for this r.
+    dFdr = 0  # dFOV/dr for this value of r.
+    for rind in range(len(Fcoeff)):
+        F += Fcoeff[rind] * (r / rmax) ** rind
+        if rind > 0:
+            dFdr += rind * Fcoeff[rind] * (r / rmax) ** (rind - 1) / rmax
+
+    GmaxFOV = 1 / (gamma * F * Ts)  # FOV limit on G
+    Gmax = min(GmaxFOV, gmax)
+
+    maxr1 = np.sqrt((gamma * Gmax) ** 2 / (1 + (2 * np.pi * F * r / N) ** 2))
+
+    if r1 > maxr1:
+        r2 = (maxr1 - r1) / T
+    else:
+        twopiFoN = 2 * np.pi * F / N
+        twopiFoN2 = twopiFoN**2
+
+        # A,B,C are coefficients of the equation which equates
+        # the slew rate calculated from r,r1,r2 with the
+        # maximum gradient slew rate.
+        #
+        # A*r2*r2 + B*r2 + C  =  0
+        #
+        # A,B,C are in terms of F,dF/dr,r,r1, N and smax.
+
+        A = 1 + twopiFoN2 * r**2
+        B = 2 * twopiFoN2 * r * r1**2 + 2 * twopiFoN2 / F * dFdr * r**2 * r1**2
+        C = (
+            twopiFoN2**2 * r**2 * r1**4
+            + 4 * twopiFoN2 * r1**4
+            + (2 * np.pi / N * dFdr) ** 2 * r**2 * r1**4
+            + 4 * twopiFoN2 / F * dFdr * r * r1**4
+            - (gamma) ** 2 * smax**2
+        )
+        # Use the big root
+        r2 = (-B + np.sqrt(B**2 - 4 * A * C)) / (2 * A)
+
+    q2 = 2 * np.pi / N * dFdr * r1**2 + 2 * np.pi * F / N * r2
+
+    return q2, r2

src/mrinufft/trajectories/trajectory2D.py

@@ -5,9 +5,17 @@
 from scipy.interpolate import CubicSpline
 
 from .gradients import patch_center_anomaly
-from .maths import R2D, compute_coprime_factors, is_from_fibonacci_sequence
+from .maths import R2D, compute_coprime_factors, is_from_fibonacci_sequence, findq2r2
 from .tools import rotate
-from .utils import KMAX, initialize_algebraic_spiral, initialize_tilt
+from .utils import (
+    DEFAULT_GMAX,
+    DEFAULT_SMAX,
+    DEFAULT_RASTER_TIME,
+    KMAX,
+    Gammas,
+    initialize_algebraic_spiral,
+    initialize_tilt,
+)
 
 #####################
 # CIRCULAR PATTERNS #
@@ -208,6 +216,97 @@ def initialize_2D_fibonacci_spiral(Nc, Ns, spiral_reduction=1, patch_center=True
     return trajectory
 
 
+def initialize_2D_vds_spiral(
+    Nc,
+    Fcoeff,
+    res,
+    oversamp=4,
+    raster_time=DEFAULT_RASTER_TIME,
+    gmax=DEFAULT_GMAX,
+    smax=DEFAULT_SMAX,
+    in_out=False,
+    gamma=Gammas.H,
+):
+    """Initialize a 2D variable-density spiral trajectory.
+
+    Based on the MATLAB Implementation of Brian Hargreaves.
+
+    Parameters
+    ----------
+    Nc: int
+        Number of Shots (interleave of the spiral)
+    Ns:
+    """
+    # Compared to the MATLAB Implementation, a few variable have been renamed
+    # MATLAB -> Python
+    # N -> Ncc
+    # T -> raster_time
+
+    To = raster_time / oversamp  # To is the period with oversampling.
+    rmax = 0.5 / res
+
+    q0 = q1 = r0 = r1 = 0.0  # initialize theta, theta', r, r'
+
+    theta = np.zeros(1_000_000)
+    r = np.zeros(1_000_000)
+    time = np.zeros(1_000_000)
+    t = 0.0
+    count = 0
+    Ncc = Nc
+    if in_out:
+        Ncc = Ncc * 2  # create twice center_out trajectories.
+
+    while r0 < rmax and count < 1_000_000 - 1:
+        q2, r2 = findq2r2(
+            r0,
+            r1,
+            To,
+            raster_time,
+            Ncc,
+            Fcoeff,
+            rmax=rmax,
+            smax=smax,
+            gmax=gmax,
+            gamma=gamma,
+        )
+
+        # Integrate for r, r', theta and theta'
+        q1 += q2 * To
+        q0 += q1 * To
+        r1 += r2 * To
+        r0 += r1 * To
+
+        t += To
+
+        # Store
+        count += 1
+        theta[count] = q0
+        r[count] = r0
+        time[count] = t
+
+    # Generate the spiral from r and theta
+    r = r[:count]
+    theta = theta[:count]
+    # Convert to cartesian coordinates
+    #
+    trajectory = r * np.exp(1j * theta)
+    trajectory = np.stack([trajectory.real, trajectory.imag], axis=-1)
+    if in_out:
+        trajectory2 = r * np.exp(1j * (theta + np.pi))
+        trajectory2 = np.stack([trajectory2.real, trajectory2.imag], axis=-1)
+        trajectory = np.concatenate(
+            [trajectory2[::-1], trajectory], axis=0
+        )  # flip and put the center in the middle
+    trajectory = np.repeat(trajectory[None, ...], Nc, axis=0)
+
+    # Rotate the first shot Nc times
+    rotation = R2D(initialize_tilt("uniform", Nc)).T
+    for i in range(1, Nc):
+        trajectory[i] = trajectory[i - 1] @ rotation
+    trajectory = KMAX * trajectory / np.max(nl.norm(trajectory, axis=-1))
+    return trajectory
+
+
 def initialize_2D_cones(Nc, Ns, tilt="uniform", in_out=False, nb_zigzags=5, width=1):
     """Initialize a 2D cone trajectory.