{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Noise estimation using PIESNO\n\nOften, one is interested in estimating the noise in the diffusion signal. One\nof the methods to do this is the Probabilistic Identification and Estimation of\nNoise (PIESNO) framework [Koay2009]_. Using this method, one can detect the\nstandard deviation of the noise from Diffusion-Weighted Imaging (DWI). PIESNO\nalso works with multiple channel DWI datasets that are acquired from N array\ncoils for both SENSE and GRAPPA reconstructions.\n\nThe PIESNO method works in two steps:\n\n1) First, it finds voxels that are most likely background voxels. Intuitively,\nthese voxels have very similar diffusion-weighted intensities (up to some\nnoise) in the fourth dimension of the DWI dataset. White matter, gray matter\nor CSF voxels have diffusion intensities that vary quite a lot across different\ndirections.\n\n2) From these estimated background voxels and the input number of coils $N$,\nPIESNO finds what sigma each Gaussian from each of the $N$ coils would have\ngenerated the observed Rician ($N = 1$) or non-central Chi ($N > 1$)\ndistributed noise profile in the DWI datasets.\n\nPIESNO makes an important assumption: the Gaussian noise standard deviation is\nassumed to be uniform. The noise is uniform across multiple slice locations or\nacross multiple images of the same location.\n\nFor the full details, please refer to the original paper.\n\nIn this example, we will demonstrate the use of PIESNO with a 3-shell data-set.\n\nWe start by importing necessary modules and functions and loading the data:\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nfrom dipy.denoise.noise_estimate import piesno\nfrom dipy.data import get_fnames\nfrom dipy.io.image import load_nifti, save_nifti\n\ndwi_fname, dwi_bval_fname, dwi_bvec_fname = get_fnames('sherbrooke_3shell')\ndata, affine = load_nifti(dwi_fname)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now that we have fetched a dataset, we must call PIESNO with the right number\nof coils used to acquire this dataset. It is also important to know what\nwas the parallel reconstruction algorithm used. Here, the data comes from a\nGRAPPA reconstruction, was acquired with a 12-elements head coil available on\nthe Tim Trio Siemens, for which the 12 coil elements are combined into 4 groups\nof 3 coil elements each. The signal is therefore received through 4 distinct\ngroups of receiver channels, yielding N = 4. Had we used a GE acquisition, we\nwould have used N=1 even if multiple channel coils are used because GE uses a\nSENSE reconstruction, which has a Rician noise nature and thus N is always 1.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sigma, mask = piesno(data, N=4, return_mask=True)\n\naxial = data[:, :, data.shape[2] // 2, 0].T\naxial_piesno = mask[:, :, data.shape[2] // 2].T\n\nimport matplotlib.pyplot as plt\n\nfig, ax = plt.subplots(1, 2)\nax[0].imshow(axial, cmap='gray', origin='lower')\nax[0].set_title('Axial slice of the b=0 data')\nax[1].imshow(axial_piesno, cmap='gray', origin='lower')\nax[1].set_title('Background voxels from the data')\nfor a in ax:\n    a.set_axis_off()\n\nplt.savefig('piesno.png', bbox_inches='tight')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: piesno.png\n   :align: center\n\n   Showing the mid axial slice of the b=0 image (left) and estimated\n   background voxels (right) used to estimate the noise standard deviation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "save_nifti('mask_piesno.nii.gz', mask.astype(np.uint8), affine)\n\nprint('The noise standard deviation is sigma = ', sigma)\nprint('The std of the background is =',\n      np.std(data[mask[..., :].astype(bool)]))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here, we obtained a noise standard deviation of 7.26. For comparison, a simple\nstandard deviation of all voxels in the estimated mask (as done in the previous\nexample `example_snr_in_cc`) gives a value of 6.1.\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## References\n\n.. [Koay2009] Koay C.G., E. Ozarslan, C. Pierpaoli. Probabilistic\n   Identification and Estimation of Noise (PIESNO): A self-consistent approach\n   and its applications in MRI. JMR, 199(1):94-103, 2009.\n\n.. include:: ../links_names.inc\n\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
}