From a271ee0e53d038b2d09041fb424e297e1cf7a7aa Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Mon, 25 Sep 2023 22:37:00 +0200 Subject: [PATCH] Add Seiffert spiral and shell trajectories --- examples/example_3D_trajectories.py | 159 +++++++++++++++-- src/mrinufft/__init__.py | 4 + src/mrinufft/trajectories/__init__.py | 4 + src/mrinufft/trajectories/trajectory3D.py | 208 +++++++++++++++++++++- 4 files changed, 354 insertions(+), 21 deletions(-) diff --git a/examples/example_3D_trajectories.py b/examples/example_3D_trajectories.py index c70b5d2a..c7a13957 100644 --- a/examples/example_3D_trajectories.py +++ b/examples/example_3D_trajectories.py @@ -66,9 +66,9 @@ def show_argument(function, arguments, one_shot, subfigure_size): # These options are used in the examples below as default values for all trajectories. # Trajectory parameters -Nc = 72 # Number of shots -Ns = 512 # Number of samples per shot -in_out = True # Choose between in-out or center-out trajectories +Nc = 100 # Number of shots +Ns = 500 # Number of samples per shot +in_out = False # Choose between in-out or center-out trajectories tilt = "uniform" # Choose the angular distance between shots nb_shells = 8 # Number of concentric shells for shell-type trajectories @@ -86,7 +86,7 @@ def show_argument(function, arguments, one_shot, subfigure_size): # and relying on different principles. # # 3D Cones -# ------ +# -------- # # A common pattern composed of 3D cones oriented all over within a sphere. # @@ -94,10 +94,12 @@ def show_argument(function, arguments, one_shot, subfigure_size): # # - ``Nc (int)``: number of individual shots # - ``Ns (int)``: number of samples per shot -# - ``tilt (str, float)``: angle between each consecutive shot (in radians). ``(default "uniform")`` -# - ``in_out (bool)``: define whether the shots should travel toward the center then outside -# (in-out) or not (center-out). ``(default False)`` -# - ``nb_zigzags (float)``: number of revolutions over a center-out shot. ``(default 5)`` +# - ``tilt (str, float)``: angle between each consecutive shot (in radians). +# ``(default "uniform")`` +# - ``in_out (bool)``: define whether the shots should travel toward +# the center then outside (in-out) or not (center-out). ``(default False)`` +# - ``nb_zigzags (float)``: number of revolutions over a center-out shot. +# ``(default 5)`` # - ``width (float)``: cone width factor, normalized to cover the k-space by default. # ``(default 1)`` # @@ -137,7 +139,7 @@ def show_argument(function, arguments, one_shot, subfigure_size): # is lengthened or the sampling rate is increased. # -arguments = [20, 50, 80, 200] +arguments = [10, 25, 40, 100] function = lambda x: mn.initialize_3D_cones(Nc, x, in_out=in_out) show_argument(function, arguments, one_shot=one_shot, subfigure_size=subfigure_size) @@ -169,8 +171,8 @@ def show_argument(function, arguments, one_shot, subfigure_size): # then going back to outer regions, often on the opposite side (radial, cones) # - center-out or center-center: when ``in_out=False`` the trajectory will start # at the center, but depending on the specific trajectory formula the path might -# end up in the outer regions (radial, spiral, cones, etc) or back to the center (rosette, -# lissajous). +# end up in the outer regions (radial, spiral, cones, etc) +# or back to the center (rosette, lissajous). # # Note that the behavior of both ``tilt`` and ``width`` are automatically adapted # to the changes to avoid having to update them too when switching ``in_out``. @@ -209,6 +211,82 @@ def show_argument(function, arguments, one_shot, subfigure_size): show_argument(function, arguments, one_shot=one_shot, subfigure_size=subfigure_size) +# %% +# Seiffert spirals / Yarnball +# --------------------------- +# +# A recent pattern with tightly controlled gradient norms using radially +# modulated Seiffert spirals, based on Jacobi elliptic functions. +# Note that Seiffert spirals more commonly refer to a curve evolving +# over a sphere surface rather than a volume, with the advantage of +# having a constant speed and angular velocity. The MR trajectory +# is obtained by increasing progressively the radius of the sphere. +# +# This implementation follows the proposition from [SMR18]_ based on +# works from [Er00]_ and [Br09]_. The pattern is also referred to as +# Yarnball by a different team [SB21]_, as a nod to the Yarn trajectory +# pictured in [IN95]_, even though both admittedly share little in common. +# +# Arguments: +# +# - ``Nc (int)``: number of individual shots. See 3D cones +# - ``Ns (int)``: number of samples per shot. See 3D cones +# - ``curve_index (float)``: Index controlling curvature from 0 (flat) to 1 (curvy). +# ``(default 0.3)`` +# - ``nb_revolutions (float)``: number of revolutions or elliptic periods. +# ``(default 1)`` +# - ``tilt (str, float)``: angle between each consecutive shot (in radians). +# ``(default "uniform")``. See 3D cones +# - ``in_out (bool)``: define whether the shots should travel toward the center +# then outside (in-out) or not (center-out). ``(default False)``. See 3D cones +# + +trajectory = mn.initialize_3D_seiffert_spiral(Nc, Ns, in_out=in_out) +display_3D_trajectory( + trajectory, + nb_repetitions=1, + Nc=Nc, + Ns=Ns, + size=figure_size, + one_shot=one_shot, + per_plane=False, +) +plt.show() + + +# %% +# ``curve_index (float)`` +# ~~~~~~~~~~~~~~~~~~~~~~~ +# +# An index defined over :math:`[0, 1)` controling the curvature, with :math:`0` +# corresponding to a planar spiral, and increasing the length and exploration of +# the curve while asymptotically approaching :math:`1`. +# + +arguments = [0, 0.3, 0.9, 0.99] +function = lambda x: mn.initialize_3D_seiffert_spiral( + Nc, Ns, in_out=in_out, curve_index=x +) +show_argument(function, arguments, one_shot=one_shot, subfigure_size=subfigure_size) + + +# %% +# ``nb_revolutions (float)`` +# ~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Number of revolutions, or simply the number of times a curve reaches its +# original orientation. For regular Seiffert spirals, it corresponds to the +# number of times the shot reaches the starting pole of the sphere. It +# subsequently defines the length of the curve. +# + +arguments = [0.5, 1, 1.5, 2] +function = lambda x: mn.initialize_3D_seiffert_spiral( + Nc, Ns, in_out=in_out, nb_revolutions=x +) +show_argument(function, arguments, one_shot=one_shot, subfigure_size=subfigure_size) + + # %% # Shell trajectories # ================== @@ -224,7 +302,7 @@ def show_argument(function, arguments, one_shot, subfigure_size): # other trajectories sharing this principle. # # This implementation follows the proposition from [YRB06]_ but the idea -# is much older and can be traced back at least to [PN95]_. +# is much older and can be traced back at least to [IN95]_. # # Arguments: # @@ -367,7 +445,7 @@ def show_argument(function, arguments, one_shot, subfigure_size): # shells. # # An angle of :math:`\pi / 2` allows reaching the best shot length homogeneity, -# and partition the spheres into several connex curves composed of exactly +# and it partitions the spheres into several connex curves composed of exactly # two shots. # @@ -378,16 +456,67 @@ def show_argument(function, arguments, one_shot, subfigure_size): show_argument(function, arguments, one_shot=one_shot, subfigure_size=subfigure_size) +# %% +# Seiffert shells +# --------------- +# +# An exclusive trajectory composed of re-arranged Seiffert spirals +# covering concentric shells. All curves have a constant speed and +# angular velocity, depending on the size of the sphere they belong to. +# +# This implementation is inspired by the propositions from [YRB06]_ and [SMR18]_, +# and also based on works from [Er00]_ and [Br09]_. +# +# Arguments: +# +# - ``Nc (int)``: number of individual shots. See 3D cones +# - ``Ns (int)``: number of samples per shot. See 3D cones +# - ``curve_index (float)``: Index controlling curvature from 0 (flat) to 1 (curvy). +# ``(default 0.3)``. See Seiffert spirals +# - ``nb_revolutions (float)``: number of revolutions or elliptic periods. +# ``(default 1)``. See Seiffert spirals +# - ``shell_tilt (str, float)``: angle between each consecutive shell (in radians). +# ``(default "intergaps")``. See helical shells +# - ``shot_tilt (str, float)``: angle between each consecutive shot +# over a sphere (in radians). ``(default "uniform")``. See helical shells +# + +trajectory = mn.initialize_3D_seiffert_shells(Nc, Ns, nb_shells) +display_3D_trajectory( + trajectory, + nb_repetitions=1, + Nc=Nc, + Ns=Ns, + size=figure_size, + one_shot=one_shot, + per_plane=False, +) +plt.show() + + # %% # References # ========== # -# .. [PN95] Irarrazabal, Pablo, and Dwight G. Nishimura. +# .. [IN95] Irarrazabal, Pablo, and Dwight G. Nishimura. # "Fast three dimensional magnetic resonance imaging." -# Magnetic Resonance in Medicine 33, no. 5 (1995): 656-662 +# Magnetic Resonance in Medicine 33, no. 5 (1995): 656-662. +# .. [Er00] Erdös, Paul. +# "Spiraling the earth with C. G. J. Jacobi." +# American Journal of Physics 68, no. 10 (2000): 888-895. # .. [YRB06] Shu, Yunhong, Stephen J. Riederer, and Matt A. Bernstein. # "Three‐dimensional MRI with an undersampled spherical shells trajectory." # Magnetic Resonance in Medicine 56, no. 3 (2006): 553-562. +# .. [Br09] Brizard, Alain J. +# "A primer on elliptic functions with applications in classical mechanics." +# European journal of physics 30, no. 4 (2009): 729. # .. [HM11] Gerlach, Henryk, and Heiko von der Mosel. # "On sphere-filling ropes." # The American Mathematical Monthly 118, no. 10 (2011): 863-876 +# .. [SMR18] Speidel, Tobias, Patrick Metze, and Volker Rasche. +# "Efficient 3D Low-Discrepancy k-Space Sampling +# Using Highly Adaptable Seiffert Spirals." +# IEEE Transactions on Medical Imaging 38, no. 8 (2018): 1833-1840. +# .. [SB21] Stobbe, Robert W., and Christian Beaulieu. +# "Three‐dimensional Yarnball k‐space acquisition for accelerated MRI." +# Magnetic Resonance in Medicine 85, no. 4 (2021): 1840-1854. diff --git a/src/mrinufft/__init__.py b/src/mrinufft/__init__.py index c8b25cc7..b85a65b8 100644 --- a/src/mrinufft/__init__.py +++ b/src/mrinufft/__init__.py @@ -24,8 +24,10 @@ initialize_2D_waves, initialize_3D_from_2D_expansion, initialize_3D_cones, + initialize_3D_seiffert_spiral, initialize_3D_helical_shells, initialize_3D_annular_shells, + initialize_3D_seiffert_shells, display_2D_trajectory, display_3D_trajectory, ) @@ -45,8 +47,10 @@ "initialize_2D_waves", "initialize_3D_from_2D_expansion", "initialize_3D_cones", + "initialize_3D_seiffert_spiral", "initialize_3D_helical_shells", "initialize_3D_annular_shells", + "initialize_3D_seiffert_shells", "display_2D_trajectory", "display_3D_trajectory", ] diff --git a/src/mrinufft/trajectories/__init__.py b/src/mrinufft/trajectories/__init__.py index 9157fbcf..f5c72282 100644 --- a/src/mrinufft/trajectories/__init__.py +++ b/src/mrinufft/trajectories/__init__.py @@ -14,8 +14,10 @@ from .trajectory3D import ( initialize_3D_from_2D_expansion, initialize_3D_cones, + initialize_3D_seiffert_spiral, initialize_3D_helical_shells, initialize_3D_annular_shells, + initialize_3D_seiffert_shells, ) from .display import ( @@ -35,8 +37,10 @@ "initialize_2D_waves", "initialize_3D_from_2D_expansion", "initialize_3D_cones", + "initialize_3D_seiffert_spiral", "initialize_3D_helical_shells", "initialize_3D_annular_shells", + "initialize_3D_seiffert_shells", "display_2D_trajectory", "display_3D_trajectory", ] diff --git a/src/mrinufft/trajectories/trajectory3D.py b/src/mrinufft/trajectories/trajectory3D.py index d0b988f0..b22dc34d 100644 --- a/src/mrinufft/trajectories/trajectory3D.py +++ b/src/mrinufft/trajectories/trajectory3D.py @@ -2,6 +2,8 @@ import numpy as np import numpy.linalg as nl +from scipy.special import ellipj, ellipk + from .expansions import ( stack_2D_to_3D_expansion, rotate_2D_to_3D_expansion, @@ -17,9 +19,9 @@ from .utils import Ry, Rz, Rv, initialize_tilt, KMAX -################################ -# 3D TRAJECTORY INITIALIZATION # -################################ +############################ +# FREEFORM 3D TRAJECTORIES # +############################ def initialize_3D_from_2D_expansion( @@ -95,7 +97,8 @@ def initialize_3D_cones(Nc, Ns, tilt="golden", in_out=False, nb_zigzags=5, width tilt : str, float, optional Tilt of the cones, by default "golden" in_out : bool, optional - Whether the cones are going in and out, by default False + Whether the curves are going in-and-out or start from the center, + by default False nb_zigzags : float, optional Number of zigzags of the cones, by default 5 width : float, optional @@ -137,13 +140,108 @@ def initialize_3D_cones(Nc, Ns, tilt="golden", in_out=False, nb_zigzags=5, width return trajectory.reshape((Nc, Ns, 3)) +def initialize_3D_seiffert_spiral( + Nc, Ns, curve_index=0.2, nb_revolutions=1, tilt="golden", in_out=False +): + """Initialize 3D trajectories with modulated Seiffert spirals. + + Initially introduced in [SMR18]_, but also proposed later as "Yarnball" + in [SB21]_ as a nod to [IN95]_. The implementation is based on work + from [Er00]_ and [Br09]_, using Jacobi elliptic functions rather than + auxiliary theta functions. + + Parameters + ---------- + Nc : int + Number of shots + Ns : int + Number of samples per shot + curve_index : float + Index controlling curve from 0 (flat) to 1 (curvy), by default 0.3 + nb_revolutions : float + Number of revolutions, i.e. times the polar angle of the curves + passes through 0, by default 1 + tilt : str, float, optional + Angle between shots around z-axis over precession, by default "golden" + in_out : bool + Whether the curves are going in-and-out or start from the center, + by default False + + Returns + ------- + array_like + 3D Seiffert spiral trajectory + + References + ---------- + .. [IN95] Irarrazabal, Pablo, and Dwight G. Nishimura. + "Fast three dimensional magnetic resonance imaging." + Magnetic Resonance in Medicine 33, no. 5 (1995): 656-662. + .. [Er00] Erdös, Paul. + "Spiraling the earth with C. G. J. jacobi." + American Journal of Physics 68, no. 10 (2000): 888-895. + .. [Br09] Brizard, Alain J. + "A primer on elliptic functions with applications in classical mechanics." + European journal of physics 30, no. 4 (2009): 729. + .. [SMR18] Speidel, Tobias, Patrick Metze, and Volker Rasche. + "Efficient 3D Low-Discrepancy k-Space Sampling + Using Highly Adaptable Seiffert Spirals." + IEEE Transactions on Medical Imaging 38, no. 8 (2018): 1833-1840. + .. [SB21] Stobbe, Robert W., and Christian Beaulieu. + "Three‐dimensional Yarnball k‐space acquisition for accelerated MRI." + Magnetic Resonance in Medicine 85, no. 4 (2021): 1840-1854. + + """ + # Normalize ellipses integrations by the requested period + trajectory = np.zeros((Nc, Ns // (1 + in_out), 3)) + period = 4 * ellipk(curve_index**2) + times = np.linspace(0, nb_revolutions * period, Ns // (1 + in_out), endpoint=False) + + # Initialize first shot + jacobi = ellipj(times, curve_index**2) + trajectory[0, :, 0] = jacobi[0] * np.cos(curve_index * times) + trajectory[0, :, 1] = jacobi[0] * np.sin(curve_index * times) + trajectory[0, :, 2] = jacobi[1] + + # Make it volumetric instead of just a sphere surface + magnitudes = np.sqrt(np.linspace(0, 1, Ns // (1 + in_out))) + trajectory = magnitudes.reshape((1, -1, 1)) * trajectory + + # Determine mostly evenly distributed points on sphere + points = np.zeros((Nc, 3)) + phi = initialize_tilt(tilt) * np.arange(Nc) + points[:, 0] = np.linspace(-1, 1, Nc) + radius = np.sqrt(1 - points[:, 0] ** 2) + points[:, 1] = np.cos(phi) * radius + points[:, 2] = np.sin(phi) * radius + + # Rotate initial shot Nc times + for i in np.arange(1, Nc)[::-1]: + v1 = np.array((1, 0, 0)) + v2 = points[i] + rotation = Rv(v1, v2, normalize=False) + trajectory[i] = (rotation @ trajectory[0].T).T + + # Handle in_out case + if in_out: + first_half_traj = np.copy(trajectory) + first_half_traj = -first_half_traj[:, ::-1] + trajectory = np.concatenate([first_half_traj, trajectory], axis=1) + return KMAX * trajectory + + +############################### +# SHELL-BASED 3D TRAJECTORIES # +############################### + + def initialize_3D_helical_shells( Nc, Ns, nb_shells, spiral_reduction=1, shell_tilt="intergaps", shot_tilt="uniform" ): """Initialize 3D trajectories with helical shells. The implementation follows the proposition from [YRB06]_ - but the idea can be traced back much further [PN95]_. + but the idea can be traced back much further [IN95]_. Parameters ---------- @@ -170,7 +268,7 @@ def initialize_3D_helical_shells( .. [YRB06] Shu, Yunhong, Stephen J. Riederer, and Matt A. Bernstein. "Three‐dimensional MRI with an undersampled spherical shells trajectory." Magnetic Resonance in Medicine 56, no. 3 (2006): 553-562. - .. [PN95] Irarrazabal, Pablo, and Dwight G. Nishimura. + .. [IN95] Irarrazabal, Pablo, and Dwight G. Nishimura. "Fast three dimensional magnetic resonance imaging." Magnetic Resonance in Medicine 33, no. 5 (1995): 656-662 @@ -333,6 +431,104 @@ def initialize_3D_annular_shells( return KMAX * trajectory +def initialize_3D_seiffert_shells( + Nc, + Ns, + nb_shells, + curve_index=0.5, + nb_revolutions=1, + shell_tilt="uniform", + shot_tilt="uniform", +): + """Initialize 3D trajectories with Seiffert shells. + + The implementation is based on work from [Er00]_ and [Br09]_, + using Jacobi elliptic functions to define Seiffert spirals + over shell/sphere surfaces. + + Parameters + ---------- + Nc : int + Number of shots + Ns : int + Number of samples per shot + nb_shells : int + Number of concentric shells/spheres + curve_index : float + Index controlling curve from 0 (flat) to 1 (curvy), by default 0.5 + nb_revolutions : float + Number of revolutions, i.e. times the curve passes through the upper-half + of the z-axis, by default 1 + shell_tilt : str, float, optional + Angle between consecutive shells along z-axis, by default "intergaps" + shot_tilt : str, float, optional + Angle between shots over a shell surface along z-axis, by default "uniform" + + Returns + ------- + array_like + 3D Seiffert shell trajectory + + References + ---------- + .. [IN95] Irarrazabal, Pablo, and Dwight G. Nishimura. + "Fast three dimensional magnetic resonance imaging." + Magnetic Resonance in Medicine 33, no. 5 (1995): 656-662. + .. [Er00] Erdös, Paul. + "Spiraling the earth with C. G. J. jacobi." + American Journal of Physics 68, no. 10 (2000): 888-895. + .. [Br09] Brizard, Alain J. + "A primer on elliptic functions with applications in classical mechanics." + European journal of physics 30, no. 4 (2009): 729. + + """ + # Check arguments validity + if Nc < nb_shells: + raise ValueError("Argument nb_shells should not be higher than Nc.") + trajectory = np.zeros((Nc, Ns, 3)) + + # Attribute shots to shells following a prescribed density + Nc_per_shell = np.ones(nb_shells).astype(int) + density = np.arange(1, nb_shells + 1) ** 2 # simplified version + for _ in range(Nc - nb_shells): + idx = np.argmax(density / Nc_per_shell) + Nc_per_shell[idx] += 1 + + # Normalize ellipses integrations by the requested period + period = 4 * ellipk(curve_index**2) + times = np.linspace(0, nb_revolutions * period, Ns, endpoint=False) + + # Create shells one by one + count = 0 + radii = (0.5 + np.arange(nb_shells)) / nb_shells + for i in range(nb_shells): + # Prepare shell parameters + Ms = Nc_per_shell[i] + k0 = radii[i] + + # Initialize first shot + jacobi = ellipj(times, curve_index**2) + trajectory[count, :, 0] = k0 * jacobi[0] * np.cos(curve_index * times) + trajectory[count, :, 1] = k0 * jacobi[0] * np.sin(curve_index * times) + trajectory[count, :, 2] = k0 * jacobi[1] + + # Rotate first shot Ms times to create the shell + rotation = Rz(initialize_tilt(shot_tilt, Ms)) + for j in range(1, Ms): + trajectory[count + j] = trajectory[count + j - 1] @ rotation + + # Apply shell tilt + rotation = Rz(i * initialize_tilt(shell_tilt, nb_shells)) + trajectory[count : count + Ms] = trajectory[count : count + Ms] @ rotation + count += Ms + return KMAX * trajectory + + +######### +# UTILS # +######### + + def duplicate_per_axes(trajectory, axes=(0, 1, 2)): """ Duplicate a trajectory along the specified axes.