{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# SNR estimation for Diffusion-Weighted Images\n\nComputing the Signal-to-Noise-Ratio (SNR) of DW images is still an open\nquestion, as SNR depends on the white matter structure of interest as well as\nthe gradient direction corresponding to each DWI.\n\nIn classical MRI, SNR can be defined as the ratio of the mean of the signal\ndivided by the standard deviation of the underlying Gaussian noise, that is\n$SNR = mean(signal) / std(noise)$. The noise standard deviation can be computed\nfrom the background in any of the DW images. How do we compute the mean of the\nsignal, and what signal?\n\nThe strategy here is to compute a 'worst-case' SNR for DWI. Several white\nmatter structures such as the corpus callosum (CC), corticospinal tract (CST),\nor the superior longitudinal fasciculus (SLF) can be easily identified from the\ncolored-FA (CFA) map. In this example, we will use voxels from the CC, which\nhave the characteristic of being highly red in the CFA map since they are\nmainly oriented in the left-right direction. We know that the DW image closest\nto the X-direction will be the one with the most attenuated diffusion signal.\nThis is the strategy adopted in several recent papers (see [Descoteaux2011]_\nand [Jones2013]_). It gives a good indication of the quality of the DWI data.\n\nFirst, we compute the tensor model in a brain mask (see the `reconst_dti`\nexample for further explanations).\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nfrom dipy.core.gradients import gradient_table\nfrom dipy.data import get_fnames\nfrom dipy.io.gradients import read_bvals_bvecs\nfrom dipy.io.image import load_nifti, save_nifti\nfrom dipy.segment.mask import median_otsu\nfrom dipy.reconst.dti import TensorModel\n\nhardi_fname, hardi_bval_fname, hardi_bvec_fname = get_fnames('stanford_hardi')\n\ndata, affine = load_nifti(hardi_fname)\nbvals, bvecs = read_bvals_bvecs(hardi_bval_fname, hardi_bvec_fname)\ngtab = gradient_table(bvals, bvecs)\n\nprint('Computing brain mask...')\nb0_mask, mask = median_otsu(data, vol_idx=[0])\n\nprint('Computing tensors...')\ntenmodel = TensorModel(gtab)\ntensorfit = tenmodel.fit(data, mask=mask)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, we set our red-green-blue thresholds to (0.6, 1) in the x axis and\n(0, 0.1) in the y and z axes respectively. These values work well in practice\nto isolate the very RED voxels of the cfa map.\n\nThen, as assurance, we want just RED voxels in the CC (there could be noisy\nred voxels around the brain mask and we don't want those). Unless the brain\nacquisition was badly aligned, the CC is always close to the mid-sagittal\nslice.\n\nThe following lines perform these two operations and then saves the\ncomputed mask.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print('Computing worst-case/best-case SNR using the corpus callosum...')\nfrom dipy.segment.mask import segment_from_cfa\nfrom dipy.segment.mask import bounding_box\n\nthreshold = (0.6, 1, 0, 0.1, 0, 0.1)\nCC_box = np.zeros_like(data[..., 0])\n\nmins, maxs = bounding_box(mask)\nmins = np.array(mins)\nmaxs = np.array(maxs)\ndiff = (maxs - mins) // 4\nbounds_min = mins + diff\nbounds_max = maxs - diff\n\nCC_box[bounds_min[0]:bounds_max[0],\n       bounds_min[1]:bounds_max[1],\n       bounds_min[2]:bounds_max[2]] = 1\n\nmask_cc_part, cfa = segment_from_cfa(tensorfit, CC_box, threshold,\n                                     return_cfa=True)\n\nsave_nifti('cfa_CC_part.nii.gz', (cfa*255).astype(np.uint8), affine)\nsave_nifti('mask_CC_part.nii.gz', mask_cc_part.astype(np.uint8), affine)\n\nimport matplotlib.pyplot as plt\nregion = 40\nfig = plt.figure('Corpus callosum segmentation')\nplt.subplot(1, 2, 1)\nplt.title(\"Corpus callosum (CC)\")\nplt.axis('off')\nred = cfa[..., 0]\nplt.imshow(np.rot90(red[region, ...]))\n\nplt.subplot(1, 2, 2)\nplt.title(\"CC mask used for SNR computation\")\nplt.axis('off')\nplt.imshow(np.rot90(mask_cc_part[region, ...]))\nfig.savefig(\"CC_segmentation.png\", bbox_inches='tight')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: CC_segmentation.png\n   :align: center\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "x-direction, we can use all the voxels to estimate the mean signal in this\nregion.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mean_signal = np.mean(data[mask_cc_part], axis=0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "computed before and invert it to catch the outside of the brain. This could\nalso be determined manually with a ROI in the background.\n[Warning: Certain MR manufacturers mask out the outside of the brain with 0's.\nOne thus has to be careful how the noise ROI is defined].\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from scipy.ndimage import binary_dilation\nmask_noise = binary_dilation(mask, iterations=10)\nmask_noise[..., :mask_noise.shape[-1]//2] = 1\nmask_noise = ~mask_noise\n\nsave_nifti('mask_noise.nii.gz', mask_noise.astype(np.uint8), affine)\n\nnoise_std = np.std(data[mask_noise, :])\nprint('Noise standard deviation sigma= ', noise_std)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "for DW images with gradient direction that lies the closest to\nthe X, Y and Z axes.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Exclude null bvecs from the search\nidx = np.sum(gtab.bvecs, axis=-1) == 0\ngtab.bvecs[idx] = np.inf\naxis_X = np.argmin(np.sum((gtab.bvecs-np.array([1, 0, 0]))**2, axis=-1))\naxis_Y = np.argmin(np.sum((gtab.bvecs-np.array([0, 1, 0]))**2, axis=-1))\naxis_Z = np.argmin(np.sum((gtab.bvecs-np.array([0, 0, 1]))**2, axis=-1))\n\nfor direction in [0, axis_X, axis_Y, axis_Z]:\n    SNR = mean_signal[direction]/noise_std\n    if direction == 0:\n        print(\"SNR for the b=0 image is :\", SNR)\n    else:\n        print(\"SNR for direction\", direction, \" \",\n              gtab.bvecs[direction], \"is :\", SNR)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## SNR for the b=0 image is : ''42.0695455758''\n SNR for direction 58  [ 0.98875  0.1177  -0.09229] is : ''5.46995373635''\n##############################################################################\n SNR for direction 57  [-0.05039  0.99871  0.0054406] is : ''23.9329492871''\n##############################################################################\n SNR for direction 126 [-0.11825  -0.039925  0.99218 ] is : ''23.9965694823''\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Since the CC is aligned with the X axis, the lowest SNR is for that gradient\ndirection. In comparison, the DW images in the perpendicular Y and Z axes have\na high SNR. The b0 still exhibits the highest SNR, since there is no signal\nattenuation.\n\nHence, we can say the Stanford diffusion data has a 'worst-case' SNR of\napproximately 5, a 'best-case' SNR of approximately 24, and a SNR of 42 on the\nb0 image.\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### References\n\n.. [Descoteaux2011] Descoteaux, M., Deriche, R., Le Bihan, D., Mangin, J.-F.,\n   and Poupon, C. Multiple q-shell diffusion propagator imaging. Medical Image\n   Analysis, 15(4), 603, 2011.\n\n.. [Jones2013] Jones, D. K., Knosche, T. R., & Turner, R. White Matter\n   Integrity, Fiber Count, and Other Fallacies: The Dos and Don'ts of Diffusion\n   MRI. NeuroImage, 73, 239, 2013.\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
}