Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add AZTEK trajectory #234

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Add AZTEK trajectory #234

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
231 changes: 231 additions & 0 deletions src/mrinufft/trajectories/trajectory3D.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,9 @@
Rz,
generate_fibonacci_circle,
)

import math

from .tools import conify, duplicate_along_axes, epify, precess, stack
from .trajectory2D import initialize_2D_radial, initialize_2D_spiral
from .utils import KMAX, Packings, Spirals, initialize_shape_norm, initialize_tilt
Expand Down Expand Up @@ -237,6 +240,234 @@ def initialize_3D_park_radial(
return trajectory


def initialize_AZTEK_trajectory(Nc: int, Ns: int, twist: float, shuffle: float, speed: int) -> NDArray:
"""
Compute the 3D radial trajectory for AZTEK imaging.

The radial shots are oriented according to an adaptive zero TE k-space trajectory,
which sustains the sequence overall 3D radial trajectory in k-space. This function
generates the amplitude values of the three gradient components Gx, Gy, Gz, and the
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here and in the function it appears that there is a misuse of the term gradients to refer to k-space locations, while gradients are closer (with unit conversion) to the derivatives of the k-space locations

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, it's a misuse of the term gradient. it simply refers to the x, y and z coordinates. Because that's what the MRI uses directly, if I'm not mistaken.

orientation sign.

Parameters
----------
Nc : int
Number of shots (spokes).
Ns : int
Number of samples per shot (spoke).
twist : float
Twist parameter for trajectory.
shuffle : float
Shuffle parameter for trajectory.
speed : int
Speed parameter for trajectory.
Comment on lines +258 to +263
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We will need more details

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I tried to explain them in an old report, but I hope the English translation doesn't take away too much of the meaning.

Twist : This parameter controls how much the spokes "twist" as they progress in k-space. In a standard trajectory, spokes follow straight radial paths. With Twist, the azimuthal angle (φ) gradually increases along the polar angle (θ), creating a spiral-like effect. A higher Twist value ensures a more even distribution of spokes, reducing gaps in k-space coverage and improving image quality.

Shuffle : This parameter changes how spokes are positioned when transitioning between arcs. Normally, arcs are placed with minimal directional change. A low Shuffle value keeps the arcs closely aligned, but a higher Shuffle applies a golden-angle shift (≈1.942 rad) to each new arc. This prevents clustering of spokes in certain areas, ensuring a more uniform distribution across k-space, which is particularly useful for motion-compensated imaging.

Speed : This parameter defines how spokes are read in a given arc. A Speed of 0 means all spokes are acquired sequentially. A higher Speed skips some spokes and returns to them in later passes, allowing faster coverage of different directions in k-space while still maintaining an even distribution. This acts as an acceleration factor, helping to balance temporal and spatial resolution.

The initial aim of this trajectory was dynamic imaging. In static, they are useless.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you, that's quite helpful !


Returns
-------
NDArray
3D matrix containing the trajectory points.

References
----------
.. [Boucneau+20] Boucneau, Tanguy, Brice Fernandez, Florent Besson,
Anne Menini, Florian Wiesinger, Emmanuel Durand, Caroline Caramella,
Luc Darrasse, and Xavier Maître.
"AZTEK: Adaptive Zero TE K-space trajectories."
Journal of Magnetic Resonance Imaging, 2020.
"""

MAX_GRADIENT = 32767
scale_factor = 1

gx, gy, gz = [0] * Nc, [0] * Nc, [0] * Nc
max_gradient = [MAX_GRADIENT]

# Constants
goldenRatioInv = 2 / (1 + math.sqrt(5)) # Inverse of the golden ratio
tmp_scale = MAX_GRADIENT / scale_factor # Temporary scaling factor

# Set maximum gradient amplitude
max_gradient[0] = MAX_GRADIENT

# Compute number of phi angles
N_phi = math.ceil(2 * math.sqrt(Nc) / (1 + twist))
N_phi_corrShift = 13 # Correction shift for phi angles

# Adjust correction shift if necessary
if N_phi < (2 * N_phi_corrShift):
N_phi_corrShift = 1
if N_phi % 2 == 1:
N_phi_corrShift += 1

# Compute golden ratio adjustments
goldenRatioAccDecim = [0.0] * (2 * N_phi_corrShift)
for n in range(N_phi_corrShift):
goldenRatioAcc = goldenRatioInv * (N_phi + 2 * n - N_phi_corrShift + 1)
goldenRatioAccDecim[n] = math.floor(goldenRatioAcc) - goldenRatioAcc
goldenRatioAccDecim[N_phi_corrShift + n] = math.ceil(goldenRatioAcc) - goldenRatioAcc

# Find the best golden ratio adjustment
m = 0
bestGoldenRatioAccDecim = goldenRatioAccDecim[m]
for n in range(1, 2 * N_phi_corrShift):
if abs(goldenRatioAccDecim[n]) < abs(goldenRatioAccDecim[m]):
m = n
bestGoldenRatioAccDecim = goldenRatioAccDecim[m]

# Adjust N_phi and compute golden ratio correction
N_phi += 2 * (m % N_phi_corrShift) - N_phi_corrShift + 1
goldenRatioCorr = bestGoldenRatioAccDecim / N_phi
Comment on lines +292 to +319
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you explain to me what is done here and how ?


# Compute mean number of theta angles per phi
mean_N_theta = Nc / N_phi
dPhi = 2 * math.pi / N_phi # Delta phi
twist_factor = (twist * N_phi) / 4 # Twist factor

# Initialize theta and phi arrays
theta = [0.0] * Nc
phi = [0.0] * Nc
nbThetaTable = [0] * N_phi # Number of theta angles per phi
nbThetaTableCumul = [0] * (N_phi + 1) # Cumulative sum of theta angles
nbThetaTableCumul[0] = 0
N_theta_cumul = 0
theta_pos_shift = 0.5 # Initial shift for theta positions

# Compute theta and phi angles
for i_phi in range(N_phi):
N_theta = math.ceil((i_phi + 1) * mean_N_theta) - math.ceil(i_phi * mean_N_theta)
for i_theta in range(N_theta):
theta_pos = (i_theta + theta_pos_shift) / N_theta
thetaTmp = math.acos(1 - 2 * theta_pos)
theta[N_theta_cumul + i_theta] = thetaTmp
phi[N_theta_cumul + i_theta] = (i_phi + twist_factor * math.cos(thetaTmp)) * dPhi

# Update theta position shift
theta_pos_shift += goldenRatioInv + goldenRatioCorr
theta_pos_shift -= math.floor(theta_pos_shift)

# Update theta table and cumulative sum
nbThetaTable[i_phi] = N_theta
nbThetaTableCumul[i_phi + 1] = nbThetaTableCumul[i_phi] + N_theta
N_theta_cumul += N_theta

# Compute shuffle step and initialize shuffle selection
shufflePhiStep = N_phi * (0.5 + shuffle * (goldenRatioInv - 0.5))
phiSelected = [0] * N_phi
i_phi_shuffle_double = -shufflePhiStep
i_phi_shuffle = 0
N_theta_cumul = 0

# Temporary arrays for gradient scales
gxSTmp = [0] * Nc
gySTmp = [0] * Nc
gzSTmp = [0] * Nc

# Compute gradient scales
for i_phi in range(N_phi):
i_phi_shuffle_double += shufflePhiStep
i_phi_shuffle = math.floor(i_phi_shuffle_double)
i_phi_shuffle %= N_phi
counter = 0
increment = 0

# Find the next unselected phi angle
while phiSelected[i_phi_shuffle] == 1:
counter += 1
increment = counter * (1 - 2 * ((counter - 1) % 2))
i_phi_shuffle += increment
if i_phi_shuffle < 0:
i_phi_shuffle += N_phi
i_phi_shuffle %= N_phi

i_phi_shuffle_double += increment
if i_phi_shuffle_double < 0:
i_phi_shuffle_double += N_phi
Comment on lines +365 to +384
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This piece can certainly be reduced, but also refactored into a function as a similar piece of code is used below for speed order

Copy link
Author

@DlWAT DlWAT Feb 11, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe something like this, i didn't check if the results was exactly the same

for i_phi in range(N_phi):
    i_phi_shuffle_double += shufflePhiStep
    i_phi_shuffle = int(i_phi_shuffle_double) % N_phi
    counter = increment = 0

    while phiSelected[i_phi_shuffle]:
        counter += 1
        increment = counter * (1 - 2 * ((counter - 1) % 2))
        i_phi_shuffle = (i_phi_shuffle + increment) % N_phi

    i_phi_shuffle_double = (i_phi_shuffle_double + increment) % N_phi


# Compute gradient scales for the selected phi angle
for i_theta_shuffle in range(nbThetaTable[i_phi_shuffle]):
gradPos = N_theta_cumul + i_theta_shuffle
if i_phi % 2 == 0:
angPos = nbThetaTableCumul[i_phi_shuffle] + i_theta_shuffle
else:
angPos = nbThetaTableCumul[i_phi_shuffle + 1] - 1 - i_theta_shuffle

gxSTmp[gradPos] = int(tmp_scale * (math.sin(theta[angPos]) * math.cos(phi[angPos])))
gySTmp[gradPos] = int(tmp_scale * (math.sin(theta[angPos]) * math.sin(phi[angPos])))
gzSTmp[gradPos] = int(tmp_scale * math.cos(theta[angPos]))
Comment on lines +394 to +396
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would like to challenge what is done here based on the fact that tmp_scale is always equal to MAX_GRADIENT=32767. This value corresponds to 2^15, so we are enforcing a 16-bit precision (15+1 sign bit). Why exactly is it done ? We are not saving computations because everything is still in 64-bit precision, it's just that information is destroyed after 16 bits.
I think we could just remove the scale and int casts


phiSelected[i_phi_shuffle] = 1
N_theta_cumul += nbThetaTable[i_phi_shuffle]

# Apply speed ordering
speed += 1
totSpkOverSpdFac = Nc // speed
totSpkModSpdFac = Nc % speed

speedOrder = [0] * speed
speedOrderSelected = [0] * speed

# Compute speed order
speedOrder_double = 0
goldenRatioSpeedStep = speed * goldenRatioInv
for n in range(speed):
speedOrder_double += goldenRatioSpeedStep
speedOrder[n] = math.floor(speedOrder_double)
speedOrder[n] %= speed
counter = 0
increment = 0

# Find the next unselected speed order
while speedOrderSelected[speedOrder[n]] == 1:
counter += 1
increment = counter * (1 - 2 * ((counter - 1) % 2))
speedOrder[n] += increment
if speedOrder[n] < 0:
speedOrder[n] += speed
speedOrder[n] %= speed

speedOrder_double += increment
if speedOrder_double < 0:
speedOrder_double += speed

speedOrderSelected[speedOrder[n]] = 1

# Apply speed correction
speedOrderCorr = [0] * speed
for n in range(speed):
if speedOrder[n] > 0:
for m in range(speedOrder[n]):
for p in range(totSpkModSpdFac):
if m == speedOrder[p]:
speedOrderCorr[n] += 1

# Reorder gradients based on speed
for n in range(Nc):
nModSpdFac = n % speed
nWithSpdFac = n // speed + speedOrder[nModSpdFac] * totSpkOverSpdFac + speedOrderCorr[nModSpdFac]

gx[nWithSpdFac] = gxSTmp[n]
gy[nWithSpdFac] = gySTmp[n]
gz[nWithSpdFac] = gzSTmp[n]

# Normalize the gradients to be within -0.5 to 0.5
max_val = max(max(abs(max(gx)), abs(min(gx))), max(abs(max(gy)), abs(min(gy))), max(abs(max(gz)), abs(min(gz))))
gx = [g / (2 * max_val) for g in gx]
gy = [g / (2 * max_val) for g in gy]
gz = [g / (2 * max_val) for g in gz]

trajectory = np.zeros((Nc, Ns, 3))

for i in range(Nc):
for j in range(Ns):
fraction = (j + 1) / Ns
trajectory[i, j, 0] = gx[i] * fraction
trajectory[i, j, 1] = gy[i] * fraction
trajectory[i, j, 2] = gz[i] * fraction
Comment on lines +460 to +465
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

fraction never starts at 0, is that something we want ? It means the spoke never starts from the center but slightly away from it


return trajectory



############################
# FREEFORM 3D TRAJECTORIES #
############################
Expand Down
Loading