Feature Request/Discussion: analytical DCF

[fzimmermann89]

Do we want to have analytical DCFs? Where to put them? Belongs somewhere close to the trajectory calculation....

[ckolbPTB]

Analytical DCFs exist only for very few trajectories (e.g. regular radial sampling schemes) which we hardly ever use. I don't think we really need that functionality. If we add it it should be in utilities because it actually does not really depend on the trajectory (e.g. for regular radial sampling it depends on the number of points along a radial line and the angular step but not on the actual k-space locations).

[fzimmermann89]

We also can use an analytical DCF for golden angle sampling. Or any trajectory that factories: The dcf for a golden angle radial is the product dcf_spoke * dcf_angle. dcf_spoke is 1/r, dcf_angle is how a single circle is cut into pieces by the spokes. these can be calculated independently as far as I can tell...

Maybe we can either include this logic in the numerical dcf calculation or as a separate analytical_dcf calculator.

The current dcf logical is a bit conservative and only factorize in kx,ky,kz. I think it should also be possible to factorize in k1,k0

[fzimmermann89]

For a golden angle radial, @mxlutz wants to try it for his stuff and will report back.

[fzimmermann89]

@mxlutz Do you by any chance have code for a golden angle dcf?

[mxlutz]

sorry not yet. There is code for MATLAB that could be used as a reference/be translated but didnt look at it yet

[fzimmermann89]

Jipity:

import torch
import math

def fibonacci_up_to(n):
    fib = [1, 2]
    while fib[-1] < n:
        fib.append(fib[-1] + fib[-2])
    return fib

def goldcmp(k):
    k = k.squeeze()
    assert k.ndim == 2, 'k must be 2D'
    num_readouts, num_spokes = k.shape
    dk = torch.abs(k[-1, 0] - k[0, 0]) / (num_readouts - 1)

    phi = (math.sqrt(5) - 1) / 2
    fibonacci = fibonacci_up_to(num_spokes)

    idx = len(fibonacci) - 1
    while fibonacci[idx] > num_spokes:
        idx -= 1
    m = num_spokes - fibonacci[idx]

    ang_comp = torch.ones(num_spokes)

    if m == 0:
        ang_comp[:fibonacci[idx-2]] = 1 + phi
        ang_comp[fibonacci[idx-2]:fibonacci[idx-1]] = 2
        ang_comp[fibonacci[idx-1]:] = 1 + phi
    else:
        ang_comp[m:num_spokes-m] = 1 + phi
        ang_comp[:m] = 1
        ang_comp[-m:] = 1
        if m <= fibonacci[idx-3]:
            ang_comp[fibonacci[idx-2]:fibonacci[idx-1]] = 2
            ang_comp[m:fibonacci[idx-2]+m] = 1 + phi
        else:
            ang_comp[m:fibonacci[idx-1]] = 1 + phi
            ang_comp[fibonacci[idx-1]:fibonacci[idx-2]+m] = 2 * phi

    ang_comp /= ang_comp.sum()

    comp = ang_comp.repeat(num_readouts, 1)
    comp[k == 0] = 1 / (4 * num_spokes)

    comp *= torch.abs(k) / dk
    comp[k == 0] = 1 / (4 * num_spokes)

    return comp