{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Suppress Gibbs oscillations\n\nMagnetic Resonance (MR) images are reconstructed from the Fourier coefficients\nof acquired k-space images. Since only a finite number of Fourier coefficients\ncan be acquired in practice, reconstructed MR images can be corrupted by Gibbs\nartefacts, which is manifested by intensity oscillations adjacent to edges of\ndifferent tissues types [1]_. Although this artefact affects MR images in\ngeneral, in the context of diffusion-weighted imaging, Gibbs oscillations\ncan be magnified in derived diffusion-based estimates [1]_, [2]_.\n\nIn the following example, we show how to suppress Gibbs artefacts of MR images.\nThis algorithm is based on an adapted version of a sub-voxel Gibbs suppression\nprocedure [3]_. Full details of the implemented algorithm can be found in\nchapter 3 of [4]_  (please cite [3]_, [4]_ if you are using this code).\n\nThe algorithm to suppress Gibbs oscillations can be imported from the denoise\nmodule of dipy:\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from dipy.denoise.gibbs import gibbs_removal"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We first apply this algorithm to a T1-weighted dataset which can be fetched\nusing the following code:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from dipy.data import get_fnames\nfrom dipy.io.image import load_nifti_data\n\n\nt1_fname, t1_denoised_fname, ap_fname = get_fnames('tissue_data')\nt1 = load_nifti_data(t1_denoised_fname)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's plot a slice of this dataset.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\n\naxial_slice = 88\nt1_slice = t1[..., axial_slice]\n\nfig = plt.figure(figsize=(15, 4))\nfig.subplots_adjust(wspace=0.2)\n\nt1_slice = np.rot90(t1_slice)\n\nplt.subplot(1, 2, 1)\nplt.imshow(t1_slice, cmap='gray', vmin=100, vmax=400)\nplt.colorbar()\nfig.savefig('structural.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: structural.png\n   :align: center\n\n   Representative slice of a T1-weighted structural image.\n\nDue to the high quality of the data, Gibbs artefacts are not visually\nevident in this dataset. Therefore, to analyse the benefits of the Gibbs\nsuppression algorithm, Gibbs artefacts are artificially introduced by\nremoving high frequencies of the image's Fourier transform.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "c = np.fft.fft2(t1_slice)\nc = np.fft.fftshift(c)\nN = c.shape[0]\nc_crop = c[64: 192, 64: 192]\nN = c_crop.shape[0]\nt1_gibbs = abs(np.fft.ifft2(c_crop)/4)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Gibbs oscillation suppression of this single data slice can be performed by\nrunning the following command:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "t1_unring = gibbs_removal(t1_gibbs, inplace=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let\u2019s plot the results:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig1, ax = plt.subplots(1, 3, figsize=(12, 6),\n                        subplot_kw={'xticks': [], 'yticks': []})\n\nax.flat[0].imshow(t1_gibbs, cmap=\"gray\", vmin=100, vmax=400)\nax.flat[0].annotate('Rings', fontsize=10, xy=(81, 70),\n                    color='red',\n                    xycoords='data', xytext=(30, 0),\n                    textcoords='offset points',\n                    arrowprops=dict(arrowstyle=\"->\",\n                                    color='red'))\nax.flat[0].annotate('Rings', fontsize=10, xy=(75, 50),\n                    color='red',\n                    xycoords='data', xytext=(30, 0),\n                    textcoords='offset points',\n                    arrowprops=dict(arrowstyle=\"->\",\n                                    color='red'))\n\nax.flat[1].imshow(t1_unring, cmap=\"gray\", vmin=100, vmax=400)\nax.flat[1].annotate('WM/GM', fontsize=10, xy=(75, 50),\n                    color='green',\n                    xycoords='data', xytext=(30, 0),\n                    textcoords='offset points',\n                    arrowprops=dict(arrowstyle=\"->\",\n                                    color='green'))\n\nax.flat[2].imshow(t1_unring - t1_gibbs, cmap=\"gray\", vmin=0, vmax=10)\nax.flat[2].annotate('Rings', fontsize=10, xy=(81, 70),\n                    color='red',\n                    xycoords='data', xytext=(30, 0),\n                    textcoords='offset points',\n                    arrowprops=dict(arrowstyle=\"->\",\n                                    color='red'))\nax.flat[2].annotate('Rings', fontsize=10, xy=(75, 50),\n                    color='red',\n                    xycoords='data', xytext=(30, 0),\n                    textcoords='offset points',\n                    arrowprops=dict(arrowstyle=\"->\",\n                                    color='red'))\nplt.show()\nfig1.savefig('Gibbs_suppression_structural.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: Gibbs_suppression_structural.png\n   :align: center\n\n   Uncorrected and corrected structural images are shown in the left\n   and middle panels, while the difference between these images is shown\n   in the right panel.\n\nThe image artificially corrupted with Gibb's artefacts is shown in the left\npanel. In this panel, the characteristic ringing profile of Gibbs artefacts\ncan be visually appreciated (see intensity oscillations pointed by the red\narrows). The corrected image is shown in the middle panel. One can appreciate\nthat artefactual oscillations are visually suppressed without compromising\nthe contrast between white and grey matter (e.g. details pointed by the green\narrow). The difference between uncorrected and corrected data is plotted in\nthe right panel which highlights the suppressed Gibbs ringing profile.\n\n\nNow let's show how to use the Gibbs suppression algorithm in diffusion-weighted\nimages. We fetch the multi-shell diffusion-weighted dataset which was kindly\nsupplied by Romain Valabr\u00e8gue, CENIR, ICM, Paris [5]_.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from dipy.data import read_cenir_multib\n\nbvals = [200, 400, 1000, 2000]\n\nimg, gtab = read_cenir_multib(bvals)\n\ndata = np.asarray(img.dataobj)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For illustration purposes, we select two slices of this dataset\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "data_slices = data[:, :, 40:42, :]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Gibbs oscillation suppression of all multi-shell data and all slices\ncan be performed in the following way:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "data_corrected = gibbs_removal(data_slices, slice_axis=2, num_processes=-1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Due to the high dimensionality of diffusion-weighted data, we recommend\nthat you specify which is the axis of data matrix that corresponds to different\nslices in the above step. This is done by using the optional parameter\n'slice_axis'.\n\nBelow we plot the results for an image acquired with b-value=0:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig2, ax = plt.subplots(1, 3, figsize=(12, 6),\n                        subplot_kw={'xticks': [], 'yticks': []})\n\nax.flat[0].imshow(data_slices[:, :, 0, 0].T, cmap='gray', origin='lower',\n                  vmin=0, vmax=10000)\nax.flat[0].set_title('Uncorrected b0')\nax.flat[1].imshow(data_corrected[:, :, 0, 0].T, cmap='gray',\n                  origin='lower', vmin=0, vmax=10000)\nax.flat[1].set_title('Corrected b0')\nax.flat[2].imshow(data_corrected[:, :, 0, 0].T - data_slices[:, :, 0, 0].T,\n                  cmap='gray', origin='lower', vmin=-500, vmax=500)\nax.flat[2].set_title('Gibbs residuals')\n\nplt.show()\nfig2.savefig('Gibbs_suppression_b0.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: Gibbs_suppression_b0.png\n   :align: center\n\n   Uncorrected (left panel) and corrected (middle panel) b-value=0 images. For\n   reference, the difference between uncorrected and corrected images is shown\n   in the right panel.\n\nThe above figure shows that the benefits of suppressing Gibbs artefacts is hard\nto observe on b-value=0 data. Therefore, diffusion derived metrics for both\nuncorrected and corrected data are computed using the mean signal diffusion\nkurtosis image technique (`example_reconst_msdki`).\n\nTo avoid unnecessary calculations on the background of the image, we also\ncompute a brain mask.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Create a brain mask\nfrom dipy.segment.mask import median_otsu\n\nmaskdata, mask = median_otsu(data_slices, vol_idx=range(10, 50),\n                             median_radius=3, numpass=1, autocrop=False,\n                             dilate=1)\n\n# Define mean signal diffusion kurtosis model\nimport dipy.reconst.msdki as msdki\n\ndki_model = msdki.MeanDiffusionKurtosisModel(gtab)\n\n# Fit the uncorrected data\ndki_fit = dki_model.fit(data_slices, mask=mask)\nMSKini = dki_fit.msk\n\n# Fit the corrected data\ndki_fit = dki_model.fit(data_corrected, mask=mask)\nMSKgib = dki_fit.msk"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's plot the results\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig3, ax = plt.subplots(1, 3, figsize=(12, 12),\n                        subplot_kw={'xticks': [], 'yticks': []})\n\nax.flat[0].imshow(MSKini[:, :, 0].T, cmap='gray', origin='lower',\n                  vmin=0, vmax=1.5)\nax.flat[0].set_title('MSK (uncorrected)')\nax.flat[0].annotate('Rings', fontsize=12, xy=(59, 63),\n                    color='red',\n                    xycoords='data', xytext=(45, 0),\n                    textcoords='offset points',\n                    arrowprops=dict(arrowstyle=\"->\",\n                                    color='red'))\n\nax.flat[1].imshow(MSKgib[:, :, 0].T, cmap='gray', origin='lower',\n                  vmin=0, vmax=1.5)\nax.flat[1].set_title('MSK (corrected)')\n\nax.flat[2].imshow(MSKgib[:, :, 0].T - MSKini[:, :, 0].T, cmap='gray',\n                  origin='lower', vmin=-0.2, vmax=0.2)\nax.flat[2].set_title('MSK (uncorrected - corrected')\nax.flat[2].annotate('Rings', fontsize=12, xy=(59, 63),\n                    color='red',\n                    xycoords='data', xytext=(45, 0),\n                    textcoords='offset points',\n                    arrowprops=dict(arrowstyle=\"->\",\n                                    color='red'))\n\nplt.show()\nfig3.savefig('Gibbs_suppression_msdki.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: Gibbs_suppression_msdki.png\n   :align: center\n\n   Uncorrected and corrected mean signal kurtosis images are shown in the left\n   and middle panels. The difference between uncorrected and corrected images\n   are show in the right panel.\n\nIn the left panel of the figure above, Gibbs artefacts can be appreciated by\nthe negative values of mean signal kurtosis (black voxels) adjacent to the\nbrain ventricle (red arrows). These negative values seem to be suppressed\nafter the `gibbs_removal` function is applied. For a better visualization of\nGibbs oscillations, the difference between corrected and uncorrected images are\nshown in the right panel.\n\n\n## References\n.. [1] Veraart, J., Fieremans, E., Jelescu, I.O., Knoll, F., Novikov, D.S.,\n       2015. Gibbs Ringing in Diffusion MRI. Magn Reson Med 76(1): 301-314.\n       https://doi.org/10.1002/mrm.25866\n.. [2] Perrone, D., Aelterman, J., Piz\u030curica, A., Jeurissen, B., Philips, W.,\n       Leemans A., 2015. The effect of Gibbs ringing artifacts on measures\n       derived from diffusion MRI. Neuroimage 120, 441-455.\n       https://doi.org/10.1016/j.neuroimage.2015.06.068.\n.. [3] Kellner, E., Dhital, B., Kiselev, V.G, Reisert, M., 2016. Gibbs\u2010ringing\n       artifact removal based on local subvoxel\u2010shifts. Magn Reson Med\n       76:1574\u20131581.\n       https://doi.org/10.1002/mrm.26054.\n.. [4] Neto Henriques, R., 2018. Advanced Methods for Diffusion MRI Data\n       Analysis and their Application to the Healthy Ageing Brain\n       (Doctoral thesis). https://doi.org/10.17863/CAM.29356\n.. [5] Valabr\u00e8gue, R. (2015). Diffusion MRI measured at multiple b-values.\n       Retrieved from:\n       https://digital.lib.washington.edu/researchworks/handle/1773/33311\n\n.. include:: ../links_names.inc\n\n"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.9.13"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}