added analytical 2D rad dcf

src/mrpro/algorithms/dcf/__init__.py

@@ -1,4 +1,5 @@
 """Density Compensation Calculation."""
 
 from mrpro.algorithms.dcf.dcf_voronoi import dcf_1d, dcf_2d3d_voronoi
-__all__ = ["dcf_1d", "dcf_2d3d_voronoi"]
+from mrpro.algorithms.dcf.dcf_radial import dcf_2dradial
+__all__ = ["dcf_1d", "dcf_2d3d_voronoi", "dcf_2dradial"]

src/mrpro/algorithms/dcf/dcf_radial.py

@@ -0,0 +1,92 @@
+"""2D radial analytical density compensation function."""
+
+import torch
+
+from mrpro.algorithms.dcf import dcf_1d
+
+
+def extract_angle_distance_along_spoke_unique(traj: torch.Tensor) -> tuple[torch.Tensor, torch.Tensor]:
+    """Extract the constant angle and unique distance along spokes using Singular Value Decomposition (SVD).
+
+    This function computes the principal spoke direction via SVD and projects each k-space
+    (kx, ky) point onto this direction to obtain its signed distance from the origin.
+    It ensures that the computed distances are identical across acquisitions and collapses
+    them into a single 1D vector.
+
+    Parameters
+    ----------
+    traj : torch.Tensor
+        k-space trajectory positions, shaped `[2, 1, acquisitions, spokes]`.
+
+    Returns
+    -------
+    angles : torch.Tensor
+        The computed angles for each acquisition, shaped `[acquisitions]`.
+
+    distances_along_spoke : torch.Tensor
+        The unique signed distances along the spoke, shaped `[spokes]`.
+
+    Raises
+    ------
+    ValueError
+        If the computed distances are not unique across acquisitions.
+    """
+    m = traj.squeeze(1).permute(1, 2, 0)
+    _, _, v = torch.svd(m)
+    principal_directions = v[:, :, 0]
+
+    angles = torch.atan2(principal_directions[:, 1], principal_directions[:, 0])
+    angles %= torch.pi
+
+    distances_along_spoke = torch.einsum('aji,ai->aj', m, principal_directions)
+
+    # Ensure all acquisitions have the same distances
+    if not torch.allclose(distances_along_spoke[0], distances_along_spoke):
+        raise ValueError('Distances along spokes are not unique across acquisitions!')
+
+    distances_along_spoke_unique = distances_along_spoke[0]
+
+    return angles, distances_along_spoke_unique
+
+
+def dcf_2dradial(traj: torch.Tensor) -> torch.Tensor:
+    """Calculate sample density compensation function for an 2D radial trajectory.
+
+    The density compensation function is calculated by calculating the fraction of the area of the k-space.
+    The dcf is calculated along the polar angle and the spoke dim and is then combined.
+    The calculation is based on the formula from
+    https://users.fmrib.ox.ac.uk/~karla/reading_group/lecture_notes/AdvRecon_Pauly_read.pdf for equidistant radial
+    sampling.
+
+    Parameters
+    ----------
+    traj
+        k-space positions `(2, 1, k1, k0)`
+
+    Returns
+    -------
+        density compensation values for analytical radial trajectory `(1, 1, k1, k0)`
+    """
+    angles, distances_along_spoke_unique = extract_angle_distance_along_spoke_unique(traj)
+
+    # get dcf along the polar angle dim which should sum to 1 to get the fraction of the areas
+    dcf_polar = dcf_1d(angles, periodicity=torch.pi)
+    dcf_polar /= dcf_polar.sum()
+
+    # get dcf along the spoke dim, this is mainly pi*r^2 of the outer ring - pi*r^2 of the inner ring
+    dcf_spoke = torch.pi * dcf_1d(distances_along_spoke_unique**2)
+
+    dcf = torch.outer(dcf_polar, dcf_spoke)
+
+    # fix center value
+    zero_indices = torch.nonzero(distances_along_spoke_unique == 0)
+
+    # Assert that there is only one zero in the array
+    if len(zero_indices) != 1:
+        raise ValueError('The array should contain exactly one zero.')
+
+    center_idx = int(zero_indices[0].item())
+
+    dcf[:, center_idx] = dcf_spoke[center_idx] / 4 / len(angles)
+
+    return dcf.unsqueeze(0).unsqueeze(0)

src/mrpro/algorithms/dcf/dcf_voronoi.py

@@ -15,7 +15,7 @@ def _volume(v: ArrayLike):
     return ConvexHull(v).volume
 
 
-def dcf_1d(traj: torch.Tensor) -> torch.Tensor:
+def dcf_1d(traj: torch.Tensor, periodicity: float | None = None) -> torch.Tensor:
     """Calculate sample density compensation function for 1D trajectory.
 
     This function operates on a single `other` sample.
@@ -25,6 +25,8 @@ def dcf_1d(traj: torch.Tensor) -> torch.Tensor:
     ----------
     traj
         k-space positions, 1D tensor
+    periodicity
+        periodicity of the trajectory, if None, the trajectory is assumed to be non-periodic
 
     Returns
     -------
@@ -48,8 +50,12 @@ def dcf_1d(traj: torch.Tensor) -> torch.Tensor:
 
     if (elements := len(traj_sorted)) >= 3:
         central_diff = torch.nn.functional.conv1d(traj_sorted[None, None, :], kernel)[0, 0]
-        first = traj_sorted[1] - traj_sorted[0]
-        last = traj_sorted[-1] - traj_sorted[-2]
+        if periodicity:
+            first = traj_sorted[1] / 2 - (traj_sorted[-1] - periodicity) / 2
+            last = traj_sorted[0] / 2 - (traj_sorted[-2] - periodicity) / 2
+        else:
+            first = traj_sorted[1] - traj_sorted[0]
+            last = traj_sorted[-1] - traj_sorted[-2]
         central_diff = torch.cat((first[None], central_diff, last[None]), -1)
     elif elements == 2:
         diff = traj_sorted[1] - traj_sorted[0]

tests/algorithms/dcf/test_dcf_radial.py

@@ -0,0 +1,60 @@
+"""Tests for algorithms to calculate the DCF with voronoi."""
+
+import pytest
+import torch
+from einops import repeat
+from mrpro.algorithms.dcf import dcf_2dradial
+from mrpro.data import KTrajectory
+
+
+def example_traj_rad_2d(n_kr, n_ka, phi0=0.0, broadcast=True):
+    """Create 2D radial trajectory with uniform angular gap."""
+    krad = repeat(
+        torch.linspace(-n_kr // 2, n_kr // 2 - 1, n_kr) / n_kr,
+        'k0 -> other coils k2 k1 k0',
+        other=1,
+        coils=1,
+        k2=1,
+        k1=1,
+    )
+    kang = repeat(
+        torch.linspace(0, n_ka - 1, n_ka) * (torch.pi / n_ka) + phi0,
+        'k1 -> other coils k2 k1 k0',
+        other=1,
+        coils=1,
+        k2=1,
+        k0=1,
+    )
+    kz = torch.zeros(1, 1, 1, 1, 1)
+    ky = torch.sin(kang) * krad
+    kx = torch.cos(kang) * krad
+    trajectory = KTrajectory(kz, ky, kx, repeat_detection_tolerance=1e-8 if broadcast else None)
+    return trajectory
+
+
+@pytest.mark.parametrize(
+    ('n_kr', 'n_ka', 'phi0', 'broadcast'),
+    [
+        (20, 20, 0, True),
+        (20, 2, 0, True),
+        (20, 20, torch.pi / 4, True),
+        (20, 2, torch.pi / 4, True),
+        (20, 2, 0, False),
+    ],
+)
+def test_dcf_2drad_analytical_equidist(n_kr, n_ka, phi0, broadcast):
+    """Compare 2d dcf calculation for 2D radial trajectory to
+    analytical solution for 2D equidistant dcf."""
+    # 2D radial trajectory
+    traj = example_traj_rad_2d(n_kr, n_ka, phi0, broadcast)
+    trajectory = traj.as_tensor()
+
+    dcf_2drad = dcf_2dradial(trajectory[1:3, 0, 0, ...])
+    krad_idx = torch.linspace(-n_kr // 2, n_kr // 2 - 1, n_kr)
+    dcf_analytical_equidist = torch.pi / n_ka * torch.abs(krad_idx) * (1 / n_kr) ** 2
+    dcf_analytical_equidist[krad_idx == 0] = 2 * torch.pi / n_ka * 1 / 8 * (1 / n_kr) ** 2
+    dcf_analytical_equidist = torch.repeat_interleave(
+        repeat(dcf_analytical_equidist, 'k0->k2 k1 k0', k1=1, k2=1), n_ka, dim=-2
+    ).unsqueeze(0)
+    # Do not test outer points because they have to be approximated and cannot be calculated exactly
+    torch.testing.assert_close(dcf_analytical_equidist[:, :, :, 1:-1], dcf_2drad[:, :, :, 1:-1])