PTB-MR / mrpro

MR image reconstruction and processing.
https://ptb-mr.github.io/mrpro/
Apache License 2.0
17 stars 2 forks source link

Dcf calculation fails for single spiral arm #84

Closed ckolbPTB closed 1 year ago

ckolbPTB commented 1 year ago

When the trajectory is a single spiral arm, then DcfData.from_traj_voronoi detects this as a single non-singleton dimension and carries out a 1D calculation. This fails for a spiral trajectory because it covers a 2D k-space.

Possible solution: Check if points lie on a straight line, otherwise use voronoi This will still fail for certain trajectories but this then needs to be handled by the user with a case by case solution.

fzimmermann89 commented 1 year ago

This could be the "rotated direction" case I mentioned in the original PR and would be a major problem with the current logic, let me think a bit about a good solution

fzimmermann89 commented 1 year ago

Also, thinking about it a bit more: I am a bit confused about the correct result. If I just have a single spiral arm, The correct DCF should be constant. As this is no different than a cartesian, single line trajectory. The density is constant over the radius.

But our "analytical result" is not. I believe this is wrong.

ckolbPTB commented 1 year ago

If I just have a single spiral arm, The correct DCF should be constant.

No, I would not say so. I mean there is no "correct" DCF in the first place but the best option is probably to make sure that the energy of the spiral k-space data is distributed uniformely in the 2D k-space the spiral covers. This means that the further apart the k-space points are the higher the DCF weighting.

ckolbPTB commented 1 year ago

I think the current implementation only works if the singleton dimension is aligned with k0, k1 or k2. If it is not, then the singleton dimension covers a 2D or 3D k-space and the DCF cannot be calculated by simply doing a 1D difference calculation.

fzimmermann89 commented 1 year ago

Ah, sorry! I misunderstood. I was thinking about a radial trajectory with a single spoke.

Can you add a test case for a spiral trajectory? Or even just post some code to create the trajectory manually?

ckolbPTB commented 1 year ago

Also a radial trajectory with a single spoke is an issue but one which I don't think we can solve (for alpha = 0° or 90° kx/ky are aligned with k0/k1 but the DCF should still be calculated based on a 2D k-space). Nevertheless, this is a very special case and needs to handled separately by the user.

I will add a spiral test case.

fzimmermann89 commented 1 year ago

Ah, sorry! I misunderstood. I was thinking about a radial trajectory with a single spoke.

Can you add a test case for a spiral trajectory? Or even just post some code to create the trajectory manually?

So, the issue is, that our dcf-code thinks there is broadcasting going on. But actually, it is just a single spiral. If there were two spirals, the index would not be 1 anymore. This should be easy to fix. Either by adding some logic to the trajectory or the dcf. Let me check which way looks better once I get a test case.

I think the current implementation only works if the singleton dimension is aligned with k0, k1 or k2. If it is not, then the singleton dimension covers a 2D or 3D k-space and the DCF cannot be calculated by simply doing a 1D difference calculation.

Almost. The current implementation only detects repeats if the repeat direction is axis-aligned, yes. But this fails the other way around: If the repeat-direction is non-axis aligned, we will have to do more vornoi than necessary. But this is an edge-case where the performance degradation should not matter to much.

ckolbPTB commented 1 year ago

But this fails the other way around: If the repeat-direction is non-axis aligned

Hmmm...I can't think of a case where this could happen...

ckolbPTB commented 1 year ago

This would be the code for a test spiral:

def example_traj_spiral_2d(nkr, nki, nka):
    """ Create 2D spiral trajectory with nkr points along each spiral arm,
    nki turns per spiral arm and nka spiral arms."""
    ang = torch.linspace(0, 2*torch.pi*nki, nkr)
    start_ang = torch.linspace(0, 2*torch.pi*(1-1/nka), nka)
    kz = torch.zeros(1, 1, 1, 1)
    kx = (ang[None,:] * torch.cos(ang[None,:] + start_ang[:,None]))[None, None, :, :]
    ky = (ang[None,:] * torch.sin(ang[None,:] + start_ang[:,None]))[None, None, :, :]
    ktraj = KTrajectory(kz, ky, kx)
    return ktraj

Not sure about a test case though. Afaik there is no anlytical solution for this...

fzimmermann89 commented 1 year ago

... Interpolating the DCF from the trajectory to a regular grid should result in a grid with almost constant values?

EDIT: "Calculating the DCF with the spirals at exactly the same positions should result in half the dcf for each spiral" catches the bug