{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Denoise images using the Marcenko-Pastur PCA algorithm\n\nThe PCA-based denoising algorithm exploits the redundancy across the\ndiffusion-weighted images [Manjon2013]_, [Veraart2016a]_. This algorithm has\nbeen shown to provide an optimal compromise between noise suppression and loss\nof anatomical information for different techniques such as DTI [Manjon2013]_,\nspherical deconvolution [Veraart2016a] and DKI [Henri2018]_.\n\nThe basic idea behind the PCA-based denoising algorithms is to remove the\ncomponents of the data that are classified as noise. The Principal Components\nclassification can be performed based on prior noise variance estimates\n[Manjon2013]_ (see `denoise_localpca`) or automatically based on the\nMarcenko-Pastur distribution [Veraa2016a]_. In addition to noise\nsuppression, the PCA algorithm can be used to get the standard deviation of\nthe noise [Veraa2016b]_.\n\nIn the following example, we show how to denoise diffusion MRI images and\nestimate the noise standard deviation using the PCA algorithm based\non the Marcenko-Pastur distribution [Veraa2016a]\n\nLet's load the necessary modules\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# load general modules\nimport numpy as np\nimport nibabel as nib\nimport matplotlib.pyplot as plt\nfrom time import time\n\n# load main pca function using Marcenko-Pastur distribution\nfrom dipy.denoise.localpca import mppca\n\n# load functions to fetch data for this example\nfrom dipy.data import get_fnames\n\n# load other dipy's functions that will be used for auxiliary analysis\nfrom dipy.core.gradients import gradient_table\nfrom dipy.io.image import load_nifti\nfrom dipy.io.gradients import read_bvals_bvecs\nfrom dipy.segment.mask import median_otsu\nimport dipy.reconst.dki as dki"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For this example, we use fetch to download a multi-shell dataset which was\nkindly provided by Hansen and Jespersen (more details about the data are\nprovided in their paper [Hansen2016]_). The total size of the downloaded data\nis 192 MBytes, however you only need to fetch it once.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dwi_fname, dwi_bval_fname, dwi_bvec_fname, _ = get_fnames('cfin_multib')\ndata, affine = load_nifti(dwi_fname)\nbvals, bvecs = read_bvals_bvecs(dwi_bval_fname, dwi_bvec_fname)\ngtab = gradient_table(bvals, bvecs)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For the sake of simplicity, we only select two non-zero b-values for this\nexample.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "bvals = gtab.bvals\n\nbvecs = gtab.bvecs\n\nsel_b = np.logical_or(np.logical_or(bvals == 0, bvals == 1000), bvals == 2000)\n\ndata = data[..., sel_b]\n\ngtab = gradient_table(bvals[sel_b], bvecs[sel_b])\n\nprint(data.shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As one can see from its shape, the selected data contains a total of 67\nvolumes of images corresponding to all the diffusion gradient directions\nof the selected b-values.\n\nThe PCA denoising using the Marcenko-Pastur distribution can be performed by\ncalling the function ``mppca``:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "t = time()\n\ndenoised_arr = mppca(data, patch_radius=2)\n\nprint(\"Time taken for local MP-PCA \", -t + time())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Internally, the ``mppca`` algorithm denoises the diffusion-weighted data\nusing a 3D sliding window which is defined by the variable patch_radius.\nIn total, this window should comprise a larger number of voxels than the number\nof diffusion-weighted volumes. Since our data has a total of 67 volumes, the\npatch_radius is set to 2 which corresponds to a 5x5x5 sliding window\ncomprising a total of 125 voxels.\n\nTo assess the performance of the Marcenko-Pastur PCA denoising algorithm,\nan axial slice of the original data, denoised data, and residuals are plotted\nbelow:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sli = data.shape[2] // 2\ngra = data.shape[3] - 1\norig = data[:, :, sli, gra]\nden = denoised_arr[:, :, sli, gra]\nrms_diff = np.sqrt((orig - den) ** 2)\n\nfig1, ax = plt.subplots(1, 3, figsize=(12, 6),\n                        subplot_kw={'xticks': [], 'yticks': []})\n\nfig1.subplots_adjust(hspace=0.3, wspace=0.05)\n\nax.flat[0].imshow(orig.T, cmap='gray', interpolation='none',\n                  origin='lower')\nax.flat[0].set_title('Original')\nax.flat[1].imshow(den.T, cmap='gray', interpolation='none',\n                  origin='lower')\nax.flat[1].set_title('Denoised Output')\nax.flat[2].imshow(rms_diff.T, cmap='gray', interpolation='none',\n                  origin='lower')\nax.flat[2].set_title('Residuals')\n\nfig1.savefig('denoised_mppca.png')\n\nprint(\"The result saved in denoised_mppca.png\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: denoised_mppca.png\n   :align: center\n\nThe noise suppression can be visually appreciated by comparing the original\ndata slice (left panel) to its denoised version (middle panel). The difference\nbetween original and denoised data showing only random noise indicates that\nthe data's structural information is preserved by the PCA denoising algorithm\n(right panel).\n\nBelow we show how the denoised data can be saved.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nib.save(nib.Nifti1Image(denoised_arr,\n                         affine), 'denoised_mppca.nii.gz')\n\nprint(\"Entire denoised data saved in denoised_mppca.nii.gz\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Additionally, we show how the PCA denoising algorithm affects different\ndiffusion measurements. For this, we run the diffusion kurtosis model\nbelow on both original and denoised versions of the data:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dkimodel = dki.DiffusionKurtosisModel(gtab)\n\nmaskdata, mask = median_otsu(data, vol_idx=[0, 1], median_radius=4, numpass=2,\n                             autocrop=False, dilate=1)\n\ndki_orig = dkimodel.fit(data, mask=mask)\ndki_den = dkimodel.fit(denoised_arr, mask=mask)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We use the following code to plot the MD, FA and MK estimates from the two data\nfits:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "FA_orig = dki_orig.fa\nFA_den = dki_den.fa\nMD_orig = dki_orig.md\nMD_den = dki_den.md\nMK_orig = dki_orig.mk(0, 3)\nMK_den = dki_den.mk(0, 3)\n\n\nfig2, ax = plt.subplots(2, 3, figsize=(10, 6),\n                        subplot_kw={'xticks': [], 'yticks': []})\n\nfig2.subplots_adjust(hspace=0.3, wspace=0.03)\n\nax.flat[0].imshow(MD_orig[:, :, sli].T, cmap='gray', vmin=0, vmax=2.0e-3,\n                  origin='lower')\nax.flat[0].set_title('MD (DKI)')\nax.flat[1].imshow(FA_orig[:, :, sli].T, cmap='gray', vmin=0, vmax=0.7,\n                  origin='lower')\nax.flat[1].set_title('FA (DKI)')\nax.flat[2].imshow(MK_orig[:, :, sli].T, cmap='gray', vmin=0, vmax=1.5,\n                  origin='lower')\nax.flat[2].set_title('AD (DKI)')\nax.flat[3].imshow(MD_den[:, :, sli].T, cmap='gray', vmin=0, vmax=2.0e-3,\n                  origin='lower')\nax.flat[3].set_title('MD (DKI)')\nax.flat[4].imshow(FA_den[:, :, sli].T, cmap='gray', vmin=0, vmax=0.7,\n                  origin='lower')\nax.flat[4].set_title('FA (DKI)')\nax.flat[5].imshow(MK_den[:, :, sli].T, cmap='gray', vmin=0, vmax=1.5,\n                  origin='lower')\nax.flat[5].set_title('AD (DKI)')\nplt.show()\nfig2.savefig('denoised_dki.png')\n\nprint(\"The result saved in denoised_dki.png\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: denoised_dki.png\n   :align: center\n\nIn the above figure, the DKI maps obtained from the original data are shown in\nthe upper panels, while the DKI maps from the denoised data are shown in the\nlower panels. Substantial improvements in measurement robustness can be\nvisually appreciated, particularly for the FA and MK estimates.\n\n## Noise standard deviation estimation using the Marcenko-Pastur PCA algorithm\n\nAs mentioned above, the Marcenko-Pastur PCA algorithm can also be used to\nestimate the image's noise standard deviation (std). The noise std\nautomatically computed from the ``mppca`` function can be returned by\nsetting the optional input parameter ``return_sigma`` to True.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "denoised_arr, sigma = mppca(data, patch_radius=2, return_sigma=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's plot the noise standard deviation estimate:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig3 = plt.figure('PCA Noise standard deviation estimation')\nplt.imshow(sigma[..., sli].T, cmap='gray', origin='lower')\nplt.axis('off')\nplt.show()\nfig3.savefig('pca_sigma.png')\n\nprint(\"The result saved in pca_sigma.png\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: pca_sigma.png\n   :align: center\n\nThe above figure shows that the Marcenko-Pastur PCA algorithm computes a 3D\nspatial varying noise level map. To obtain the mean noise std across all\nvoxels, you can use the following lines of code:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mean_sigma = np.mean(sigma[mask])\n\nprint(mean_sigma)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Below we use this mean noise level estimate to compute the nominal SNR of the data\n(i.e. SNR at b-value=0):\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "b0 = denoised_arr[..., 0]\n\nmean_signal = np.mean(b0[mask])\n\nsnr = mean_signal / mean_sigma\n\nprint(snr)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### References\n\n.. [Manjon2013] Manjon JV, Coupe P, Concha L, Buades A, Collins DL \"Diffusion\n                Weighted Image Denoising Using Overcomplete Local PCA\" (2013).\n                PLoS ONE 8(9): e73021. doi:10.1371/journal.pone.0073021.\n\n.. [Veraa2016a] Veraart J, Fieremans E, Novikov DS. 2016. Diffusion MRI noise\n                mapping using random matrix theory. Magnetic Resonance in\n                Medicine. doi: 10.1002/mrm.26059.\n\n.. [Henri2018] 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\n.. [Veraa2016b] Veraart J, Novikov DS, Christiaens D, Ades-aron B, Sijbers,\n                Fieremans E, 2016. Denoising of Diffusion MRI using random\n                matrix theory. Neuroimage 142:394-406.\n                doi: 10.1016/j.neuroimage.2016.08.016\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
}