-
Notifications
You must be signed in to change notification settings - Fork 437
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
WIP: New PIESNO example and small corrections #390
Changes from 7 commits
c97c7ad
90934d6
8a571bc
d4b1903
8c31108
6f930ff
3328583
6349089
4757856
ed23a57
5bc2ad9
eee2a9d
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -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,97 @@ 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) | ||
|
||
# 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) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. So this is where you are suggesting returning There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Note to myself - change this to return the full array, as it's just plain misleading to return a single value on real data. |
||
#print(sigma, num) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. this is probably a debug statement (?). You can remove it or mark it with a |
||
|
||
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): | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Another nitpick by myself : The real version of the paper is used on 2D images as opposed to 3D volume like it is now made in the example and the new function. The mode/median part is just a hack I added myself for simplicuity and does not rely on any real theoretical background. I try to specify that more explicitely in the functions, for example by removing the _ in the piesno_3D function. Anyway, is it more of a convenience than anything else and seems to work so far on all our tests, but for spinal cord imaging for example, it is better to use the slice by slice version because of the huge gaps in the Z axis, which leads to a high variability in the standard deviation estimate of each slice. |
||
""" | ||
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. | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Remark from line 32 to 36 could also be added here. |
||
|
||
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 | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Just out of curiosity (and the kind of data I have...): if this is not the case, for example if the data in the background gets masked out by the scanner, is it possible to use some of the brain data for this? For example, is it possible to use the ventricles? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Best way to know is to try it, the method is fully automatic so it will pick what it can. Sometimes it will find the skull and being a region of noise, but not the ventricules from what I have seen. Also, if there is no background, it might just not find anyhting qualifying as background/noise either. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Another question about this: have you run this on the HCP data? What N are you putting in for that? cc: @klchan13 There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Not personally for the public released datasets (I'm playing with a spinal cord one currently). Do they say the hardware specs of their scanner somewhere? They have a ton of papers about that, for example http://onlinelibrary.wiley.com/doi/10.1002/mrm.24623/abstract;jsessionid=1A58CFCAB23677366853EE1BC46F011A.d03t02 They used a SENSE based one on a 32 channel coils, so if they always used that one, and this is the right human connectome project (there are two of them), that would mean N=1 because of SENSE. Another way to check that is just draw a manual roi, then check the nosie profile. It's easy to see if it's rician or not by checking the curve. In the latter case, it's harder to guess the value of N (anyway I can't). |
||
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 +197,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 +234,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] |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,100 @@ | ||
""" | ||
============================= | ||
Noise estimation using PIESNO | ||
============================= | ||
|
||
Using PIESNO [Koay2009]_ one can detect the noise's standard deviation from | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. => "the standard deviation of the noise" |
||
diffusion-weighted imaging (DWI). PIESNO also works from multiple channel | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. => "...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 mathetically quite involved. | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. typo : "mathematically" |
||
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, as opposed to tissue voxels that have diffuison intensities | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. diffuison -> diffusion |
||
that vary quite a lot across different directions. | ||
2) From these esimated background voxels and | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. typo |
||
the input number of coils N, PIESNO finds what sigma each Gaussian distributed | ||
image profile 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 glory details. | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I think you mean "gory", but I would just say "all the details". |
||
|
||
PIESNO makes an important assumption: the | ||
Gaussian noise standard deviation is assumed to be uniform either | ||
across multiple slice locations or across multiple images of the same location, | ||
e.g., if the readout bandwidth is maintained at the same level for all the images. | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Could you break this sentence up a little bit? It's a bit hard to parse. Imagine that you are talking to a psychologist, and you want them to understand what this means :-) |
||
|
||
""" | ||
|
||
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 N | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. the right number |
||
of coil used to acquire this dataset. It is also important to know what | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. typos: with the right number of coils |
||
was the parallel reconstruction algorithm used. | ||
Here, the data comes from a GRAPPA reconstruction from | ||
a 12-element head coil available on the Tim Trio Siemens, for which | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. 12-elementS |
||
the 12 coil elements are combined into 4 groups of 3 coil elements | ||
each. These groups are received through 4 distinct receiver channels, | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Are the group received? Seems weird, it's more like the signal is received through 4 distinct groups of receiver channel or something like that no? |
||
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 i in range(2): | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. for a in ax: |
||
ax[i].set_axis_off() | ||
|
||
#plt.show() | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. If it's commented out, maybe it doesn't need to be here? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I like to keep it there for people not familiar with fvtk. They see what is needed to actually see the window without saving the PNG. I have added a comment. But I can remove if people don't like it. |
||
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) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. What previous example? Could you please link to that one with the |
||
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 | ||
|
||
|
||
""" |
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -30,11 +30,15 @@ Brain Extraction | |
|
||
- :ref:`example_brain_extraction_dwi` | ||
|
||
SNR estimation | ||
~~~~~~~~~~~~~~ | ||
Basic SNR estimation | ||
~~~~~~~~~~~~~~~~~~~~ | ||
|
||
- :ref:`example_snr_in_cc` | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Is this the example you are referring to? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Yes. |
||
|
||
PIESNO | ||
~~~~~~ | ||
- :ref:`example_piesno_sherbrooke` | ||
|
||
Denoising | ||
~~~~~~~~~ | ||
|
||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It seems to me that the mode might fail in some cases, depending on the tolerance of the mode calculation. Is this what Koay et al. did? How about the median, for robustness?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I just used the mode because most of the time you have something like 60%-75% of the slices with exactly the same sigma. The median could also be used I guess, I personaly sued it to double check that the mode worked correctly.
I don't think they actually to use the same sigma value for the whole volume, it's just really less of an hassle than to consider each slice separately for any purpose one might have afterward.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
One again, the slice by slice in a for loop is mostly for the end user convenience, and not backed in any way by theory. One can call the _piesno_3D function himself to get a slcie by slice estimate (as desired in spinal cord dMRI), so I'll think a a way to write that as a note in the docstring.