Spiral patch (#105)

* Change spiral parameterization, add method to patch abnormal gradients

* Add 2D Fibonacci spiral, add examples for all spirals, improve gradient anomaly fix

* Apply ruff/black formatting

* Fix null power case

* Update and clean gif examples

* Make finding closest Fibonacci constant time search, update documentation

* Apply black & ruff formatting

---------

Co-authored-by: Guillaume DAVAL-FREROT <gdavalfrerot@gmail.com>

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")
 

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
 # -----

examples/example_3D_trajectories.py

@@ -29,7 +29,6 @@
 
 # Internal
 import mrinufft as mn
-
 from mrinufft import display_2D_trajectory, display_3D_trajectory
 
 

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
 # ------------------

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

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,63 +69,61 @@
     ("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
     ("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 * 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,15 +170,14 @@ 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(
     "mrinufft_2D_traj.gif",
     save_all=True,
     append_images=imgs[1:],
     optimize=False,
-    duration=200,
+    duration=duration,
     loop=0,
 )