Add AZTEK trajectory

[DlWAT]

Addition of a function to create AZTEK trajectories

[Daval-G]

@DlWAT Thank you for this contribution
Could you edit the PR title to make it more explicit, and could you provide a more detailed PR description ?
The code style seems to be old school, is it taken from somewhere and/or directly adapted from another language ?
I will try to provide a review this week, I expect that the code can be considerably reduced using native numpy functions.

[DlWAT]

The code was originally in C, but I translated it. I had a more python-friendly version, but I was less sure of the result.
The code originates from this github. But it only gives the end of the spokes, i.e. the max gradient (here 32768). I've added an equal distribution of points on each spoke starting from the center, and put the point coordinates between [-0.5 ; 0.5].
I can't really explain the three parameters speed, shuffle and twist, as I haven't worked on them.
https://github.com/BioMaps-MRI/AZTEK
You can find more explanation in this article : AZTEK: Adaptive zero TE k-space trajectories.

[Daval-G]

Thank you for those informations. AZTEK is a trajectory we have been wanting to add for quite some time now (listed in issue #106).
Let me know to which extent you would like to contribute to its integration. For example, we usually add examples to the gallery for each new trajectory, but I can take care of that. Similarly, do you wish to make changes yourself during the review process as assignee or should we make the changes and validate them with you ? (to ensure fidelity to the paper, as from my understanding you are involved with the author team)

[Daval-G]

@DlWAT I made a few comments to question the code, but I also tried to reduce it myself without understanding much of the content. This version is reduced from ~190 lines to ~150 lines, but I am quite certain we can go below ~100 lines if I understand what is done. Note that I tested this code with hundreds of random speed/twist/shuffle combinations and both versions are identical up to floating precision (strictly identical if we avoid using np.linspace at the end). I suggest to resume from that version:

def initialize_3D_aztek(Nc: int, Ns: int, twist: float, shuffle: float, speed: int):
    # Constants
    scale = 32767 # FIXME: to challenge
    goldenRatioInv = 2 / (1 + math.sqrt(5))  # Inverse of the golden ratio

    # 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 then pick the best one
    goldenRatioAccDecim = np.zeros(N_phi_corrShift)
    for n in range(N_phi_corrShift):
        goldenRatioAcc = goldenRatioInv * (N_phi + 2 * n - N_phi_corrShift + 1)
        goldenRatioAccDecim[n] = np.abs(np.round(goldenRatioAcc) - goldenRatioAcc)
    m, bestGoldenRatioAccDecim = np.argmin(goldenRatioAccDecim), np.min(goldenRatioAccDecim)

    # Adjust N_phi and compute golden ratio correction
    N_phi += 2 * (m % N_phi_corrShift) - N_phi_corrShift + 1
    goldenRatioCorr = bestGoldenRatioAccDecim / N_phi

    # 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 = np.zeros(Nc)
    phi = np.zeros(Nc)
    nbThetaTable = np.zeros(N_phi, dtype=int)  # Number of theta angles per phi
    nbThetaTableCumul = np.zeros(N_phi + 1, dtype=int)  # Cumulative sum of theta angles
    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
    N_theta_cumul = 0

    # Compute locations over the k-spaces sphere
    max_locations = np.zeros((Nc, 3))
    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

        # Compute locations for the selected phi angle
        for i_theta_shuffle in range(nbThetaTable[i_phi_shuffle]):
            locPos = 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

            max_locations[locPos, 0] = int(scale * np.sin(theta[angPos]) * np.cos(phi[angPos]))
            max_locations[locPos, 1] = int(scale * np.sin(theta[angPos]) * np.sin(phi[angPos]))
            max_locations[locPos, 2] = int(scale * np.cos(theta[angPos]))

        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
    for n in range(speed):
        speedOrder_double += speed * goldenRatioInv
        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):
        for m in range(speedOrder[n]):
            for p in range(totSpkModSpdFac):
                if m == speedOrder[p]:
                    speedOrderCorr[n] += 1

    # Reorder locations based on speed
    order = [n // speed + speedOrder[n % speed] * totSpkOverSpdFac + speedOrderCorr[n % speed] for n in range(Nc)]
    max_locations = max_locations[np.argsort(order)]
    max_locations = max_locations / (2 * np.max(np.abs(max_locations)))  # normalize within -0.5 to 0.5

    # Make full shots from the center to the k-space border
    trajectory = np.linspace(max_locations / Ns, max_locations, Ns)
    trajectory = np.moveaxis(trajectory, 0, 1)

    return trajectory

[Daval-G]

We will need more details

[DlWAT]

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.

[Daval-G]

Thank you, that's quite helpful !

[Daval-G]

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

[DlWAT]

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.

[Daval-G]

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

[DlWAT]

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

[Daval-G]

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

[Daval-G]

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

[Daval-G]

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

[DlWAT]

Thinking back to last week, using AZTEK as a function that generates a trajectory is not the best use of it. You can think of it more as a function that distributes points on a sphere dynamically in a pseudo-random way. And so instead of using simple rays starting from the centre as I was doing, you can also use other shapes: cones etc...
And therefore use many more different ways of filling the Fourier space in radial.

[Daval-G]

Thinking back to last week, using AZTEK as a function that generates a trajectory is not the best use of it. You can think of it more as a function that distributes points on a sphere dynamically in a pseudo-random way. And so instead of using simple rays starting from the centre as I was doing, you can also use other shapes: cones etc... And therefore use many more different ways of filling the Fourier space in radial.

Fair point, that makes sense. I would suggest to start with this function, but then it is also always nice to have a reproduction of the paper (so fully radial trajectory) before adding extensions. We have something similar with the 3D cones trajectory but also with 3D Seiffert Spirals where we take a base shot (radial, cone, seiffert spiral) and then orientate it according to a location over the k-space sphere.

The dynamic aspect is new to me though, do you mean that we should not be working with a fixed amount of shots/spokes but yield them ? What do you have in mind ?

Besides that, could you have a look at my previous comments (on the code but also my question before) ? Thanks