Density Compensation Routines

Examples of differents density compensation methods.

Density compensation depends on the sampling trajectory,and is apply before the adjoint operation to act as preconditioner, and should make the lipschitz constant of the operator roughly equal to 1.

Imports

import brainweb_dl as bwdl
import matplotlib.pyplot as plt
import numpy as np

from mrinufft import get_density, get_operator
from mrinufft.trajectories import initialize_2D_radial
from mrinufft.trajectories.display import display_2D_trajectory

Create sample data

mri_2D = np.flipud(bwdl.get_mri(4, "T1")[80, ...]).astype(np.float32)

print(mri_2D.shape)

traj = initialize_2D_radial(192, 192).astype(np.float32)

nufft = get_operator("finufft")(traj, mri_2D.shape, density=False)
kspace = nufft.op(mri_2D)
adjoint = nufft.adj_op(kspace)

fig, axs = plt.subplots(1, 3, figsize=(15, 5))
axs[0].imshow(abs(mri_2D))
display_2D_trajectory(traj, subfigure=axs[1])
axs[2].imshow(abs(adjoint))

(256, 256)
/volatile/github-ci-mind-inria/gpu_runner/_work/_tool/Python/3.10.15/x64/lib/python3.10/site-packages/mrinufft/_utils.py:94: UserWarning: Samples will be rescaled to [-pi, pi), assuming they were in [-0.5, 0.5)
  warnings.warn(
/volatile/github-ci-mind-inria/gpu_runner/_work/_tool/Python/3.10.15/x64/lib/python3.10/site-packages/finufft/_interfaces.py:329: UserWarning: Argument `data` does not satisfy the following requirement: C. Copying array (this may reduce performance)
  warnings.warn(f"Argument `{name}` does not satisfy the following requirement: {prop}. Copying array (this may reduce performance)")

<matplotlib.image.AxesImage object at 0x741aa355c850>

As you can see, the radial sampling pattern as a strong concentration of sampling point in the center, resulting in a low-frequency biased adjoint reconstruction.

Geometry based methods

Voronoi

Voronoi Parcellation attribute a weights to each k-space coordinate, inversely proportional to its voronoi cell area.

# .. warning::
#    The current implementation of voronoi parcellation is CPU only, and is thus
#    **very** slow in 3D ( > 1h).

voronoi_weights = get_density("voronoi", traj)

nufft_voronoi = get_operator("finufft")(
    traj, shape=mri_2D.shape, density=voronoi_weights
)
adjoint_voronoi = nufft_voronoi.adj_op(kspace)
fig, axs = plt.subplots(1, 3, figsize=(15, 5))
axs[0].imshow(abs(mri_2D))
axs[0].set_title("Ground Truth")
axs[1].imshow(abs(adjoint))
axs[1].set_title("no density compensation")
axs[2].imshow(abs(adjoint_voronoi))
axs[2].set_title("Voronoi density compensation")

/volatile/github-ci-mind-inria/gpu_runner/_work/_tool/Python/3.10.15/x64/lib/python3.10/site-packages/mrinufft/_utils.py:94: UserWarning: Samples will be rescaled to [-pi, pi), assuming they were in [-0.5, 0.5)
  warnings.warn(

Text(0.5, 1.0, 'Voronoi density compensation')

Cell Counting

Cell Counting attributes weights based on the number of trajectory point lying in a same k-space nyquist voxel. This can be viewed as an approximation to the voronoi neth

# .. note::
#    Cell counting is faster than voronoi (especially in 3D), but is less precise.

# The size of the niquist voxel can be tweak by using the osf parameter. Typically as the NUFFT (and by default in MRI-NUFFT) is performed at an OSF of 2

cell_count_weights = get_density("cell_count", traj, shape=mri_2D.shape, osf=2.0)

nufft_cell_count = get_operator("finufft")(
    traj, shape=mri_2D.shape, density=cell_count_weights, upsampfac=2.0
)
adjoint_cell_count = nufft_cell_count.adj_op(kspace)
fig, axs = plt.subplots(1, 3, figsize=(15, 5))
axs[0].imshow(abs(mri_2D))
axs[0].set_title("Ground Truth")
axs[1].imshow(abs(adjoint))
axs[1].set_title("no density compensation")
axs[2].imshow(abs(adjoint_cell_count))
axs[2].set_title("cell_count density compensation")

/volatile/github-ci-mind-inria/gpu_runner/_work/_tool/Python/3.10.15/x64/lib/python3.10/site-packages/mrinufft/_utils.py:94: UserWarning: Samples will be rescaled to [-pi, pi), assuming they were in [-0.5, 0.5)
  warnings.warn(

Text(0.5, 1.0, 'cell_count density compensation')

Manual Density Estimation

For some analytical trajectory it is also possible to determine the density compensation vector directly. In radial trajectory for instance, a sample’s weight can be determined from its distance to the center.

flat_traj = traj.reshape(-1, 2)
weights = np.sqrt(np.sum(flat_traj**2, axis=1))
nufft = get_operator("finufft")(traj, shape=mri_2D.shape, density=weights)
adjoint_manual = nufft.adj_op(kspace)
fig, axs = plt.subplots(1, 3, figsize=(15, 5))
axs[0].imshow(abs(mri_2D))
axs[0].set_title("Ground Truth")
axs[1].imshow(abs(adjoint))
axs[1].set_title("no density compensation")
axs[2].imshow(abs(adjoint_manual))
axs[2].set_title("manual density compensation")

/volatile/github-ci-mind-inria/gpu_runner/_work/_tool/Python/3.10.15/x64/lib/python3.10/site-packages/mrinufft/_utils.py:94: UserWarning: Samples will be rescaled to [-pi, pi), assuming they were in [-0.5, 0.5)
  warnings.warn(

Text(0.5, 1.0, 'manual density compensation')

Operator-based method

Pipe’s Method

Pipe’s method is an iterative scheme, that use the interpolation and spreading kernel operator for computing the density compensation.

Warning

If this method is widely used in the literature, there exists no convergence guarantees for it.

Note

The Pipe method is currently only implemented for gpuNUFFT.

flat_traj = traj.reshape(-1, 2)
nufft = get_operator("gpunufft")(
    traj, shape=mri_2D.shape, density={"name": "pipe", "osf": 2}
)
adjoint_manual = nufft.adj_op(kspace)
fig, axs = plt.subplots(1, 3, figsize=(15, 5))
axs[0].imshow(abs(mri_2D))
axs[0].set_title("Ground Truth")
axs[1].imshow(abs(adjoint))
axs[1].set_title("no density compensation")
axs[2].imshow(abs(adjoint_manual))
axs[2].set_title("Pipe density compensation")
print(nufft.density)

[0.00154967 0.00707229 0.01425817 ... 0.56895894 0.46830535 0.60650784]