{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Denoise images using Adaptive Soft Coefficient Matching (ASCM)\n\nThe adaptive soft coefficient matching (ASCM) as described in [Coupe11]_ is a\nimproved extension of non-local means (NLMEANS) denoising. ASCM gives a better\ndenoised images from two standard non-local means denoised versions of the\noriginal data with different degrees sharpness. Here, one denoised input is\nmore \"smooth\" than the other (the easiest way to achieve this denoising is use\n``non_local_means`` with two different patch radii).\n\nASCM involves these basic steps\n\n* Computes wavelet decomposition of the noisy as well as denoised inputs\n\n* Combines the wavelets for the output image in a way that it takes it's\n smoothness (low frequency components) from the input with larger smoothing,\n and the sharp features (high frequency components) from the input with\n less smoothing.\n\nThis way ASCM gives us a well denoised output while preserving the sharpness\nof the image features.\n\nLet us load the necessary modules\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\nimport matplotlib.pyplot as plt\nfrom dipy.core.gradients import gradient_table\nfrom dipy.data import get_fnames\nfrom dipy.denoise.noise_estimate import estimate_sigma\nfrom dipy.io.image import load_nifti, save_nifti\nfrom dipy.io.gradients import read_bvals_bvecs\nfrom time import time\nfrom dipy.denoise.non_local_means import non_local_means\nfrom dipy.denoise.adaptive_soft_matching import adaptive_soft_matching" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Choose one of the data from the datasets in dipy_\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dwi_fname, dwi_bval_fname, dwi_bvec_fname = get_fnames('sherbrooke_3shell')\ndata, affine = load_nifti(dwi_fname)\nbvals, bvecs = read_bvals_bvecs(dwi_bval_fname, dwi_bvec_fname)\ngtab = gradient_table(bvals, bvecs)\n\nmask = data[..., 0] > 80\ndata = data[..., 1]\n\nprint(\"vol size\", data.shape)\n\nt = time()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order to generate the two pre-denoised versions of the data we will use the\n``non_local_means`` denoising. For ``non_local_means`` first we need to\nestimate the standard deviation of the noise. We use N=4 since the Sherbrooke\ndataset was acquired on a 1.5T Siemens scanner with a 4 array head coil.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "sigma = estimate_sigma(data, N=4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the denoised version of the original data which preserves sharper features,\nwe perform non-local means with smaller patch size.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "den_small = non_local_means(\n data,\n sigma=sigma,\n mask=mask,\n patch_radius=1,\n block_radius=1,\n rician=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the denoised version of the original data that implies more smoothing, we\nperform non-local means with larger patch size.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "den_large = non_local_means(\n data,\n sigma=sigma,\n mask=mask,\n patch_radius=2,\n block_radius=1,\n rician=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we perform the adaptive soft coefficient matching. Empirically we set the\nadaptive parameter in ascm to be the average of the local noise variance,\nin this case the sigma itself.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "den_final = adaptive_soft_matching(data, den_small, den_large, sigma[0])\n\nprint(\"total time\", time() - t)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To access the quality of this denoising procedure, we plot the an axial slice\nof the original data, it's denoised output and residuals.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "axial_middle = data.shape[2] // 2\n\noriginal = data[:, :, axial_middle].T\nfinal_output = den_final[:, :, axial_middle].T\ndifference = np.abs(final_output.astype(np.float64) - original.astype(np.float64))\ndifference[~mask[:, :, axial_middle].T] = 0\n\nfig, ax = plt.subplots(1, 3)\nax[0].imshow(original, cmap='gray', origin='lower')\nax[0].set_title('Original')\nax[1].imshow(final_output, cmap='gray', origin='lower')\nax[1].set_title('ASCM output')\nax[2].imshow(difference, cmap='gray', origin='lower')\nax[2].set_title('Residual')\nfor i in range(3):\n ax[i].set_axis_off()\n\nplt.savefig('denoised_ascm.png', bbox_inches='tight')\n\nprint(\"The ascm result saved in denoised_ascm.png\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: denoised_ascm.png\n :align: center\n\n Showing the axial slice without (left) and with (middle) ASCM denoising.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From the above figure we can see that the residual is really uniform in nature\nwhich dictates that ASCM denoises the data while preserving the sharpness of\nthe features.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "save_nifti('denoised_ascm.nii.gz', den_final, affine)\n\nprint(\"Saving the entire denoised output in denoised_ascm.nii.gz\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For comparison propose we also plot the outputs of the ``non_local_means``\n(both with the larger as well as with the smaller patch radius) with the ASCM\noutput.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 4)\nax[0].imshow(original, cmap='gray', origin='lower')\nax[0].set_title('Original')\nax[1].imshow(den_small[..., axial_middle].T, cmap='gray', origin='lower',\n interpolation='none')\nax[1].set_title('NLMEANS small')\nax[2].imshow(den_large[..., axial_middle].T, cmap='gray', origin='lower',\n interpolation='none')\nax[2].set_title('NLMEANS large')\nax[3].imshow(final_output, cmap='gray', origin='lower', interpolation='none')\nax[3].set_title('ASCM ')\nfor i in range(4):\n ax[i].set_axis_off()\n\nplt.savefig('ascm_comparison.png', bbox_inches='tight')\n\nprint(\"The comparison result saved in ascm_comparison.png\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: ascm_comparison.png\n :align: center\n\n Comparing outputs of the NLMEANS and ASCM.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From the above figure, we can observe that the information of two pre-denoised\nversions of the raw data, ASCM outperforms standard non-local means in\nsuppressing noise and preserving feature sharpness.\n\n## References\n\n.. [Coupe11] Pierrick Coupe, Jose Manjon, Montserrat Robles, Louis Collins.\n Adaptive Multiresolution Non-Local Means Filter for 3D MR Image Denoising.\n IET Image Processing, Institution of Engineering and Technology,\n 2011. <00645538>\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 }