Propeller

examples/example_2D_trajectories.py

@@ -67,6 +67,7 @@ def show_trajectory(trajectory, one_shot, figure_size):
 Ns = 256  # Number of samples per shot
 in_out = True  # Choose between in-out or center-out trajectories
 tilt = "uniform"  # Choose the angular distance between shots
+nb_repetitions = 6  # Number of strips when relevant
 
 # Display parameters
 figure_size = 6  # Figure size for trajectory plots
@@ -301,6 +302,46 @@ def show_trajectory(trajectory, one_shot, figure_size):
 show_trajectory(trajectory, figure_size=figure_size, one_shot=one_shot)
 
 
+# %%
+# PROPELLER
+# ---------
+#
+# The PROPELLER trajectory is generally used along a specific
+# reconstruction pipeline described in [Pip99]_ to correct for
+# motion artifacts.
+#
+# The acronym PROPELLER stands for Periodically Rotated
+# Overlapping ParallEL Lines with Enhanced Reconstruction,
+# and the method is also commonly known under other aliases
+# depending on the vendor, with some variations: BLADE,
+# MulitVane, RADAR, JET.
+#
+# Arguments:
+#
+# - ``Nc (int)``: number of individual shots. See radial
+# - ``Ns (int)``: number of samples per shot. See radial
+# - ``nb_strips (int)``: number of strips covering the k-space.
+#   ``(default "uniform")``. See radial
+#
+
+trajectory = mn.initialize_2D_propeller(Nc, Ns, nb_strips=nb_repetitions)
+show_trajectory(trajectory, figure_size=figure_size, one_shot=one_shot)
+
+
+# %%
+# ``nb_strips (int)``
+# ~~~~~~~~~~~~~~~~~~~
+#
+# The number of individual strips dividing the k-space circle. It must divide
+# the number of shots ``Nc``, and it is recommended to choose it such that the
+# ratio is even to cover the center.
+#
+
+arguments = [2, 3, 4, 6]
+function = lambda x: mn.initialize_2D_propeller(Nc, Ns, nb_strips=x)
+show_argument(function, arguments, one_shot=one_shot, subfigure_size=subfigure_size)
+
+
 # %%
 # Rings
 # -------
@@ -565,3 +606,12 @@ def show_trajectory(trajectory, one_shot, figure_size):
 arguments = [1, 1.5, 2, 3]
 function = lambda x: mn.initialize_2D_lissajous(Nc, Ns, density=x)
 show_argument(function, arguments, one_shot=one_shot, subfigure_size=subfigure_size)
+
+
+# %%
+# References
+# ==========
+#
+# .. [Pip99] Pipe, James G. "Motion correction with PROPELLER MRI:
+#    application to head motion and free‐breathing cardiac imaging."
+#    Magnetic Resonance in Medicine 42, no. 5 (1999): 963-969.

src/mrinufft/__init__.py

@@ -17,6 +17,7 @@
     initialize_2D_spiral,
     initialize_2D_cones,
     initialize_2D_sinusoide,
+    initialize_2D_propeller,
     initialize_2D_rosette,
     initialize_2D_rings,
     initialize_2D_polar_lissajous,
@@ -42,6 +43,7 @@
     "initialize_2D_spiral",
     "initialize_2D_cones",
     "initialize_2D_sinusoide",
+    "initialize_2D_propeller",
     "initialize_2D_rosette",
     "initialize_2D_rings",
     "initialize_2D_polar_lissajous",

src/mrinufft/trajectories/__init__.py

@@ -4,6 +4,7 @@
     initialize_2D_spiral,
     initialize_2D_cones,
     initialize_2D_sinusoide,
+    initialize_2D_propeller,
     initialize_2D_rosette,
     initialize_2D_rings,
     initialize_2D_polar_lissajous,
@@ -32,6 +33,7 @@
     "initialize_2D_spiral",
     "initialize_2D_cones",
     "initialize_2D_sinusoide",
+    "initialize_2D_propeller",
     "initialize_2D_rosette",
     "initialize_2D_rings",
     "initialize_2D_polar_lissajous",

src/mrinufft/trajectories/trajectory2D.py

@@ -8,6 +8,7 @@
     initialize_spiral,
     compute_coprime_factors,
 )
+from .tools import rotate
 
 
 #####################
@@ -32,14 +33,14 @@ def initialize_2D_radial(Nc, Ns, tilt="uniform", in_out=False):
     # Initialize a first shot
     segment = np.linspace(-1 if (in_out) else 0, 1, Ns)
     radius = KMAX * segment
-    trajectory2D = np.zeros((Nc, Ns, 2))
-    trajectory2D[0, :, 0] = radius
+    trajectory = np.zeros((Nc, Ns, 2))
+    trajectory[0, :, 0] = radius
 
     # Rotate the first shot Nc times
     rotation = R2D(initialize_tilt(tilt, Nc) / (1 + in_out)).T
     for i in range(1, Nc):
-        trajectory2D[i] = trajectory2D[i - 1] @ rotation
-    return trajectory2D
+        trajectory[i] = trajectory[i - 1] @ rotation
+    return trajectory
 
 
 def initialize_2D_spiral(
@@ -73,15 +74,15 @@ def initialize_2D_spiral(
     angles = 2 * np.pi * nb_revolutions * (np.abs(segment) ** initialize_spiral(spiral))
 
     # Convert the first shot to Cartesian coordinates
-    trajectory2D = np.zeros((Nc, Ns, 2))
-    trajectory2D[0, :, 0] = radius * np.cos(angles)
-    trajectory2D[0, :, 1] = radius * np.sin(angles)
+    trajectory = np.zeros((Nc, Ns, 2))
+    trajectory[0, :, 0] = radius * np.cos(angles)
+    trajectory[0, :, 1] = radius * np.sin(angles)
 
     # Rotate the first shot Nc times
     rotation = R2D(initialize_tilt(tilt, Nc) / (1 + in_out)).T
     for i in range(1, Nc):
-        trajectory2D[i] = trajectory2D[i - 1] @ rotation
-    return trajectory2D
+        trajectory[i] = trajectory[i - 1] @ rotation
+    return trajectory
 
 
 def initialize_2D_cones(Nc, Ns, tilt="uniform", in_out=False, nb_zigzags=5, width=1):
@@ -112,15 +113,15 @@ def initialize_2D_cones(Nc, Ns, tilt="uniform", in_out=False, nb_zigzags=5, widt
     segment = np.linspace(-1 if (in_out) else 0, 1, Ns)
     radius = KMAX * segment
     angles = 2 * np.pi * nb_zigzags * np.abs(segment)
-    trajectory2D = np.zeros((Nc, Ns, 2))
-    trajectory2D[0, :, 0] = radius
-    trajectory2D[0, :, 1] = radius * np.sin(angles) * width * np.pi / Nc / (1 + in_out)
+    trajectory = np.zeros((Nc, Ns, 2))
+    trajectory[0, :, 0] = radius
+    trajectory[0, :, 1] = radius * np.sin(angles) * width * np.pi / Nc / (1 + in_out)
 
     # Rotate the first shot Nc times
     rotation = R2D(initialize_tilt(tilt, Nc) / (1 + in_out)).T
     for i in range(1, Nc):
-        trajectory2D[i] = trajectory2D[i - 1] @ rotation
-    return trajectory2D
+        trajectory[i] = trajectory[i - 1] @ rotation
+    return trajectory
 
 
 def initialize_2D_sinusoide(
@@ -153,15 +154,63 @@ def initialize_2D_sinusoide(
     segment = np.linspace(-1 if (in_out) else 0, 1, Ns)
     radius = KMAX * segment
     angles = 2 * np.pi * nb_zigzags * segment
-    trajectory2D = np.zeros((Nc, Ns, 2))
-    trajectory2D[0, :, 0] = radius
-    trajectory2D[0, :, 1] = KMAX * np.sin(angles) * width * np.pi / Nc / (1 + in_out)
+    trajectory = np.zeros((Nc, Ns, 2))
+    trajectory[0, :, 0] = radius
+    trajectory[0, :, 1] = KMAX * np.sin(angles) * width * np.pi / Nc / (1 + in_out)
 
     # Rotate the first shot Nc times
     rotation = R2D(initialize_tilt(tilt, Nc) / (1 + in_out)).T
     for i in range(1, Nc):
-        trajectory2D[i] = trajectory2D[i - 1] @ rotation
-    return trajectory2D
+        trajectory[i] = trajectory[i - 1] @ rotation
+    return trajectory
+
+
+def initialize_2D_propeller(Nc, Ns, nb_strips):
+    """Initialize a 2D PROPELLER trajectory, as proposed in [Pip99]_.
+
+    The PROPELLER trajectory is generally used along a specific
+    reconstruction pipeline described in [Pip99]_ to correct for
+    motion artifacts.
+
+    The acronym PROPELLER stands for Periodically Rotated
+    Overlapping ParallEL Lines with Enhanced Reconstruction,
+    and the method is also commonly known under other aliases
+    depending on the vendor, with some variations: BLADE,
+    MulitVane, RADAR, JET.
+
+    Parameters
+    ----------
+    Nc : int
+        Number of shots
+    Ns : int
+        Number of samples per shot
+    nb_strips : int
+        Number of rotated strips, must divide ``Nc``
+
+    References
+    ----------
+    .. [Pip99] Pipe, James G. "Motion correction with PROPELLER MRI:
+       application to head motion and free‐breathing cardiac imaging."
+       Magnetic Resonance in Medicine 42, no. 5 (1999): 963-969.
+    """
+    # Check for value errors
+    if Nc % nb_strips != 0:
+        raise ValueError("Nc should be divisible by nb_strips.")
+
+    # Initialize single shot
+    Nc_per_strip = Nc // nb_strips
+    trajectory = np.linspace(-1, 1, Ns).reshape((1, Ns, 1))
+
+    # Convert single shot to single strip
+    trajectory = np.tile(trajectory, reps=(Nc_per_strip, 1, 2))
+    y_axes = np.pi / 2 / nb_strips * np.linspace(-1, 1, Nc_per_strip)
+    trajectory[:, :, 1] = y_axes[:, None]
+
+    # Rotate single strip into multiple strips
+    trajectory = rotate(trajectory, nb_rotations=nb_strips, z_tilt=np.pi / nb_strips)
+    trajectory = trajectory[..., :2]  # Remove dim added by rotate
+
+    return KMAX * trajectory
 
 
 def initialize_2D_rings(Nc, Ns, nb_rings):
@@ -245,10 +294,10 @@ def initialize_2D_rosette(Nc, Ns, in_out=False, coprime_index=0):
     radius = KMAX * np.sin(Nc / (2 - odd) * angles + shift)
 
     # Convert polar to Cartesian coordinates
-    trajectory2D = np.zeros((Nc, Ns, 2))
-    trajectory2D[:, :, 0] = (radius * np.cos(angles * coprime)).reshape((Nc, Ns))
-    trajectory2D[:, :, 1] = (radius * np.sin(angles * coprime)).reshape((Nc, Ns))
-    return trajectory2D
+    trajectory = np.zeros((Nc, Ns, 2))
+    trajectory[:, :, 0] = (radius * np.cos(angles * coprime)).reshape((Nc, Ns))
+    trajectory[:, :, 1] = (radius * np.sin(angles * coprime)).reshape((Nc, Ns))
+    return trajectory
 
 
 def initialize_2D_polar_lissajous(Nc, Ns, in_out=False, nb_segments=1, coprime_index=0):
@@ -289,15 +338,15 @@ def initialize_2D_polar_lissajous(Nc, Ns, in_out=False, nb_segments=1, coprime_i
     )
 
     # Convert polar to Cartesian coordinates for one segment
-    trajectory2D = np.zeros((Nc * nb_segments, Ns, 2))
-    trajectory2D[:Nc, :, 0] = (radius * np.cos(angles)).reshape((Nc, Ns))
-    trajectory2D[:Nc, :, 1] = (radius * np.sin(angles)).reshape((Nc, Ns))
+    trajectory = np.zeros((Nc * nb_segments, Ns, 2))
+    trajectory[:Nc, :, 0] = (radius * np.cos(angles)).reshape((Nc, Ns))
+    trajectory[:Nc, :, 1] = (radius * np.sin(angles)).reshape((Nc, Ns))
 
     # Duplicate and rotate each segment
     rotation = R2D(initialize_tilt("uniform", (1 + in_out) * nb_segments))
     for i in range(Nc, Nc * nb_segments):
-        trajectory2D[i] = trajectory2D[i - Nc] @ rotation
-    return trajectory2D
+        trajectory[i] = trajectory[i - Nc] @ rotation
+    return trajectory
 
 
 #########################
@@ -327,12 +376,12 @@ def initialize_2D_lissajous(Nc, Ns, density=1):
     angles = np.pi / 2 * np.sign(segment) * np.abs(segment)
 
     # Define each shot independenty
-    trajectory2D = np.zeros((Nc, Ns, 2))
+    trajectory = np.zeros((Nc, Ns, 2))
     tilt = initialize_tilt("uniform", Nc)
     for i in range(Nc):
-        trajectory2D[i, :, 0] = KMAX * np.sin(angles)
-        trajectory2D[i, :, 1] = KMAX * np.sin(angles * density + i * tilt)
-    return trajectory2D
+        trajectory[i, :, 0] = KMAX * np.sin(angles)
+        trajectory[i, :, 1] = KMAX * np.sin(angles * density + i * tilt)
+    return trajectory
 
 
 def initialize_2D_waves(Nc, Ns, nb_zigzags=5, width=1):
@@ -361,9 +410,9 @@ def initialize_2D_waves(Nc, Ns, nb_zigzags=5, width=1):
     line = KMAX * segment
 
     # Define each shot independently
-    trajectory2D = np.zeros((Nc, Ns, 2))
+    trajectory = np.zeros((Nc, Ns, 2))
     delta = 2 * KMAX / (Nc + width)
     for i in range(Nc):
-        trajectory2D[i, :, 0] = line
-        trajectory2D[i, :, 1] = curl + delta * (i + 0.5) - (KMAX - width / Nc / 2)
-    return trajectory2D
+        trajectory[i, :, 0] = line
+        trajectory[i, :, 1] = curl + delta * (i + 0.5) - (KMAX - width / Nc / 2)
+    return trajectory