{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Denoise images using Local PCA via empirical thresholds\n\nPCA-based denoising algorithms are effective denoising methods because they\nexplore the redundancy of the multi-dimensional information of\ndiffusion-weighted datasets. In this example, we will show how to\nperform the PCA-based denoising using the algorithm proposed by Manjon et al.\n[Manjon2013]_.\n\nThis algorithm involves the following steps:\n\n* First, we estimate the local noise variance at each voxel.\n\n* Then, we apply PCA in local patches around each voxel over the gradient\n  directions.\n\n* Finally, we threshold the eigenvalues based on the local estimate of sigma\n  and then do a PCA reconstruction\n\n\nTo perform PCA denoising without a prior noise standard deviation estimate\nplease see the following example instead: `denoise_mppca`\n\nLet's load the necessary modules\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport nibabel as nib\nimport matplotlib.pyplot as plt\nfrom time import time\nfrom dipy.core.gradients import gradient_table\nfrom dipy.denoise.localpca import localpca\nfrom dipy.denoise.pca_noise_estimate import pca_noise_estimate\nfrom dipy.data import get_fnames\nfrom dipy.io.image import load_nifti\nfrom dipy.io.gradients import read_bvals_bvecs"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Load one of the datasets. These data were acquired with 63 gradients and 1\nnon-diffusion (b=0) image.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dwi_fname, dwi_bval_fname, dwi_bvec_fname = get_fnames('isbi2013_2shell')\ndata, affine = load_nifti(dwi_fname)\nbvals, bvecs = read_bvals_bvecs(dwi_bval_fname, dwi_bvec_fname)\ngtab = gradient_table(bvals, bvecs)\n\nprint(\"Input Volume\", data.shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Estimate the noise standard deviation\n\nWe use the ``pca_noise_estimate`` method to estimate the value of sigma to be\nused in the local PCA algorithm proposed by Manjon et al. [Manjon2013]_.\nIt takes both data and the gradient table object as input and returns an\nestimate of local noise standard deviation as a 3D array. We return a smoothed\nversion, where a Gaussian filter with radius 3 voxels has been applied to the\nestimate of the noise before returning it.\n\nWe correct for the bias due to Rician noise, based on an equation developed by\nKoay and Basser [Koay2006]_.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "t = time()\nsigma = pca_noise_estimate(data, gtab, correct_bias=True, smooth=3)\nprint(\"Sigma estimation time\", time() - t)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Perform the localPCA using the function `localpca`\n\nThe localpca algorithm takes into account the multi-dimensional information of\nthe diffusion MR data. It performs PCA on a local 4D patch and\nthen removes the noise components by thresholding the lowest eigenvalues.\nThe eigenvalue threshold will be computed from the local variance estimate\nperformed by the ``pca_noise_estimate`` function, if this is inputted in the\nmain ``localpca`` function. The relationship between the noise variance\nestimate and the eigenvalue threshold can be adjusted using the input parameter\n``tau_factor``. According to Manjon et al. [Manjon2013]_, this parameter is set\nto 2.3.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "t = time()\n\ndenoised_arr = localpca(data, sigma, tau_factor=2.3, patch_radius=2)\n\nprint(\"Time taken for local PCA (slow)\", -t + time())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The ``localpca`` function returns the denoised data which is plotted below\n(middle panel) together with the original version of the data (left panel) and\nthe algorithm residual image (right panel) .\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sli = data.shape[2] // 2\ngra = data.shape[3] // 2\norig = data[:, :, sli, gra]\nden = denoised_arr[:, :, sli, gra]\nrms_diff = np.sqrt((orig - den) ** 2)\n\nfig, ax = plt.subplots(1, 3)\nax[0].imshow(orig, cmap='gray', origin='lower', interpolation='none')\nax[0].set_title('Original')\nax[0].set_axis_off()\nax[1].imshow(den, cmap='gray', origin='lower', interpolation='none')\nax[1].set_title('Denoised Output')\nax[1].set_axis_off()\nax[2].imshow(rms_diff, cmap='gray', origin='lower', interpolation='none')\nax[2].set_title('Residual')\nax[2].set_axis_off()\nplt.savefig('denoised_localpca.png', bbox_inches='tight')\n\nprint(\"The result saved in denoised_localpca.png\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: denoised_localpca.png\n   :align: center\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_localpca.nii.gz')\n\nprint(\"Entire denoised data saved in denoised_localpca.nii.gz\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. [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.. [Koay2006]  Koay CG, Basser PJ (2006). \"Analytically exact correction scheme\n               for signal extraction from noisy magnitude MR signals\". JMR 179:\n               317-322.\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
}