{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Patch2Self: Self-Supervised Denoising via Statistical Independence\n\nPatch2Self [Fadnavis20]_ is  a self-supervised learning method for denoising\nDWI data, which uses the entire volume to learn a full-rank locally linear\ndenoiser for that volume. By taking advantage of the oversampled q-space of DWI\ndata, Patch2Self can separate structure from noise without requiring an\nexplicit model for either.\n\nClassical denoising algorithms such as Local PCA [Manjon2013]_, [Veraa2016a]_,\nNon-local Means [Coupe08]_, Total Variation Norm [Knoll11]_, etc. assume\ncertain properties on the signal structure. Patch2Self *does not* make any such\nassumption on the signal, instead using the fact that the noise across\ndifferent 3D volumes of the DWI signal originates from random fluctuations in\nthe acquired signal.\n\nSince Patch2Self only relies on the randomness of the noise, it can be applied\nat any step in the pre-processing pipeline. The design of Patch2Self is such\nthat it can work on any type of diffusion data/ any body part without\nrequiring a noise estimation or assumptions on the type of noise (such as its\ndistribution).\n\nThe Patch2Self Framework:\n\n.. figure:: https://github.com/dipy/dipy_data/blob/master/Patch2Self_Framework.PNG?raw=true\n   :scale: 60 %\n   :align: center\n\nThe above figure demonstrates the working of Patch2Self. The idea is to build\na new regressor for denoising each 3D volume of the 4D diffusion data. This is\ndone in the following 2 phases:\n\n(A) Self-supervised training: First, we extract 3D Patches from all the \u2018n\u2019\nvolumes and hold out the target volume to denoise. Each patch from the rest of\nthe \u2018(n-1)\u2019 volumes predicts the center voxel of the corresponding patch in the\ntarget volume.\n\nThis is done by using the self-supervised loss:\n$\\mathcal{L}\\left(\\Phi_{J}\right)=\\mathbb{E}\\left\\|\\Phi_{J}\\left(Y_{*, *,-j}\right)-Y_{*, 0, j}\right\\|^{2}$\n\n(B) Prediction: The same 'n-1' volumes which were used in the training are now\nfed into the regressor $\\Phi$ built in phase (A). The prediction is a\ndenoised version of held-out volume.\n\nNote: The volume to be denoised is merely used as the target in the training\nphase. But is not used in the training set for (A) nor is used to predict the\ndenoised output in (B).\n\nLet's load the necessary modules:\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nfrom dipy.data import get_fnames\nfrom dipy.io.image import load_nifti, save_nifti\nimport matplotlib.pyplot as plt\n\nfrom dipy.denoise.patch2self import patch2self"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now let's load an example dataset and denoise it with Patch2Self. Patch2Self\ndoes not require noise estimation and should work with any kind of diffusion\ndata.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "hardi_fname, hardi_bval_fname, hardi_bvec_fname = get_fnames('stanford_hardi')\ndata, affine = load_nifti(hardi_fname)\nbvals = np.loadtxt(hardi_bval_fname)\ndenoised_arr = patch2self(data, bvals, model='ols', shift_intensity=True,\n                          clip_negative_vals=False, b0_threshold=50)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The above parameters should give optimal denoising performance for Patch2Self.\nThe ordinary least squares regression (model='ols') tends to be a little slower\ndepending on the size of the data. In that case, please consider switching to\nridge regression (model='ridge').\n\nPlease do note that sometimes using ridge regression can hamper the\nperformance of Patch2Self. If so, please use model='ols'.\n\nThe array `denoised_arr` contains the denoised output obtained from Patch2Self.\n\n*Note:* Depending on the acquisition, b0 may exhibit signal attenuation or\nother artefacts that are not ideal for any denoising algorithm. We therefore\nprovide an option to skip denoising b0 volumes in the data. This can be done\nby using the option `b0_denoising=False` within Patch2Self.\n\nPlease set `shift_intensity=True` and `clip_negative_vals=False` by default to\navoid negative values in the denoised output.\n\nThe `b0_threshold` is used to separate the b0 volumes from the DWI volumes.\nChanging the value of the b0 threshold is needed if the b0 volumes in the\n`bval` file have a value greater than the default `b0_threshold`.\n\nThe default value of `b0_threshold` in DIPY is set to 50. If using data\nsuch as HCP 7T, the b0 volumes tend to have a higher b-value (>=50)\nassociated with them in the `bval` file. Please check the b-values for b0s and\nadjust the `b0_threshold` accordingly.\n\nNow let's visualize the output and the residuals obtained from the denoising.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Gets the center slice and the middle volume of the 4D diffusion data.\nsli = data.shape[2] // 2\ngra = 60  # pick out a random volume for a particular gradient direction\n\norig = data[:, :, sli, gra]\nden = denoised_arr[:, :, sli, gra]\n\n# computes the residuals\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_patch2self.png')\n\nprint(\"The result saved in denoised_patch2self.png\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: denoised_patch2self.png\n   :align: center\n\nPatch2Self preserved anatomical detail. This can be visually verified by\ninspecting the residuals obtained above. Since we do not see any structure in\nthe difference residuals, it is clear that it preserved the underlying signal\nstructure and got rid of the stochastic noise.\n\nBelow we show how the denoised data can be saved.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "save_nifti('denoised_patch2self.nii.gz', denoised_arr, affine)\n\nprint(\"Entire denoised data saved in denoised_patch2self.nii.gz\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Lastly, one can also use Patch2Self in batches if the number of gradient\ndirections is very high (>=200 volumes). For instance, if the data has 300\nvolumes, one can split the data into 2 batches, (150 directions each) and still\nget the same denoising performance. One can simply run Patch2Self using::\n\n   denoised_batch1 = patch2self(data[..., :150], bvals[:150])\n   denoised_batch2 = patch2self(data[..., 150:], bvals[150:])\n\nAfter doing this, the 2 denoised batches can be merged as follows::\n\n   denoised_p2s = np.concatenate((denoised_batch1, denoised_batch2), axis=3)\n\nOne can also consider using the above batching approach to denoise each\nshell separately if working with multi-shell data.\n\n\n## References\n\n.. [Fadnavis20] S. Fadnavis, J. Batson, E. Garyfallidis, Patch2Self:\n                Denoising Diffusion MRI with Self-supervised Learning,\n                Advances in Neural Information Processing Systems 33 (2020)\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.. [Coupe08] P. Coupe, P. Yger, S. Prima, P. Hellier, C. Kervrann, C.\n             Barillot, An Optimized Blockwise Non Local Means Denoising\n             Filter for 3D Magnetic Resonance Images, IEEE Transactions on\n             Medical Imaging, 27(4):425-441, 2008\n\n.. [Knoll11] F. Knoll, K. Bredies, T. Pock, R. Stollberger, Second order total\n             generalized variation (TGV) for MRI. Magnetic resonance in\n             medicine, 65(2), pp.480-491.\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
}