\n\n# Need GPU warning\n\nRunning this mri-nufft example requires a GPU, and hence is NOT possible on binder currently We request you to kindly run this notebook on Google Colab by clicking the link below. Additionally, please make sure to set the runtime on Colab to use a GPU and install the below libraries before running.\n
\n "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n# Density Compensation Routines\n\nExamples of differents density compensation methods.\n\nDensity compensation depends on the sampling trajectory,and is apply before the\nadjoint operation to act as preconditioner, and should make the lipschitz constant\nof the operator roughly equal to 1.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Install libraries\n!pip install mri-nufft[gpunufft] finufft\npip install brainweb-dl # Required for data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Imports\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import brainweb_dl as bwdl\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nfrom mrinufft import get_density, get_operator\nfrom mrinufft.trajectories import initialize_2D_radial\nfrom mrinufft.trajectories.display import display_2D_trajectory"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Create sample data\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"mri_2D = np.flipud(bwdl.get_mri(4, \"T1\")[80, ...]).astype(np.float32)\n\nprint(mri_2D.shape)\n\ntraj = initialize_2D_radial(192, 192).astype(np.float32)\n\nnufft = get_operator(\"finufft\")(traj, mri_2D.shape, density=False)\nkspace = nufft.op(mri_2D)\nadjoint = nufft.adj_op(kspace)\n\nfig, axs = plt.subplots(1, 3, figsize=(15, 5))\naxs[0].imshow(abs(mri_2D))\ndisplay_2D_trajectory(traj, subfigure=axs[1])\naxs[2].imshow(abs(adjoint))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Geometry based methods\n\n## Voronoi\n\nVoronoi Parcellation attribute a weights to each k-space coordinate, inversely\nproportional to its voronoi cell area.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# .. warning::\n# The current implementation of voronoi parcellation is CPU only, and is thus\n# **very** slow in 3D ( > 1h)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"voronoi_weights = get_density(\"voronoi\", traj)\n\nnufft_voronoi = get_operator(\"finufft\")(\n traj, shape=mri_2D.shape, density=voronoi_weights\n)\nadjoint_voronoi = nufft_voronoi.adj_op(kspace)\nfig, axs = plt.subplots(1, 3, figsize=(15, 5))\naxs[0].imshow(abs(mri_2D))\naxs[0].set_title(\"Ground Truth\")\naxs[1].imshow(abs(adjoint))\naxs[1].set_title(\"no density compensation\")\naxs[2].imshow(abs(adjoint_voronoi))\naxs[2].set_title(\"Voronoi density compensation\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Cell Counting\n\nCell Counting attributes weights based on the number of trajectory point lying in a same k-space nyquist voxel.\nThis can be viewed as an approximation to the voronoi neth\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# .. note::\n# Cell counting is faster than voronoi (especially in 3D), but is less precise.\n\n# 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_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"cell_count_weights = get_density(\"cell_count\", traj, shape=mri_2D.shape, osf=2.0)\n\nnufft_cell_count = get_operator(\"finufft\")(\n traj, shape=mri_2D.shape, density=cell_count_weights, upsampfac=2.0\n)\nadjoint_cell_count = nufft_cell_count.adj_op(kspace)\nfig, axs = plt.subplots(1, 3, figsize=(15, 5))\naxs[0].imshow(abs(mri_2D))\naxs[0].set_title(\"Ground Truth\")\naxs[1].imshow(abs(adjoint))\naxs[1].set_title(\"no density compensation\")\naxs[2].imshow(abs(adjoint_cell_count))\naxs[2].set_title(\"cell_count density compensation\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Manual Density Estimation\n\nFor some analytical trajectory it is also possible to determine the density compensation vector directly.\nIn radial trajectory for instance, a sample's weight can be determined from its distance to the center.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"flat_traj = traj.reshape(-1, 2)\nweights = np.sqrt(np.sum(flat_traj**2, axis=1))\nnufft = get_operator(\"finufft\")(traj, shape=mri_2D.shape, density=weights)\nadjoint_manual = nufft.adj_op(kspace)\nfig, axs = plt.subplots(1, 3, figsize=(15, 5))\naxs[0].imshow(abs(mri_2D))\naxs[0].set_title(\"Ground Truth\")\naxs[1].imshow(abs(adjoint))\naxs[1].set_title(\"no density compensation\")\naxs[2].imshow(abs(adjoint_manual))\naxs[2].set_title(\"manual density compensation\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Operator-based method\n\n## Pipe's Method\nPipe's method is an iterative scheme, that use the interpolation and spreading kernel operator for computing the density compensation.\n\n
Warning
If this method is widely used in the literature, there exists no convergence guarantees for it.
\n\n
Note
The Pipe method is currently only implemented for gpuNUFFT.
![\"Open](\"https://colab.research.google.com/assets/colab-badge.svg\"){kind=icon}
\n \n