diff --git a/examples/conftest.py b/examples/conftest.py index 395ee8cb..9eb34cc7 100644 --- a/examples/conftest.py +++ b/examples/conftest.py @@ -12,11 +12,11 @@ """ -from pathlib import Path import runpy -import pytest -import matplotlib as mpl +from pathlib import Path +import matplotlib as mpl +import pytest mpl.use("agg") diff --git a/examples/example_2D_trajectories.py b/examples/example_2D_trajectories.py index efe2c6c3..b4498f13 100644 --- a/examples/example_2D_trajectories.py +++ b/examples/example_2D_trajectories.py @@ -23,7 +23,7 @@ # Internal import mrinufft as mn - +import mrinufft.trajectories.maths as mntm from mrinufft import display_2D_trajectory @@ -168,10 +168,10 @@ def show_trajectory(trajectory, one_shot, figure_size): # Spiral # ------ # -# A generalized function that generates spirals defined through the -# :math:`r = a \theta^{1/n}` equality, with :math:`r` the radius and -# :math:`\theta` the polar angle. Note that the most common spirals, -# Archimedes and Fermat, are subcases of this equation. +# A generalized function that generates algebraic spirals defined +# through the :math:`r = a \theta^n` equation, with :math:`r` the radius, +# :math:`\theta` the polar angle and :math:`n` the spiral power. +# Common algebraic spirals include Archimedes, Fermat and Galilean spirals. # # Arguments: # @@ -210,18 +210,125 @@ def show_trajectory(trajectory, one_shot, figure_size): # ``spiral (str, float)`` # ~~~~~~~~~~~~~~~~~~~~~~~ # +# The algebraic spiral power defined through :math:`n` in the +# :math:`r = a \theta^n` equality, with :math:`r` the radius and +# :math:`\theta` the polar angle. It defines the gradient behavior, +# and therefore the distance between consecutive points and the shape +# of the spiral. It does not affect the number of revolutions, but +# rather the curve length and point distribution. Spirals with small +# :math:`n` (close to 0) tend to have radial behaviors +# around the center, and dedicate more points towards curved edges. # -# The shape of the spiral defined through :math:`n` in the -# :math:`r = a \theta^{1/n}` equality, with :math:`r` the radius and -# :math:`\theta` the polar angle. Both ``"archimedes"`` and ``"fermat"`` -# are available as string options for convenience. +# ``"archimedes"`` (1), ``"fermat"`` (0.5) and ``"galilean"`` (2) are available +# as string options for convenience. Algebraic spirals with negative powers, +# such as hyperbolic or lithuus spirals, are not considered relevant because +# of their asymptotic behavior around the center. # -arguments = ["archimedes", "fermat", 0.5, 1.5] +arguments = ["galilean", "archimedes", "fermat", 1 / 4] function = lambda x: mn.initialize_2D_spiral(Nc, Ns, tilt=tilt, spiral=x, in_out=in_out) show_argument(function, arguments, one_shot=one_shot, subfigure_size=subfigure_size) +# %% +# ``patch_center (float)`` +# ~~~~~~~~~~~~~~~~~~~~~~~~ +# +# A slew rate anomaly is present at the center of algebraic spirals +# when their power is inferior to 1 (e.g. Fermat's) and parameterized +# through their angles in the above equation. +# +# To fix this problem, points at the center are re-arranged along +# the spiral until the gradients are monotically increasing from +# the center to the edges. This correction can be deactivated, +# but it is generally preferred to keep it. +# +# The spiral path is not changed, but the density can be altered +# over the first few samples. However the difference is extremely +# subtle, as shown below. +# + +arguments = [False, True] +function = lambda x: mn.initialize_2D_spiral( + Nc, + Ns, + patch_center=x, +) +show_argument(function, arguments, one_shot=one_shot, subfigure_size=subfigure_size) + + +# %% +# Fibonacci spiral +# ---------------- +# +# A non-algebraic spiral trajectory based on the Fibonacci sequence, +# reproducing the proposition from [CA99]_ in order to generate +# a uniform distribution with center-out shots. +# +# The number of shots is required to belong to the Fibonacci +# sequence for the trajectory definition to be relevant. +# +# Arguments: +# +# - ``Nc (int)``: number of individual shots. See radial +# - ``Ns (int)``: number of samples per shot. See radial +# - ``spiral_reduction (float)``: factor used to reduce the automatic spiral length. ``(default 1)`` +# - ``patch_center (bool)``: whether the spiral anomaly at the center should be patched. +# ``(default True)`` +# + +Nc_fibonacci = mntm.get_closest_fibonacci_number(Nc) +trajectory = mn.initialize_2D_fibonacci_spiral(Nc_fibonacci, Ns) +show_trajectory(trajectory, figure_size=figure_size, one_shot=one_shot) + + +# %% +# ``spiral_reduction (float)`` +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Factor used to reduce the automatic spiral length. In opposition to +# ``initialize_2D_spiral``, the number of spiral revolutions here +# is automatically determined from ``Ns`` and ``Nc`` to match a uniform +# density over the k-space sphere. It can lead to unrealistically +# strong gradients, and therefore we provide this factor to reduce the +# spiral length, which makes k-space denser along the shorter shots. +# + +arguments = [0.5, 1, 2, 3] +function = lambda x: mn.initialize_2D_fibonacci_spiral( + Nc_fibonacci, + Ns, + spiral_reduction=x, +) +show_argument(function, arguments, one_shot=one_shot, subfigure_size=subfigure_size) + + +# %% +# ``patch_center (float)`` +# ~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Similarly to algebraic spirals from ``initialize_2D_spiral``, +# the trajectory definition creates small anomalies at the center +# that makes slew rate requirements needlessly high. +# +# It is here related to the uniform density that requires central +# samples to be more strongly spaced than anywhere else because +# most shots start close to the center. +# +# The spiral path can be altered over the first few samples, +# but generally the difference is extremely subtle, as shown +# below. +# + +arguments = [False, True] +function = lambda x: mn.initialize_2D_fibonacci_spiral( + Nc_fibonacci, + Ns, + patch_center=x, +) +show_argument(function, arguments, one_shot=one_shot, subfigure_size=subfigure_size) + + # %% # Cones # ----- diff --git a/examples/example_3D_trajectories.py b/examples/example_3D_trajectories.py index 7f33a434..047e47e7 100644 --- a/examples/example_3D_trajectories.py +++ b/examples/example_3D_trajectories.py @@ -29,7 +29,6 @@ # Internal import mrinufft as mn - from mrinufft import display_2D_trajectory, display_3D_trajectory diff --git a/examples/example_density.py b/examples/example_density.py index 98950aee..20f45b3c 100644 --- a/examples/example_density.py +++ b/examples/example_density.py @@ -12,14 +12,13 @@ """ import brainweb_dl as bwdl - import matplotlib.pyplot as plt import numpy as np -from mrinufft import get_density, get_operator, check_backend + +from mrinufft import check_backend, get_density, get_operator from mrinufft.trajectories import initialize_2D_radial from mrinufft.trajectories.display import display_2D_trajectory - # %% # Create sample data # ------------------ diff --git a/examples/example_display_config.py b/examples/example_display_config.py index 925d40b9..1bf03679 100644 --- a/examples/example_display_config.py +++ b/examples/example_display_config.py @@ -8,16 +8,14 @@ You can tune these parameters to your own taste and needs. """ +import matplotlib as mpl +import matplotlib.pyplot as plt + # %% import numpy as np -import matplotlib.pyplot as plt -import matplotlib as mpl -from mrinufft import displayConfig, display_2D_trajectory, display_3D_trajectory -from mrinufft.trajectories import ( - initialize_2D_spiral, - conify, -) +from mrinufft import display_2D_trajectory, display_3D_trajectory, displayConfig +from mrinufft.trajectories import conify, initialize_2D_spiral # Trajectory parameters Nc = 120 # Number of shots diff --git a/examples/example_gif_2D.py b/examples/example_gif_2D.py index 5ae7039d..31133ea6 100644 --- a/examples/example_gif_2D.py +++ b/examples/example_gif_2D.py @@ -7,39 +7,30 @@ """ -# %% +import joblib import matplotlib.pyplot as plt -import mrinufft.trajectories.display as mtd -import mrinufft.trajectories.trajectory2D as mtt import numpy as np - -import joblib from PIL import Image, ImageSequence - +import mrinufft.trajectories.display as mtd +import mrinufft.trajectories.trajectory2D as mtt from mrinufft.trajectories.display import displayConfig +# Options -# %% [markdown] -# # Options - -# %% Nc = 16 Ns = 200 - nb_repetitions = 8 -# %% one_shot = 0 figsize = 4 -nb_fps = 12 -name_slow = 1 +nb_frames = 3 +duration = 150 # seconds -# %% [markdown] -# # Generation -# %% +# Generation + # Initialize trajectory function functions = [ # Radial @@ -52,11 +43,12 @@ ("Radial", lambda x: mtt.initialize_2D_radial(Nc, Ns, tilt=x)), ("Radial", lambda x: mtt.initialize_2D_radial(Nc, Ns, tilt="uniform")), # Spiral + ("Spiral", lambda x: mtt.initialize_2D_spiral(Nc, Ns, nb_revolutions=x)), ("Spiral", lambda x: mtt.initialize_2D_spiral(Nc, Ns, spiral=x)), ("Spiral", lambda x: mtt.initialize_2D_spiral(Nc, Ns, spiral=x)), ("Spiral", lambda x: mtt.initialize_2D_spiral(Nc, Ns, nb_revolutions=x)), ("Spiral", lambda x: mtt.initialize_2D_spiral(Nc, Ns, nb_revolutions=x)), - ("Spiral", lambda x: mtt.initialize_2D_spiral(Nc, Ns, nb_revolutions=0)), + ("Spiral", lambda x: mtt.initialize_2D_spiral(Nc, Ns, nb_revolutions=1e-5)), # Cones ("Cones", lambda x: mtt.initialize_2D_cones(Nc, Ns, nb_zigzags=x)), ("Cones", lambda x: mtt.initialize_2D_cones(Nc, Ns, width=x)), @@ -77,7 +69,6 @@ ("Rings", lambda x: mtt.initialize_2D_rings(x, Ns, nb_rings=nb_repetitions)[::-1]), ("Rings", lambda x: mtt.initialize_2D_rings(Nc, Ns, nb_rings=nb_repetitions)[::-1]), # Rosette - ("Rosette", lambda x: mtt.initialize_2D_rosette(Nc, Ns, coprime_index=0)), ("Rosette", lambda x: mtt.initialize_2D_rosette(Nc, Ns, coprime_index=x)), ("Rosette", lambda x: mtt.initialize_2D_rosette(Nc, Ns, coprime_index=30)), # Waves @@ -85,55 +76,54 @@ ("Waves", lambda x: mtt.initialize_2D_waves(Nc, Ns, nb_zigzags=6 * x, width=x)), ("Waves", lambda x: mtt.initialize_2D_waves(Nc, Ns, nb_zigzags=6, width=1)), # Lissajous - ("Lissajous", lambda x: mtt.initialize_2D_lissajous(Nc, Ns, density=1)), ("Lissajous", lambda x: mtt.initialize_2D_lissajous(Nc, Ns, density=x)), ("Lissajous", lambda x: mtt.initialize_2D_lissajous(Nc, Ns, density=10)), ] -# %% # Initialize trajectory arguments arguments = [ # Radial - np.around(np.linspace(1, Nc, nb_fps)).astype(int), # Nc - np.linspace(2 * np.pi / Nc, np.pi / Nc, nb_fps // 2), # tilt - np.around(np.sin(np.linspace(0, 2 * np.pi, nb_fps // 2))).astype(bool), # in_out - np.linspace(np.pi / Nc, 2 * np.pi / Nc, nb_fps // 2), # tilt - [None] * (nb_fps // 4), # None + np.around(np.linspace(1, Nc, 4 * nb_frames)).astype(int), # Nc + np.linspace(2 * np.pi / Nc, np.pi / Nc, 2 * nb_frames), # tilt + np.around(np.sin(np.linspace(0, 2 * np.pi, 2 * nb_frames))).astype(bool), # in_out + np.linspace(np.pi / Nc, 2 * np.pi / Nc, 2 * nb_frames), # tilt + [None] * nb_frames, # None # Spiral - np.linspace(0, np.sqrt(3), nb_fps // 2) ** 2, # spiral - np.linspace(3, 1, nb_fps // 3), # spiral - np.linspace(1, 2, nb_fps // 2), # nb_revolutions - np.linspace(2, 0, nb_fps), # nb_revolutions - [None] * (nb_fps // 4), # None + np.linspace(1e-5, 1, 2 * nb_frames), # nb_revolutions + np.linspace(1, np.sqrt(1 / 3), 2 * nb_frames) ** 2, # spiral + np.linspace(1 / 3, 1, 2 * nb_frames), # spiral + np.linspace(1, 3, 2 * nb_frames), # nb_revolutions + np.linspace(3, 1e-5, 4 * nb_frames), # nb_revolutions + [None] * nb_frames, # None # Cones - np.linspace(0, 5, nb_fps // 2), # nb_zigzags - np.linspace(1, 2, nb_fps // 4), # width - np.linspace(2, 0, nb_fps // 2), # width - [None] * (nb_fps // 4), # None + np.linspace(0, 5, 2 * nb_frames), # nb_zigzags + np.linspace(1, 2, nb_frames), # width + np.linspace(2, 0, 2 * nb_frames), # width + [None] * nb_frames, # None # Sinusoids - np.linspace(0, 1, nb_fps // 2), # width & nb_zigzags - np.linspace(1, 0, nb_fps // 2), # width & nb_zigzags - [None] * (nb_fps // 4), # None + np.linspace(0, 1, 2 * nb_frames), # width & nb_zigzags + np.linspace(1, 0, 2 * nb_frames), # width & nb_zigzags + [None] * nb_frames, # None # Rings - np.around(np.linspace(1, nb_repetitions, nb_fps)).astype(int), # Nc & nb_rings - np.around(np.linspace(nb_repetitions, Nc, nb_fps // 2)).astype(int), # Nc - [None] * (nb_fps // 2), # None + np.around(np.linspace(1, nb_repetitions, 4 * nb_frames)).astype( + int + ), # Nc & nb_rings + np.around(np.linspace(nb_repetitions, Nc, 2 * nb_frames)).astype(int), # Nc + [None] * nb_frames, # None # Rosette - [None] * (nb_fps // 2), # None - np.around(np.linspace(0, 30, nb_fps)).astype(int), # coprime_index - [None] * (nb_fps // 2), # None + np.around(np.linspace(0, np.sqrt(30), 6 * nb_frames) ** 2).astype( + int + ), # coprime_index + [None] * nb_frames, # None # Waves - np.linspace(0, 2, nb_fps // 2), # width & nb_zigzags - np.linspace(2, 1, nb_fps // 2), # width & nb_zigzags - [None] * (nb_fps // 4), # None + np.linspace(0, 2, 4 * nb_frames), # width & nb_zigzags + np.linspace(2, 1, 2 * nb_frames), # width & nb_zigzags + [None] * nb_frames, # None # Lissajous - [None] * (nb_fps // 2), # None - np.linspace(1, 10, 2 * nb_fps), # density - [None] * (nb_fps // 2), # None + np.linspace(1, 10, 6 * nb_frames), # density + [None] * nb_frames, # None ] -# %% - frame_setup = [ (f, i, name, arg) @@ -180,7 +170,6 @@ def draw_frame(func, index, name, arg, save_dir="/tmp/"): ) -# %% # Make a GIF of all images. imgs = [Image.open(img) for img in image_files] imgs[0].save( @@ -188,7 +177,7 @@ def draw_frame(func, index, name, arg, save_dir="/tmp/"): save_all=True, append_images=imgs[1:], optimize=False, - duration=200, + duration=duration, loop=0, ) diff --git a/examples/example_gif_3D.py b/examples/example_gif_3D.py index 1b5e5cfc..4965b22d 100644 --- a/examples/example_gif_3D.py +++ b/examples/example_gif_3D.py @@ -7,97 +7,45 @@ """ -# %% +import joblib import matplotlib.pyplot as plt -import mrinufft.trajectories.display as mtd -import mrinufft.trajectories.trajectory3D as mtt import numpy as np -import joblib from PIL import Image, ImageSequence + +import mrinufft.trajectories.display as mtd +import mrinufft.trajectories.trajectory3D as mtt from mrinufft.trajectories.display import displayConfig -from tqdm.notebook import tqdm -# %% [markdown] -# # Options +# Options -# %% Nc = 8 * 8 Ns = 200 - nb_repetitions = 8 -# %% one_shot = 0 figsize = 5 -nb_fps = 12 -name_slow = 1 +nb_frames = 3 +duration = 150 # seconds -# %% [markdown] -# # Generation -# %% +# Generation + # Initialize trajectory function functions = [ # 3D Cones - ("3D Cones", lambda x: mtt.initialize_3D_cones(Nc, Ns, nb_zigzags=0)[::-1]), - ("3D Cones", lambda x: mtt.initialize_3D_cones(Nc, Ns, nb_zigzags=x)[::-1]), ("3D Cones", lambda x: mtt.initialize_3D_cones(Nc, Ns, width=x)[::-1]), ("3D Cones", lambda x: mtt.initialize_3D_cones(Nc, Ns, width=x)[::-1]), + ("3D Cones", lambda x: mtt.initialize_3D_cones(Nc, Ns, nb_zigzags=x)[::-1]), + ("3D Cones", lambda x: mtt.initialize_3D_cones(Nc, Ns, nb_zigzags=x)[::-1]), ("3D Cones", lambda x: mtt.initialize_3D_cones(Nc, Ns)[::-1]), # FLORET - ( - "FLORET", - lambda x: mtt.initialize_3D_floret( - Nc // nb_repetitions, - Ns, - nb_revolutions=1, - max_angle=np.pi / 2, - axes=[0], - cone_tilt=None, - ), - ), - ( - "FLORET", - lambda x: mtt.initialize_3D_floret( - x * Nc // nb_repetitions, - Ns, - nb_revolutions=x, - max_angle=np.pi / 2, - axes=[0], - cone_tilt=None, - ), - ), - ( - "FLORET", - lambda x: mtt.initialize_3D_floret( - Nc, Ns, nb_revolutions=nb_repetitions, max_angle=x, axes=[0], cone_tilt=None - ), - ), - ( - "FLORET", - lambda x: mtt.initialize_3D_floret( - Nc, Ns, nb_revolutions=nb_repetitions, max_angle=x, axes=[0], cone_tilt=None - ), - ), - # ("FLORET", lambda x: mtt.initialize_3D_floret(Nc, Ns, nb_cones=nb_repetitions, max_angle=np.pi / 2, nb_revolutions=x / 2, spiral=x, axes=[0], cone_tilt=None)), - # ("FLORET", lambda x: mtt.initialize_3D_floret(Nc, Ns, nb_cones=nb_repetitions, max_angle=np.pi / 2, nb_revolutions=x / 2, spiral=x, axes=[0], cone_tilt=None)), - ( - "FLORET", - lambda x: mtt.initialize_3D_floret( - Nc, - Ns, - nb_revolutions=nb_repetitions, - max_angle=np.pi / 2, - axes=[0], - cone_tilt=None, - ), - ), + ("FLORET", lambda x: mtt.initialize_3D_floret(Nc, Ns, nb_revolutions=x)), + ("FLORET", lambda x: mtt.initialize_3D_floret(Nc, Ns, nb_revolutions=x)), + ("FLORET", lambda x: mtt.initialize_3D_floret(Nc, Ns, max_angle=x)), + ("FLORET", lambda x: mtt.initialize_3D_floret(Nc, Ns, max_angle=x)), + ("FLORET", lambda x: mtt.initialize_3D_floret(Nc, Ns)), # Seiffert spirals - ( - "Seiffert spiral / Yarnball", - lambda x: mtt.initialize_3D_seiffert_spiral(Nc, Ns, curve_index=0), - ), ( "Seiffert spiral / Yarnball", lambda x: mtt.initialize_3D_seiffert_spiral(Nc, Ns, curve_index=x), @@ -105,19 +53,19 @@ ( "Seiffert spiral / Yarnball", lambda x: mtt.initialize_3D_seiffert_spiral( - Nc, Ns, curve_index=0.9, nb_revolutions=x + Nc, Ns, curve_index=0.7, nb_revolutions=x ), ), ( "Seiffert spiral / Yarnball", lambda x: mtt.initialize_3D_seiffert_spiral( - Nc, Ns, curve_index=0.9, nb_revolutions=x + Nc, Ns, curve_index=0.7, nb_revolutions=x ), ), ( "Seiffert spiral / Yarnball", lambda x: mtt.initialize_3D_seiffert_spiral( - Nc, Ns, curve_index=0.9, nb_revolutions=1 + Nc, Ns, curve_index=0.7, nb_revolutions=1 ), ), # Helical shells @@ -139,18 +87,6 @@ Nc, Ns, nb_shells=nb_repetitions, spiral_reduction=3 )[::-1], ), - # Annular shells - # ("Annular shells", lambda x: mtt.initialize_3D_annular_shells(x * Nc // nb_repetitions, Ns, nb_shells=x, ring_tilt=0)[::-1]), - # ("Annular shells", lambda x: mtt.initialize_3D_annular_shells(Nc, Ns, nb_shells=nb_repetitions, ring_tilt=x)[::-1]), - # ("Annular shells", lambda x: mtt.initialize_3D_annular_shells(Nc, Ns, nb_shells=nb_repetitions, ring_tilt=x)[::-1]), - # ("Annular shells", lambda x: mtt.initialize_3D_annular_shells(Nc, Ns, nb_shells=nb_repetitions, ring_tilt=np.pi / 2)[::-1]), - # Seiffert shells - # ("Seiffert shells", lambda x: mtt.initialize_3D_seiffert_shells(x * Nc // nb_repetitions, Ns, nb_shells=x, curve_index=0)[::-1]), - # ("Seiffert shells", lambda x: mtt.initialize_3D_seiffert_shells(Nc, Ns, nb_shells=nb_repetitions, curve_index=x)[::-1]), - # ("Seiffert shells", lambda x: mtt.initialize_3D_seiffert_shells(Nc, Ns, nb_shells=nb_repetitions, curve_index=x)[::-1]), - # ("Seiffert shells", lambda x: mtt.initialize_3D_seiffert_shells(Nc, Ns, nb_shells=nb_repetitions, nb_revolutions=x)[::-1]), - # ("Seiffert shells", lambda x: mtt.initialize_3D_seiffert_shells(Nc, Ns, nb_shells=nb_repetitions, nb_revolutions=x)[::-1]), - # ("Seiffert shells", lambda x: mtt.initialize_3D_seiffert_shells(Nc, Ns, nb_shells=nb_repetitions, nb_revolutions=1)[::-1]), # Wave-CAIPI ( "Wave-CAIPI", @@ -164,58 +100,38 @@ 2 * Nc, Ns, nb_revolutions=5 * x, width=x ), ), - # ("Wave-CAIPI", lambda x: mtt.initialize_3D_wave_caipi(2 * Nc, Ns, packing=["triangular", "square", "circular", "triangular"][x])), ("Wave-CAIPI", lambda x: mtt.initialize_3D_wave_caipi(2 * Nc, Ns)), ] -# %% # Initialize trajectory arguments arguments = [ # 3D Cones - [None] * (nb_fps // 4), # None - np.linspace(0, np.sqrt(5), nb_fps) ** 2, # nb_zigzags - np.linspace(1, 2, nb_fps // 2), # width - np.linspace(2, 1, nb_fps // 2), # width - [None] * (nb_fps // 4), # None + np.linspace(0, 2, 4 * nb_frames), # width + np.linspace(2, 1, 2 * nb_frames), # width + np.linspace(np.sqrt(5), 1, 4 * nb_frames) ** 2, # nb_zigzags + np.linspace(1, np.sqrt(5), 2 * nb_frames) ** 2, # nb_zigzags + [None] * nb_frames, # None # FLORET - [None] * (nb_fps // 4), # None - np.around(np.linspace(1, nb_repetitions, nb_fps)).astype(int), # nb_cones - np.linspace(np.pi / 2, np.pi / 4, nb_fps // 2), # max_angle - np.linspace(np.pi / 4, np.pi / 2, nb_fps // 2), # max_angle - # np.linspace(2, 1, nb_fps // 2), # spiral & nb_revolutions - # np.linspace(1, 2, nb_fps // 2), # spiral & nb_revolutions - [None] * (nb_fps // 4), # None + np.linspace(1, 3, 4 * nb_frames), # nb_revolutions + np.linspace(3, 1, 2 * nb_frames), # nb_revolutions + np.linspace(np.pi / 2, np.pi / 4, 4 * nb_frames), # max_angle + np.linspace(np.pi / 4, np.pi / 2, 2 * nb_frames), # max_angle + [None] * nb_frames, # None # Seiffert spiral - [None] * (nb_fps // 4), # None - np.linspace(0, 0.9, nb_fps), # curve_index - np.linspace(1, 3, nb_fps), # nb_revolutions - np.linspace(3, 1, nb_fps // 2), # nb_revolutions - [None] * (nb_fps // 4), # None + np.linspace(0, 0.7, 4 * nb_frames), # curve_index + np.linspace(1, 2, 4 * nb_frames), # nb_revolutions + np.linspace(2, 1, 2 * nb_frames), # nb_revolutions + [None] * nb_frames, # None # Helical shells - np.around(np.linspace(1, nb_repetitions, nb_fps)).astype(int), # nb_cones - np.linspace(1, 3, nb_fps), # spiral_reduction - [None] * (nb_fps // 4), # None - # Annular shells - # np.around(np.linspace(1, nb_repetitions, nb_fps)).astype(int), # nb_cones - # np.linspace(0, np.pi, nb_fps), # ring_tilt - # np.linspace(np.pi, np.pi / 2, nb_fps // 2), # ring_tilt - # [None] * (nb_fps // 4), # None - # Seiffert shells - # np.around(np.linspace(1, nb_repetitions, nb_fps)).astype(int), # nb_cones - # np.linspace(0, 0.95, nb_fps), # curve_index - # np.linspace(0.95, 0.5, nb_fps), # curve_index - # np.linspace(1, 3, nb_fps), # nb_revolutions - # np.linspace(3, 1, nb_fps // 2), # nb_revolutions - # [None] * (nb_fps // 4), # None + np.around(np.linspace(1, nb_repetitions, 4 * nb_frames)).astype(int), # nb_cones + np.linspace(1, 3, 4 * nb_frames), # spiral_reduction + [None] * nb_frames, # None # Wave-CAIPI - np.linspace(0, 2, nb_fps), # nb_revolutions & width - np.linspace(2, 1, nb_fps // 2), # nb_revolutions & width - # np.around(np.linspace(0, 3, nb_fps)).astype(int), # packing - [None] * (nb_fps // 4), # None + np.linspace(0, 2, 4 * nb_frames), # nb_revolutions & width + np.linspace(2, 1, 2 * nb_frames), # nb_revolutions & width + [None] * nb_frames, # None ] -# %% - frame_setup = [ (f, i, name, arg) @@ -259,7 +175,6 @@ def draw_frame(func, index, name, arg, save_dir="/tmp/"): ) -# %% # Make a GIF of all images. imgs = [Image.open(img) for img in image_files] imgs[0].save( @@ -267,7 +182,7 @@ def draw_frame(func, index, name, arg, save_dir="/tmp/"): save_all=True, append_images=imgs[1:], optimize=False, - duration=200, + duration=duration, loop=0, ) diff --git a/examples/example_readme.py b/examples/example_readme.py index 1b9b58d8..d1341f26 100644 --- a/examples/example_readme.py +++ b/examples/example_readme.py @@ -6,11 +6,12 @@ """ import matplotlib.pyplot as plt -from scipy.datasets import face import numpy as np +from scipy.datasets import face + import mrinufft -from mrinufft.trajectories import display from mrinufft.density import voronoi +from mrinufft.trajectories import display # Create a 2D Radial trajectory for demo samples_loc = mrinufft.initialize_2D_radial(Nc=100, Ns=500) diff --git a/examples/example_stacked.py b/examples/example_stacked.py index d7f5f746..b70a31f8 100644 --- a/examples/example_stacked.py +++ b/examples/example_stacked.py @@ -13,6 +13,7 @@ import matplotlib.pyplot as plt import numpy as np + from mrinufft import display_2D_trajectory plt.rcParams["image.cmap"] = "gray" diff --git a/examples/example_trajectory_tools.py b/examples/example_trajectory_tools.py index a8f08a60..5cdab60c 100644 --- a/examples/example_trajectory_tools.py +++ b/examples/example_trajectory_tools.py @@ -28,8 +28,7 @@ # Internal import mrinufft as mn import mrinufft.trajectories.tools as tools - -from mrinufft import displayConfig, display_2D_trajectory, display_3D_trajectory +from mrinufft import display_2D_trajectory, display_3D_trajectory, displayConfig from mrinufft.trajectories.utils import KMAX diff --git a/src/mrinufft/__init__.py b/src/mrinufft/__init__.py index f468a91b..a74a28c3 100644 --- a/src/mrinufft/__init__.py +++ b/src/mrinufft/__init__.py @@ -15,6 +15,7 @@ from .trajectories import ( initialize_2D_radial, initialize_2D_spiral, + initialize_2D_fibonacci_spiral, initialize_2D_cones, initialize_2D_sinusoide, initialize_2D_propeller, @@ -51,6 +52,7 @@ "list_backends", "initialize_2D_radial", "initialize_2D_spiral", + "initialize_2D_fibonacci_spiral", "initialize_2D_cones", "initialize_2D_sinusoide", "initialize_2D_propeller", @@ -75,6 +77,7 @@ "shellify", "duplicate_along_axes", "radialize_center", + "patch_center_anomaly", "displayConfig", "display_2D_trajectory", "display_3D_trajectory", diff --git a/src/mrinufft/trajectories/__init__.py b/src/mrinufft/trajectories/__init__.py index 5df218cf..11a2e07d 100644 --- a/src/mrinufft/trajectories/__init__.py +++ b/src/mrinufft/trajectories/__init__.py @@ -3,6 +3,7 @@ from .trajectory2D import ( initialize_2D_radial, initialize_2D_spiral, + initialize_2D_fibonacci_spiral, initialize_2D_cones, initialize_2D_sinusoide, initialize_2D_propeller, @@ -40,9 +41,12 @@ display_3D_trajectory, ) +from .gradients import patch_center_anomaly + __all__ = [ "initialize_2D_radial", "initialize_2D_spiral", + "initialize_2D_fibonacci_spiral", "initialize_2D_cones", "initialize_2D_sinusoide", "initialize_2D_propeller", @@ -70,4 +74,5 @@ "displayConfig", "display_2D_trajectory", "display_3D_trajectory", + "patch_center_anomaly", ] diff --git a/src/mrinufft/trajectories/display.py b/src/mrinufft/trajectories/display.py index bc924bd4..50f20344 100644 --- a/src/mrinufft/trajectories/display.py +++ b/src/mrinufft/trajectories/display.py @@ -1,4 +1,4 @@ -"""Display function for trajectories.""" +"""Display functions for trajectories.""" import itertools diff --git a/src/mrinufft/trajectories/gradients.py b/src/mrinufft/trajectories/gradients.py new file mode 100644 index 00000000..39ac934c --- /dev/null +++ b/src/mrinufft/trajectories/gradients.py @@ -0,0 +1,118 @@ +"""Functions to improve/modify gradients.""" + +import numpy as np +import numpy.linalg as nl +from scipy.interpolate import CubicSpline + + +def patch_center_anomaly( + shot_or_params, + update_shot=None, + update_parameters=None, + in_out=False, + learning_rate=1e-1, +): + """Re-position samples to avoid center anomalies. + + Some trajectories behave slightly differently from expected when + approaching definition bounds, most often the k-space center as + for spirals in some cases. + + This function enforces non-strictly increasing monoticity of + sample distances from the center, effectively reducing slew + rates and smoothing gradient transitions locally. + + Shots can be updated with provided functions to keep fitting + their strict definition, or it can be smoothed using cubic + spline approximations. + + Parameters + ---------- + shot_or_params : np.array, list + Either a single shot of shape (Ns, Nd), or a list of arbitrary + arguments used by ``update_shot`` to initialize a single shot. + update_shot : function, optional + Function used to initialize a single shot based on parameters + provided by ``update_parameters``. If None, cubic splines are + used as an approximation instead, by default None + update_parameters : function, optional + Function used to update shot parameters when provided in + ``shot_or_params`` from an updated shot and parameters. + If None, cubic spline parameterization is used instead, + by default None + in_out : bool, optional + Whether the shot is going in-and-out or start from the center, + by default False + learning_rate : float, optional + Learning rate used in the iterative optimization process, + by default 1e-1 + + Returns + ------- + array_like + N-D trajectory based on ``shot_or_params`` if a shot or + update_shot otherwise. + list + Updated parameters either in the ``shot_or_params`` format + if params, or cubic spline parameterization as an array of + floats between 0 and 1 otherwise. + """ + if update_shot is None or update_parameters is None: + single_shot = shot_or_params + else: + parameters = shot_or_params + single_shot = update_shot(*parameters) + Ns = single_shot.shape[0] + + old_index = 0 + old_x_axis = np.zeros(Ns) + x_axis = np.linspace(0, 1, Ns) + + if update_shot is None or update_parameters is None: + + def _default_update_parameters(shot, *parameters): + return parameters + + update_parameters = _default_update_parameters + parameters = [x_axis] + update_shot = CubicSpline(x_axis, single_shot) + + while nl.norm(x_axis - old_x_axis) / Ns > 1e-10: + # Determine interval to fix + single_shot = update_shot(*parameters) + gradient_norm = nl.norm( + np.diff(single_shot[(Ns // 2) * in_out :], axis=0), axis=-1 + ) + arc_length = np.cumsum(gradient_norm) + + index = gradient_norm[1:] - arc_length[1:] / ( + 1 + np.arange(1, len(gradient_norm)) + ) + index = np.where(index < 0, np.inf, index) + index = max(old_index, 2 + np.argmin(index)) + start_index = (Ns // 2 + (Ns % 2) - index) * in_out + final_index = (Ns // 2) * in_out + index + + # Balance over interval + cbs = CubicSpline(x_axis, single_shot.T, axis=1) + gradient_norm = nl.norm( + np.diff(single_shot[start_index:final_index], axis=0), axis=-1 + ) + + old_x_axis = np.copy(x_axis) + new_x_axis = np.cumsum(np.diff(x_axis[start_index:final_index]) / gradient_norm) + new_x_axis *= (x_axis[final_index - 1] - x_axis[start_index]) / new_x_axis[ + -1 + ] # Normalize + new_x_axis += x_axis[start_index] + x_axis[start_index + 1 : final_index - 1] += ( + new_x_axis[:-1] - x_axis[start_index + 1 : final_index - 1] + ) * learning_rate + + # Update loop variables + old_index = index + single_shot = cbs(x_axis).T + parameters = update_parameters(single_shot, *parameters) + + single_shot = single_shot = update_shot(*parameters) + return single_shot, parameters diff --git a/src/mrinufft/trajectories/maths.py b/src/mrinufft/trajectories/maths.py index d8ad192f..a413df2c 100644 --- a/src/mrinufft/trajectories/maths.py +++ b/src/mrinufft/trajectories/maths.py @@ -210,6 +210,55 @@ def Ra(vector, theta): ############# +def is_from_fibonacci_sequence(n): + """Check if an integer belongs to the Fibonacci sequence. + + An integer belongs to the Fibonacci sequence if either + :math:`5*n²+4` or :math:`5*n²-4` is a perfect square + (`Wikipedia `_). + + Parameters + ---------- + n : int + Integer to check. + + Returns + ------- + bool + Whether or not ``n`` belongs to the Fibonacci sequence. + """ + + def _is_perfect_square(n): + r = int(np.sqrt(n)) + return r * r == n + + return _is_perfect_square(5 * n**2 + 4) or _is_perfect_square(5 * n**2 - 4) + + +def get_closest_fibonacci_number(x): + """Provide the closest Fibonacci number. + + Parameters + ---------- + x : float + Value to match. + + Returns + ------- + int + Closest number from the Fibonacci sequence. + """ + # Find the power such that x ~= phi ** power + phi = (1 + np.sqrt(5)) / 2 + power = np.ceil(np.log(x) / np.log(phi)) + 1 + + # Check closest between the ones below and above n + lower_xf = int(np.around(phi ** (power) / np.sqrt(5))) + upper_xf = int(np.around(phi ** (power + 1) / np.sqrt(5))) + xf = lower_xf if (x - lower_xf) < (upper_xf - x) else upper_xf + return xf + + def generate_fibonacci_lattice(nb_points, epsilon=0.25): """Generate 2D Cartesian coordinates using the Fibonacci lattice. diff --git a/src/mrinufft/trajectories/tools.py b/src/mrinufft/trajectories/tools.py index a3feebcd..316190e9 100644 --- a/src/mrinufft/trajectories/tools.py +++ b/src/mrinufft/trajectories/tools.py @@ -1,12 +1,9 @@ -"""Functions to manipulate trajectories.""" +"""Functions to manipulate/modify trajectories.""" import numpy as np from .maths import Rv, Rx, Ry, Rz -from .utils import ( - KMAX, - initialize_tilt, -) +from .utils import KMAX, initialize_tilt ################ # DIRECT TOOLS # diff --git a/src/mrinufft/trajectories/trajectory2D.py b/src/mrinufft/trajectories/trajectory2D.py index bf0193ec..f87c6362 100644 --- a/src/mrinufft/trajectories/trajectory2D.py +++ b/src/mrinufft/trajectories/trajectory2D.py @@ -1,18 +1,13 @@ -"""2D trajectory initializations.""" +"""Functions to initialize 2D trajectories.""" import numpy as np +import numpy.linalg as nl +from scipy.interpolate import CubicSpline -from .maths import ( - R2D, - compute_coprime_factors, -) -from .utils import ( - KMAX, - initialize_tilt, - initialize_spiral, -) +from .gradients import patch_center_anomaly +from .maths import R2D, compute_coprime_factors, is_from_fibonacci_sequence from .tools import rotate - +from .utils import KMAX, initialize_algebraic_spiral, initialize_tilt ##################### # CIRCULAR PATTERNS # @@ -47,9 +42,20 @@ def initialize_2D_radial(Nc, Ns, tilt="uniform", in_out=False): def initialize_2D_spiral( - Nc, Ns, tilt="uniform", in_out=False, nb_revolutions=1, spiral="archimedes" + Nc, + Ns, + tilt="uniform", + in_out=False, + nb_revolutions=1, + spiral="archimedes", + patch_center=True, ): - """Initialize a 2D spiral trajectory. + """Initialize a 2D algebraic spiral trajectory. + + A generalized function that generates algebraic spirals defined + through the :math:`r = a O^n` equality, with :math:`r` the radius, + :math:`O` the polar angle and :math:`n` the spiral power. + Common algebraic spirals include Archimedes, Fermat and Galilean spirals. Parameters ---------- @@ -64,27 +70,141 @@ def initialize_2D_spiral( nb_revolutions : int, optional Number of revolutions, by default 1 spiral : str, float, optional - Spiral type, by default "archimedes" + Spiral type or algebraic power, by default "archimedes" + patch_center : bool, optional + Whether the spiral anomaly at the center should be patched + or not for spirals with `spiral` :math:`>2`, by default True Returns ------- array_like 2D spiral trajectory + + Raises + ------ + ValueError + If `spiral` is negative. + + Notes + ----- + Algebraic spirals with negative powers, like hyperbolic or + lithuus spirals, show asymptotic behaviors around the center. + It makes them irrelevant for MRI and therefore negative powers + are not allowed as an argument. """ - # Initialize a first shot in polar coordinates - segment = np.linspace(-1 if (in_out) else 0, 1, Ns) - radius = KMAX * segment - angles = 2 * np.pi * nb_revolutions * (np.abs(segment) ** initialize_spiral(spiral)) + # Check spiral power is not negative + spiral_power = initialize_algebraic_spiral(spiral) + if spiral_power <= 0: + raise ValueError(f"Negative spiral definition is invalid (spiral={spiral}).") - # Convert the first shot to Cartesian coordinates - trajectory = np.zeros((Nc, Ns, 2)) - trajectory[0, :, 0] = radius * np.cos(angles) - trajectory[0, :, 1] = radius * np.sin(angles) + # Initialize a first shot in polar coordinates + angles = 2 * np.pi * nb_revolutions * np.linspace(-1 if (in_out) else 0, 1, Ns) + radius = np.abs(angles) ** spiral_power + + # Algebraic spirals with power coefficients superior to 1 + # have a non-monotonic gradient norm when varying the angle + # over [0, +inf) + def _update_shot(angles, radius, *args): + shot = np.sign(angles) * np.abs(radius) * np.exp(1j * np.abs(angles)) + return np.stack([shot.real, shot.imag], axis=-1) + + def _update_parameters(single_shot, angles, radius, spiral_power): + radius = nl.norm(single_shot, axis=-1) + angles = np.sign(angles) * np.abs(radius) ** (1 / spiral_power) + return angles, radius, spiral_power + + if spiral_power < 1 and patch_center: + parameters = (angles, radius, spiral_power) + learning_rate = min( + 1, spiral_power + ) # because low spiral power requires higher accuracy + _, parameters = patch_center_anomaly( + parameters, + update_shot=_update_shot, + update_parameters=_update_parameters, + in_out=in_out, + learning_rate=learning_rate, + ) + angles, radius, _ = parameters + + # Convert the first shot from polar to Cartesian coordinates + trajectory = np.zeros((Nc, len(angles), 2)) + trajectory[0, :] = _update_shot(angles, radius) # Rotate the first shot Nc times rotation = R2D(initialize_tilt(tilt, Nc) / (1 + in_out)).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_fibonacci_spiral(Nc, Ns, spiral_reduction=1, patch_center=True): + """Initialize a 2D Fibonacci spiral trajectory. + + A non-algebraic spiral trajectory based on the Fibonacci sequence, + reproducing the proposition from [CA99]_ in order to generate + a uniform distribution with center-out shots. + + The number of shots is required to belong to the Fibonacci + sequence for the trajectory definition to be relevant. + + Parameters + ---------- + Nc : int + Number of shots + Ns : int + Number of samples per shot + spiral_reduction : float, optional + Factor used to reduce the automatic spiral length, by default 1 + patch_center : bool, optional + Whether the spiral anomaly at the center should be patched + or not, by default True + + Returns + ------- + array_like + 2D Fibonacci spiral trajectory + + References + ---------- + .. [CA99] Cline, Harvey E., and Thomas R. Anthony. + "Uniform k-space sampling with an interleaved + Fibonacci spiral acquisition." + In Proceedings of the 7th Annual Meeting of ISMRM, + Philadelphia, USA, vol. 1657. 1999. + """ + # Check if Nc is in the Fibonacci sequence + if not is_from_fibonacci_sequence(Nc): + raise ValueError("Nc should belong to the Fibonacci sequence.") + + # Initialize all shots + Ns_reduced = int(np.around(Ns / spiral_reduction)) + inter_range = np.arange(Nc).reshape((-1, 1)) + intra_range = np.arange(Ns_reduced).reshape((1, -1)) + phi_bonacci = (np.sqrt(5) - 1) / 2 + radius = np.sqrt((intra_range + (inter_range / Nc)) / (Nc * Ns_reduced)) + angles = 2j * np.pi * phi_bonacci * np.around(Nc * intra_range + inter_range) + trajectory = radius * np.exp(angles) + + # Put Ns samples along reduced spirals if relevant + if spiral_reduction != 1: + reduced_x_axis = np.linspace(0, 1, Ns_reduced) + normal_x_axis = np.linspace(0, 1, Ns) + cbs = CubicSpline(reduced_x_axis, trajectory, axis=1) + trajectory = cbs(normal_x_axis) + + # Normalize and reformat trajectory + trajectory *= KMAX / np.max(np.abs(trajectory)) + trajectory = np.stack([trajectory.real, trajectory.imag], axis=-1) + + # Patch center anomaly if requested + if patch_center: + patched_trajectory = [] + for i in range(Nc): + patched_shot, _ = patch_center_anomaly(trajectory[i], in_out=False) + patched_trajectory.append(patched_shot) + trajectory = np.array(patched_trajectory) return trajectory diff --git a/src/mrinufft/trajectories/trajectory3D.py b/src/mrinufft/trajectories/trajectory3D.py index aab30197..b67d7833 100644 --- a/src/mrinufft/trajectories/trajectory3D.py +++ b/src/mrinufft/trajectories/trajectory3D.py @@ -1,16 +1,15 @@ -"""3D Trajectory initialization functions.""" +"""Functions to initialize 3D trajectories.""" + +from functools import partial import numpy as np import numpy.linalg as nl - -from functools import partial from scipy.special import ellipj, ellipk -from .maths import Ry, Rz, Ra, generate_fibonacci_circle, CIRCLE_PACKING_DENSITY -from .tools import precess, conify, duplicate_along_axes +from .maths import CIRCLE_PACKING_DENSITY, Ra, Ry, Rz, generate_fibonacci_circle +from .tools import conify, duplicate_along_axes, precess from .trajectory2D import initialize_2D_spiral -from .utils import initialize_tilt, initialize_shape_norm, KMAX, Packings - +from .utils import KMAX, Packings, initialize_shape_norm, initialize_tilt ############################ # FREEFORM 3D TRAJECTORIES # @@ -411,8 +410,8 @@ def initialize_3D_helical_shells( Number of samples per shot nb_shells : int Number of concentric shells/spheres - spiral_reduction : float - Factor used to reduce the automatic spiral curvature, by default 1 + spiral_reduction : float, optional + Factor used to reduce the automatic spiral length, by default 1 shell_tilt : str, float, optional Angle between consecutive shells along z-axis, by default "intergaps" shot_tilt : str, float, optional diff --git a/src/mrinufft/trajectories/utils.py b/src/mrinufft/trajectories/utils.py index ec4c1008..748a5270 100644 --- a/src/mrinufft/trajectories/utils.py +++ b/src/mrinufft/trajectories/utils.py @@ -1,4 +1,4 @@ -"""Utility functions for the trajectory design.""" +"""Utility functions in general.""" from enum import Enum, EnumMeta from numbers import Real @@ -66,14 +66,16 @@ class Gammas(FloatEnum): class Spirals(FloatEnum): - """Enumerate spiral types.""" + """Enumerate algebraic spiral types.""" ARCHIMEDES = 1 ARITHMETIC = ARCHIMEDES - FERMAT = 2 + GALILEAN = 2 + GALILEO = GALILEAN + FERMAT = 0.5 PARABOLIC = FERMAT HYPERBOLIC = -1 - LITHUUS = -2 + LITHUUS = -0.5 class NormShapes(FloatEnum): @@ -485,8 +487,8 @@ def initialize_tilt(tilt, nb_partitions=1): raise NotImplementedError(f"Unknown tilt name: {tilt}") -def initialize_spiral(spiral): - """Initialize the spiral type. +def initialize_algebraic_spiral(spiral): + """Initialize the algebraic spiral type. Parameters ----------