WIP: New PIESNO example and small corrections

dipy/denoise/noise_estimate.py

@@ -3,6 +3,7 @@
 import numpy as np
 
 from scipy.special import gammainccinv
+from scipy.stats import mode
 
 
 def _inv_nchi_cdf(N, K, alpha):
@@ -11,7 +12,7 @@ def _inv_nchi_cdf(N, K, alpha):
     return gammainccinv(N * K, 1 - alpha) / K
 
 
-def piesno(data, N=1, alpha=0.01, l=100, itermax=100, eps=1e-5, return_mask=False):
+def piesno(data, N, alpha=0.01, l=100, itermax=100, eps=1e-5, return_mask=False):
     """
     Probabilistic Identification and Estimation of Noise (PIESNO)
     A routine for finding the underlying gaussian distribution standard
@@ -28,6 +29,11 @@ def piesno(data, N=1, alpha=0.01, l=100, itermax=100, eps=1e-5, return_mask=Fals
 
     N : int
         The number of phase array coils of the MRI scanner.
+        If your scanner does a SENSE reconstruction, ALWAYS use N=1, as the noise
+        profile is always Rician. 
+        If your scanner does a GRAPPA reconstruction, set N as the number
+        of phase array coils.
+
 
     alpha : float
         Probabilistic estimation threshold for the gamma function.
@@ -42,7 +48,7 @@ def piesno(data, N=1, alpha=0.01, l=100, itermax=100, eps=1e-5, return_mask=Fals
     eps : float
         Tolerance for the convergence criterion. Convergence is
         reached if two subsequent estimates are smaller than eps.
-    
+
     return_mask : bool
         If True, return a mask identyfing all the pure noise voxel
         that were found.
@@ -75,6 +81,103 @@ def piesno(data, N=1, alpha=0.01, l=100, itermax=100, eps=1e-5, return_mask=Fals
     Journal of Magnetic Resonance 2009; 197: 108-119.
     """
 
+    # This method works on a 2D array with repetitions as the third dimension,
+    # so process the dataset slice by slice.
+
+    if data.ndim < 3:
+        raise ValueError("This function only works on datasets of at least 3 dimensions.")
+
+    if data.ndim == 4:
+
+        sigma = np.zeros(data.shape[-2], dtype=np.float32)
+        mask_noise = np.zeros(data.shape[:-1], dtype=np.bool)
+
+        for idx in range(data.shape[-2]):
+            sigma[idx], mask_noise[..., idx] = _piesno_3D(data[..., idx, :], N, 
+                                                          alpha=alpha, l=l, itermax=itermax, eps=eps)
+
+        #debug    
+        # Take the mode of all the sigmas from each slice as the best estimate,
+        # this should be stable with more or less 50% of the guesses at the same value.
+        #print(sigma)
+
+        sigma, num = mode(sigma, axis=None)
+        #debug
+        #print(sigma, num)
+
+    else:
+        sigma, mask_noise = _piesno_3D(data, N, alpha=alpha, l=l, itermax=itermax, eps=eps)
+
+    if return_mask:
+        return sigma, mask_noise
+
+    return sigma
+
+
+def _piesno_3D(data, N, alpha=0.01, l=100, itermax=100, eps=1e-5):
+    """
+    Probabilistic Identification and Estimation of Noise (PIESNO)
+    This is the slice by slice version.
+
+    Parameters
+    -----------
+    data : ndarray
+        The magnitude signals to analyse. The last dimension must contain the
+        same realisation of the volume, such as dMRI or fMRI data.
+
+    N : int
+        The number of phase array coils of the MRI scanner.
+        If your scanner does a SENSE reconstruction, ALWAYS use N=1, as the noise
+        profile is always Rician. 
+        If your scanner does a GRAPPA reconstruction, set N as the number
+        of phase array coils.
+
+    alpha : float
+        Probabilistic estimation threshold for the gamma function.
+
+    l : int
+        number of initial estimates for sigma to try.
+
+    itermax : int
+        Maximum number of iterations to execute if convergence
+        is not reached.
+
+    eps : float
+        Tolerance for the convergence criterion. Convergence is
+        reached if two subsequent estimates are smaller than eps.
+
+    return_mask : bool
+        If True, return a mask identyfing all the pure noise voxel
+        that were found.
+
+    Returns
+    --------
+    sigma : float
+        The estimated standard deviation of the gaussian noise.
+
+    mask : ndarray
+        A boolean mask indicating the voxels identified as pure noise.
+
+    Note
+    ------
+    This function assumes two things : 1. The data has a noisy, non-masked
+    background and 2. The data is a repetition of the same measurements
+    along the last axis, i.e. dMRI or fMRI data, not structural data like T1/T2.
+
+    References
+    ------------
+
+    .. [1] Koay CG, Ozarslan E and Pierpaoli C.
+    "Probabilistic Identification and Estimation of Noise (PIESNO):
+    A self-consistent approach and its applications in MRI."
+    Journal of Magnetic Resonance 2009; 199: 94-103.
+
+    .. [2] Koay CG, Ozarslan E and Basser PJ.
+    "A signal transformational framework for breaking the noise floor
+    and its applications in MRI."
+    Journal of Magnetic Resonance 2009; 197: 108-119.
+    """
+
     # Get optimal quantile for N if available, else use the median.
     opt_quantile = {1: 0.79681213002002,
                     2: 0.7306303027491917,
@@ -100,7 +203,7 @@ def piesno(data, N=1, alpha=0.01, l=100, itermax=100, eps=1e-5, return_mask=Fals
     K = data.shape[-1]
     sum_m2 = np.sum(data**2, axis=-1)
 
-    sigma = np.zeros_like(phi)
+    sigma = np.zeros(phi.shape, dtype=phi.dtype)
     mask = np.zeros(phi.shape + data.shape[:-1])
 
     lambda_minus = _inv_nchi_cdf(N, K, alpha/2)
@@ -137,13 +240,9 @@ def piesno(data, N=1, alpha=0.01, l=100, itermax=100, eps=1e-5, return_mask=Fals
         # Remember the biggest omega array as giving the optimal
         # sigma amongst all initial estimates from phi
         if omega_size > max_length_omega:
-            pos, max_length_omega, best_mask = num, omega_size, idx
+            pos, max_length_omega = num, omega_size
 
         sigma[num] = sig
         mask[num] = idx
 
-    if return_mask:
-        return sigma[pos], mask[pos]
-
-    return sigma[pos]
-
+    return sigma[pos], mask[pos]

doc/examples/denoise_nlmeans.py

@@ -20,7 +20,7 @@
 img, gtab = read_sherbrooke_3shell()
 
 data = img.get_data()
-aff = img.get_affine()
+affine = img.get_affine()
 
 mask = data[..., 0] > 80
 
@@ -71,7 +71,7 @@
    **Showing the middle axial slice without (left) and with (right) NLMEANS denoising**.
 """
 
-nib.save(nib.Nifti1Image(den, aff), 'denoised.nii.gz')
+nib.save(nib.Nifti1Image(den, affine), 'denoised.nii.gz')
 
 """
 

doc/examples/piesno.py

@@ -0,0 +1,100 @@
+"""
+=============================
+Noise estimation using PIESNO
+=============================
+
+Using PIESNO [Koay2009]_ one can detect the standard deviation of the noise from
+diffusion-weighted imaging (DWI). PIESNO also works with multiple channel
+DWI datasets that are acquired from N array coils for both SENSE and
+GRAPPA reconstructions. 
+
+The PIESNO paper [Koay2009]_ is mathematically quite involved. 
+PIESNO works in two steps. 1) First, it finds voxels that are most likely
+background voxels. Intuitively, these voxels have very similar 
+diffusion-weighted intensities (up to some noise) in the fourth dimension
+of the DWI dataset. White matter, gray matter or CSF voxels have diffusion intensities 
+that vary quite a lot across different directions. 
+2) From these esimated background voxels and 
+the input number of coils N, PIESNO finds what sigma each Gaussian 
+from each of the N coils would have generated the observed Rician (N=1) 
+or non-central Chi (N>1) distributed noise profile in the DWI datasets.
+[Koay2009]_ gives all the details.
+
+PIESNO makes an important assumption: the
+Gaussian noise standard deviation is assumed to be uniform. The noise
+is uniform across multiple slice locations or across multiple images of the same location.  
+
+"""
+
+import nibabel as nib
+import numpy as np
+from dipy.denoise.noise_estimate import piesno
+from dipy.data import fetch_sherbrooke_3shell, read_sherbrooke_3shell
+
+
+fetch_sherbrooke_3shell()
+img, gtab = read_sherbrooke_3shell()
+data = img.get_data()
+
+"""
+
+Now that we have fetched a dataset, we must call PIESNO with right number
+of coil used to acquire this dataset. It is also important to know what 
+was the parallel reconstruction algorithm used. 
+Here, the data comes from a GRAPPA reconstruction from
+a 12-elements head coil available on the Tim Trio Siemens, for which
+the 12 coil elements are combined into 4 groups of 3 coil elements
+each. The signal is received through 4 distinct groups of receiver channels,
+yielding N = 4. Had we used a GE acquisition, we would have used N=1 even
+if multiple channel coils are used because GE uses a SENSE reconstruction, 
+which has a Rician noise nature and thus N is always 1.
+
+"""
+
+sigma, mask = piesno(data, N=4, return_mask=True)
+
+axial = data[:, :, data.shape[2] / 2, 0].T
+axial_piesno = mask[:, :, data.shape[2] / 2].T
+
+import matplotlib.pyplot as plt
+fig, ax = plt.subplots(1, 2)
+ax[0].imshow(axial, cmap='gray', origin='lower')
+ax[0].set_title('Axial slice of the b=0 data')
+ax[1].imshow(axial_piesno, cmap='gray', origin='lower')
+ax[1].set_title('Background voxels from the data')
+for a in ax:
+    a.set_axis_off()
+
+# Uncomment the coming line if you want a window to show the images. 
+#plt.show()
+plt.savefig('piesno.png', bbox_inches='tight')
+
+"""
+.. figure:: piesno.png
+   :align: center
+
+   **Showing the mid axial slice of the b=0 image (left) and estimated background voxels (right) used to estimate the noise standard deviation**.
+"""
+
+nib.save(nib.Nifti1Image(mask, img.get_affine(), img.get_header()), 
+         'mask_piesno.nii.gz')
+
+print('The noise standard deviation is sigma= ', sigma)
+print('The std of the background is =', np.std(data[mask[...,None].astype(np.bool)]))
+
+"""
+
+Here, we obtained a noise standard deviation of 7.26. For comparison, a simple 
+standard deviation of all voxels in the estimated mask (as done in the previous 
+example :ref:`example_snr_in_cc`) gives a value of 6.1.
+
+"""
+
+"""
+
+.. [Koay2009] Koay C.G., E. Ozarslan, C. Pierpaoli. Probabilistic Identification and Estimation of Noise (PIESNO): A self-consistent approach and its applications in MRI. JMR, 199(1):94-103, 2009.
+
+.. include:: ../links_names.inc
+
+
+"""

doc/examples/snr_in_cc.py

@@ -130,6 +130,7 @@
 nib.save(mask_noise_img, 'mask_noise.nii.gz')
 
 noise_std = np.std(data[mask_noise, :])
+print('Noise standard deviation sigma= ', noise_std)
 
 """We can now compute the SNR for each DWI. For example, report SNR
 for DW images with gradient direction that lies the closest to

doc/examples_index.rst

@@ -30,11 +30,15 @@ Brain Extraction
 
 - :ref:`example_brain_extraction_dwi`
 
-SNR estimation
-~~~~~~~~~~~~~~
+Basic SNR estimation
+~~~~~~~~~~~~~~~~~~~~
 
 - :ref:`example_snr_in_cc`
 
+PIESNO
+~~~~~~
+- :ref:`example_piesno_sherbrooke`
+
 Denoising
 ~~~~~~~~~