From f6f9e1c9c141a11e6af51fb742a01a85fa0ab2d2 Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Tue, 21 Nov 2023 22:51:50 +0100 Subject: [PATCH 1/5] Add 2D propeller trajectory --- src/mrinufft/trajectories/trajectory2D.py | 93 ++++++++++++++--------- 1 file changed, 57 insertions(+), 36 deletions(-) diff --git a/src/mrinufft/trajectories/trajectory2D.py b/src/mrinufft/trajectories/trajectory2D.py index 4fd801ce..38656f90 100644 --- a/src/mrinufft/trajectories/trajectory2D.py +++ b/src/mrinufft/trajectories/trajectory2D.py @@ -32,14 +32,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 +73,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 +112,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 +153,36 @@ 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, tilt="uniform"): + # Check for value errors + if (Nc % nb_strips != 0): + raise ValueError("Nc should be divisible by nb_strips.") + + # Initialize single shot + Nc_per_band = 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_band, 1, 2)) + y_axes = np.pi / 2 / nb_strips * np.linspace(-1, 1, Nc_per_band) + trajectory[:, :, 1] = y_axes[:, None] + + # Rotate single band 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 +266,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 +310,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 +348,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 +382,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 From d1a16530093f7ba5bd559fc42068bb41035cce64 Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Tue, 21 Nov 2023 22:57:22 +0100 Subject: [PATCH 2/5] Fix imports & format --- src/mrinufft/trajectories/trajectory2D.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/mrinufft/trajectories/trajectory2D.py b/src/mrinufft/trajectories/trajectory2D.py index 38656f90..f33b67cc 100644 --- a/src/mrinufft/trajectories/trajectory2D.py +++ b/src/mrinufft/trajectories/trajectory2D.py @@ -8,6 +8,7 @@ initialize_spiral, compute_coprime_factors, ) +from .tools import rotate ##################### @@ -166,7 +167,7 @@ def initialize_2D_sinusoide( def initialize_2D_propeller(Nc, Ns, nb_strips, tilt="uniform"): # Check for value errors - if (Nc % nb_strips != 0): + if Nc % nb_strips != 0: raise ValueError("Nc should be divisible by nb_strips.") # Initialize single shot From ba0154746d6d045501a5cbc4d2e699bf6bc39f01 Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Tue, 21 Nov 2023 23:20:46 +0100 Subject: [PATCH 3/5] Add documentation --- src/mrinufft/trajectories/trajectory2D.py | 31 ++++++++++++++++++++--- 1 file changed, 27 insertions(+), 4 deletions(-) diff --git a/src/mrinufft/trajectories/trajectory2D.py b/src/mrinufft/trajectories/trajectory2D.py index f33b67cc..f450ca0d 100644 --- a/src/mrinufft/trajectories/trajectory2D.py +++ b/src/mrinufft/trajectories/trajectory2D.py @@ -166,20 +166,43 @@ def initialize_2D_sinusoide( def initialize_2D_propeller(Nc, Ns, nb_strips, tilt="uniform"): + """Initialize a 2D PROPELLER trajectory, as proposed in []_. + + The PROPELLER trajectory is generally used along a specific + reconstruction pipeline described in []_ 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``. + tilt : str, float, optional + Tilt of the strips, by default "uniform" + """ # Check for value errors if Nc % nb_strips != 0: raise ValueError("Nc should be divisible by nb_strips.") # Initialize single shot - Nc_per_band = Nc // nb_strips + 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_band, 1, 2)) - y_axes = np.pi / 2 / nb_strips * np.linspace(-1, 1, Nc_per_band) + 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 band into multiple strips + # 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 From 21462f076a8e72f11ffbac3c177316964237c913 Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Wed, 22 Nov 2023 00:16:18 +0100 Subject: [PATCH 4/5] Add example --- examples/example_2D_trajectories.py | 51 +++++++++++++++++++++++ src/mrinufft/__init__.py | 2 + src/mrinufft/trajectories/__init__.py | 2 + src/mrinufft/trajectories/trajectory2D.py | 12 ++++-- 4 files changed, 64 insertions(+), 3 deletions(-) diff --git a/examples/example_2D_trajectories.py b/examples/example_2D_trajectories.py index 12d58010..46433808 100644 --- a/examples/example_2D_trajectories.py +++ b/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,47 @@ 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. +# - ``tilt (str, float)``: angle between each strip (in radians) +# ``(default "uniform")``. See radial +# + +trajectory = mn.initialize_2D_propeller(Nc, Ns, nb_strips=nb_repetitions, tilt=tilt) +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, tilt=tilt) +show_argument(function, arguments, one_shot=one_shot, subfigure_size=subfigure_size) + + # %% # Rings # ------- @@ -565,3 +607,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. diff --git a/src/mrinufft/__init__.py b/src/mrinufft/__init__.py index a6ed0179..69503718 100644 --- a/src/mrinufft/__init__.py +++ b/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", diff --git a/src/mrinufft/trajectories/__init__.py b/src/mrinufft/trajectories/__init__.py index f8760f91..0672d157 100644 --- a/src/mrinufft/trajectories/__init__.py +++ b/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", diff --git a/src/mrinufft/trajectories/trajectory2D.py b/src/mrinufft/trajectories/trajectory2D.py index f450ca0d..3ff4bbaa 100644 --- a/src/mrinufft/trajectories/trajectory2D.py +++ b/src/mrinufft/trajectories/trajectory2D.py @@ -166,10 +166,10 @@ def initialize_2D_sinusoide( def initialize_2D_propeller(Nc, Ns, nb_strips, tilt="uniform"): - """Initialize a 2D PROPELLER trajectory, as proposed in []_. + """Initialize a 2D PROPELLER trajectory, as proposed in [Pip99]_. The PROPELLER trajectory is generally used along a specific - reconstruction pipeline described in []_ to correct for + reconstruction pipeline described in [Pip99]_ to correct for motion artifacts. The acronym PROPELLER stands for Periodically Rotated @@ -185,9 +185,15 @@ def initialize_2D_propeller(Nc, Ns, nb_strips, tilt="uniform"): Ns : int Number of samples per shot nb_strips : int - Number of rotated strips, must divide ``Nc``. + Number of rotated strips, must divide ``Nc`` tilt : str, float, optional Tilt of the strips, by default "uniform" + + 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: From b0561b58e59f1851aee6037ebdb0ffc560f6d831 Mon Sep 17 00:00:00 2001 From: Guillaume DAVAL-FREROT Date: Wed, 22 Nov 2023 09:54:25 +0100 Subject: [PATCH 5/5] Remove unused argument tilt --- examples/example_2D_trajectories.py | 5 ++--- src/mrinufft/trajectories/trajectory2D.py | 4 +--- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/examples/example_2D_trajectories.py b/examples/example_2D_trajectories.py index 46433808..e5bfc8cc 100644 --- a/examples/example_2D_trajectories.py +++ b/examples/example_2D_trajectories.py @@ -321,11 +321,10 @@ def show_trajectory(trajectory, one_shot, figure_size): # - ``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. -# - ``tilt (str, float)``: angle between each strip (in radians) # ``(default "uniform")``. See radial # -trajectory = mn.initialize_2D_propeller(Nc, Ns, nb_strips=nb_repetitions, tilt=tilt) +trajectory = mn.initialize_2D_propeller(Nc, Ns, nb_strips=nb_repetitions) show_trajectory(trajectory, figure_size=figure_size, one_shot=one_shot) @@ -339,7 +338,7 @@ def show_trajectory(trajectory, one_shot, figure_size): # arguments = [2, 3, 4, 6] -function = lambda x: mn.initialize_2D_propeller(Nc, Ns, nb_strips=x, tilt=tilt) +function = lambda x: mn.initialize_2D_propeller(Nc, Ns, nb_strips=x) show_argument(function, arguments, one_shot=one_shot, subfigure_size=subfigure_size) diff --git a/src/mrinufft/trajectories/trajectory2D.py b/src/mrinufft/trajectories/trajectory2D.py index 3ff4bbaa..8beee398 100644 --- a/src/mrinufft/trajectories/trajectory2D.py +++ b/src/mrinufft/trajectories/trajectory2D.py @@ -165,7 +165,7 @@ def initialize_2D_sinusoide( return trajectory -def initialize_2D_propeller(Nc, Ns, nb_strips, tilt="uniform"): +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 @@ -186,8 +186,6 @@ def initialize_2D_propeller(Nc, Ns, nb_strips, tilt="uniform"): Number of samples per shot nb_strips : int Number of rotated strips, must divide ``Nc`` - tilt : str, float, optional - Tilt of the strips, by default "uniform" References ----------