Source code for mrinufft.trajectories.trajectory2D

"""Functions to initialize 2D trajectories."""

import numpy as np
import numpy.linalg as nl
from scipy.interpolate import CubicSpline

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


[docs] def initialize_2D_radial(Nc, Ns, tilt="uniform", in_out=False): """Initialize a 2D radial trajectory. Parameters ---------- Nc : int Number of shots Ns : int Number of samples per shot tilt : str, float, optional Tilt of the shots, by default "uniform" in_out : bool, optional Whether to start from the center or not, by default False """ # Initialize a first shot segment = np.linspace(-1 if (in_out) else 0, 1, Ns) radius = KMAX * segment 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): trajectory[i] = trajectory[i - 1] @ rotation return trajectory
[docs] def initialize_2D_spiral( Nc, Ns, tilt="uniform", in_out=False, nb_revolutions=1, spiral="archimedes", patch_center=True, ): """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 ---------- Nc : int Number of shots Ns : int Number of samples per shot tilt : str, float, optional Tilt of the shots, by default "uniform" in_out : bool, optional Whether to start from the center or not, by default False nb_revolutions : int, optional Number of revolutions, by default 1 spiral : str, float, optional 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. """ # 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}).") # 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
[docs] 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
[docs] def initialize_2D_cones(Nc, Ns, tilt="uniform", in_out=False, nb_zigzags=5, width=1): """Initialize a 2D cone trajectory. Parameters ---------- Nc : int Number of shots Ns : int Number of samples per shot tilt : str, optional Tilt of the shots, by default "uniform" in_out : bool, optional Whether to start from the center or not, by default False nb_zigzags : float, optional Number of zigzags, by default 5 width : float, optional Width of the cone, by default 1 Returns ------- array_like 2D cone trajectory """ # Initialize a first shot segment = np.linspace(-1 if (in_out) else 0, 1, Ns) radius = KMAX * segment angles = 2 * np.pi * nb_zigzags * np.abs(segment) 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): trajectory[i] = trajectory[i - 1] @ rotation return trajectory
[docs] def initialize_2D_sinusoide( Nc, Ns, tilt="uniform", in_out=False, nb_zigzags=5, width=1 ): """Initialize a 2D sinusoide trajectory. Parameters ---------- Nc : int Number of shots Ns : int Number of samples per shot tilt : str, float, optional Tilt of the shots, by default "uniform" in_out : bool, optional Whether to start from the center or not, by default False nb_zigzags : float, optional Number of zigzags, by default 5 width : float, optional Width of the sinusoide, by default 1 Returns ------- array_like 2D sinusoide trajectory """ # Initialize a first shot segment = np.linspace(-1 if (in_out) else 0, 1, Ns) radius = KMAX * segment angles = 2 * np.pi * nb_zigzags * segment 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): trajectory[i] = trajectory[i - 1] @ rotation return trajectory
[docs] 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
[docs] def initialize_2D_rings(Nc, Ns, nb_rings): """Initialize a 2D ring trajectory, as proposed in [HHN08]_. Parameters ---------- Nc : int Number of shots Ns : int Number of samples per shot nb_rings : int Number of rings partitioning the k-space. Returns ------- array_like 2D ring trajectory References ---------- .. [HHN08] Wu, Hochong H., Jin Hyung Lee, and Dwight G. Nishimura. "MRI using a concentric rings trajectory." Magnetic Resonance in Medicine 59, no. 1 (2008): 102-112. """ if Nc < nb_rings: raise ValueError("Argument nb_rings should not be higher than Nc.") # Choose number of shots per rings nb_shots_per_rings = np.ones(nb_rings).astype(int) rings_radius = (0.5 + np.arange(nb_rings)) / nb_rings for _ in range(nb_rings, Nc): longest_shot = np.argmax(rings_radius / nb_shots_per_rings) nb_shots_per_rings[longest_shot] += 1 # Decompose each ring into shots trajectory = [] for rid in range(nb_rings): ring = np.zeros(((nb_shots_per_rings[rid]) * Ns, 2)) angles = np.linspace(0, 2 * np.pi, Ns * nb_shots_per_rings[rid]) ring[:, 0] = rings_radius[rid] * np.cos(angles) ring[:, 1] = rings_radius[rid] * np.sin(angles) for i in range(nb_shots_per_rings[rid]): trajectory.append(ring[i * Ns : (i + 1) * Ns]) return KMAX * np.array(trajectory)
[docs] def initialize_2D_rosette(Nc, Ns, in_out=False, coprime_index=0): """Initialize a 2D rosette trajectory. Parameters ---------- Nc : int Number of shots Ns : int Number of samples per shot in_out : bool, optional Whether to start from the center or not, by default False coprime_index : int, optional Index of the coprime factor, by default 0 Returns ------- array_like 2D rosette trajectory """ # Prepare to parametrize with coprime factor according to Nc parity odd = Nc % 2 coprime = compute_coprime_factors( Nc // (2 - odd), coprime_index + 1, start=1 if odd else (Nc // 2) % 2 + 1, update=2, )[-1] # Define the whole curve in polar coordinates angles = np.pi * np.linspace(-1, 1, Nc * Ns) / (1 + odd) shift = np.pi * (odd - in_out) / 2 radius = KMAX * np.sin(Nc / (2 - odd) * angles + shift) # Convert polar to Cartesian coordinates 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
[docs] def initialize_2D_polar_lissajous(Nc, Ns, in_out=False, nb_segments=1, coprime_index=0): """Initialize a 2D polar Lissajous trajectory. Parameters ---------- Nc : int Number of shots Ns : int Number of samples per shot in_out : bool, optional Whether to start from the center or not, by default False nb_segments : int, optional Number of segments, by default 1 coprime_index : int, optional Index of the coprime factor, by default 0 Returns ------- array_like 2D polar Lissajous trajectory """ # Adapt the parameters to subcases nb_segments = nb_segments * (2 - in_out) Nc = Nc // nb_segments # Define the whole curve in polar coordinates segment = np.pi / 2 * np.linspace(-1, 1, Nc * Ns) shift = np.pi * (Nc % 2 - in_out) / 2 radius = KMAX * np.sin(Nc * segment + shift) coprime_factors = compute_coprime_factors(Nc, coprime_index + 1, start=Nc % 2 + 1) angles = ( np.pi / (1 + in_out) / nb_segments * np.sin((Nc - coprime_factors[-1]) * segment) ) # Convert polar to Cartesian coordinates for one segment 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): trajectory[i] = trajectory[i - Nc] @ rotation return trajectory
######################### # NON-CIRCULAR PATTERNS # #########################
[docs] def initialize_2D_lissajous(Nc, Ns, density=1): """Initialize a 2D Lissajous trajectory. Parameters ---------- Nc : int Number of shots Ns : int Number of samples per shot density : float, optional Density of the trajectory, by default 1 Returns ------- array_like 2D Lissajous trajectory """ # Define the whole curve in Cartesian coordinates segment = np.linspace(-1, 1, Ns) angles = np.pi / 2 * np.sign(segment) * np.abs(segment) # Define each shot independenty trajectory = np.zeros((Nc, Ns, 2)) tilt = initialize_tilt("uniform", Nc) for i in range(Nc): trajectory[i, :, 0] = KMAX * np.sin(angles) trajectory[i, :, 1] = KMAX * np.sin(angles * density + i * tilt) return trajectory
[docs] def initialize_2D_waves(Nc, Ns, nb_zigzags=5, width=1): """Initialize a 2D waves trajectory. Parameters ---------- Nc : int Number of shots Ns : int Number of samples per shot nb_zigzags : float, optional Number of zigzags, by default 5 width : float, optional Width of the trajectory, by default 1 Returns ------- array_like 2D waves trajectory """ # Initialize a first shot segment = np.linspace(-1, 1, Ns) segment = np.sign(segment) * np.abs(segment) curl = KMAX * width / Nc * np.cos(nb_zigzags * np.pi * segment) line = KMAX * segment # Define each shot independently trajectory = np.zeros((Nc, Ns, 2)) delta = 2 * KMAX / (Nc + width) for i in range(Nc): trajectory[i, :, 0] = line trajectory[i, :, 1] = curl + delta * (i + 0.5) - (KMAX - width / Nc / 2) return trajectory