Skip to content
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.

Sign up for GitHub

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

Jump to bottom

FIX: PEP8 in denoise #957

Merged
merged 5 commits into from Mar 18, 2016
Merged

FIX: PEP8 in denoise #957

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
98024a1
PEP8 fixed in denoise
manu-tej Mar 1, 2016
b4b7829
Update test_noise_estimate.py
manu-tej Mar 1, 2016
00cd4e4
Update test_noise_estimate.py
manu-tej Mar 7, 2016
82d6f0f
Update test_noise_estimate.py
manu-tej Mar 8, 2016
c0c55de
Update noise_estimate.py
manu-tej Mar 12, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
46 changes: 27 additions & 19 deletions dipy/denoise/noise_estimate.py
View file
Open in desktop
Expand Up @@ -24,7 +24,8 @@ def _inv_nchi_cdf(N, K, alpha):
128: 0.5322811923303339}


def piesno(data, N, 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).

Expand All @@ -36,8 +37,8 @@ def piesno(data, N, alpha=0.01, l=100, itermax=100, eps=1e-5, return_mask=False)

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 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.

Expand Down Expand Up @@ -71,7 +72,8 @@ def piesno(data, N, alpha=0.01, l=100, itermax=100, eps=1e-5, return_mask=False)
------
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.
along the last axis, i.e. dMRI or fMRI data, not structural data like
T1/T2.

This function processes the data slice by slice, as originally designed in
the paper. Use it to get a slice by slice estimation of the noise, as in
Expand Down Expand Up @@ -103,15 +105,17 @@ def piesno(data, N, alpha=0.01, l=100, itermax=100, eps=1e-5, return_mask=False)
q = 0.5

# Initial estimation of sigma
initial_estimation = np.percentile(data, q * 100) / np.sqrt(2 * _inv_nchi_cdf(N, 1, q))
initial_estimation = (np.percentile(data, q * 100) /
np.sqrt(2 * _inv_nchi_cdf(N, 1, q)))

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,
sigma[idx], mask_noise[..., idx] = _piesno_3D(data[..., idx, :],
N,
alpha=alpha,
l=l,
itermax=itermax,
Expand All @@ -120,7 +124,8 @@ def piesno(data, N, alpha=0.01, l=100, itermax=100, eps=1e-5, return_mask=False)
initial_estimation=initial_estimation)

else:
sigma, mask_noise = _piesno_3D(data, N,
sigma, mask_noise = _piesno_3D(data,
N,
alpha=alpha,
l=l,
itermax=itermax,
Expand Down Expand Up @@ -185,7 +190,8 @@ def _piesno_3D(data, N, alpha=0.01, l=100, itermax=100, eps=1e-5,
------
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.
along the last axis, i.e. dMRI or fMRI data, not structural data like
T1/T2.

References
------------
Expand Down Expand Up @@ -231,11 +237,11 @@ def _piesno_3D(data, N, alpha=0.01, l=100, itermax=100, eps=1e-5,
lambda_minus = _inv_nchi_cdf(N, K, alpha/2)
lambda_plus = _inv_nchi_cdf(N, K, 1 - alpha/2)


for sigma_init in phi:

s = sum_m2 / (2 * K * sigma_init**2)
found_idx = np.sum(np.logical_and(lambda_minus <= s, s <= lambda_plus), dtype=np.int16)
found_idx = np.sum(np.logical_and(lambda_minus <= s, s <= lambda_plus),
dtype=np.int16)

if found_idx > prev_idx:
sigma = sigma_init
Expand Down Expand Up @@ -280,7 +286,8 @@ def estimate_sigma(arr, disable_background_masking=False, N=0):
Number of coils of the receiver array. Use N = 1 in case of a SENSE
reconstruction (Philips scanners) or the number of coils for a GRAPPA
reconstruction (Siemens and GE). Use 0 to disable the correction factor,
as for example if the noise is Gaussian distributed. See [1] for more information.
as for example if the noise is Gaussian distributed. See [1] for more
information.

Returns
-------
Expand All @@ -295,19 +302,19 @@ def estimate_sigma(arr, disable_background_masking=False, N=0):
(see [1]_, equation 18) with theta = 0.
Since this function was introduced in [2]_ for T1 imaging,
it is expected to perform ok on diffusion MRI data, but might oversmooth
some regions and leave others un-denoised for spatially varying noise profiles.
Consider using :func:`piesno` to estimate sigma instead if visual inacuracies
are apparent in the denoised result.
some regions and leave others un-denoised for spatially varying noise
profiles. Consider using :func:`piesno` to estimate sigma instead if visual
inacuracies are apparent in the denoised result.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Typo (was there before, but maybe we fix it here?):

"inacuracies"=> "inaccuracies"


Reference
-------
.. [1] Koay, C. G., & Basser, P. J. (2006). Analytically exact correction
scheme for signal extraction from noisy magnitude MR signals.
Journal of Magnetic Resonance), 179(2), 317-22.

.. [2] Coupe, P., Yger, P., Prima, S., Hellier, P., Kervrann, C., Barillot, C., 2008.
An optimized blockwise nonlocal means denoising filter for 3-D magnetic
resonance images, IEEE Trans. Med. Imaging 27, 425-41.
.. [2] Coupe, P., Yger, P., Prima, S., Hellier, P., Kervrann, C., Barillot,
C., 2008. An optimized blockwise nonlocal means denoising filter for 3-D
magnetic resonance images, IEEE Trans. Med. Imaging 27, 425-41.
"""
k = np.zeros((3, 3, 3), dtype=np.int8)

Expand All @@ -318,8 +325,9 @@ def estimate_sigma(arr, disable_background_masking=False, N=0):
k[1, 1, 0] = 1
k[1, 1, 2] = 1

# Precomputed factor from Koay 2006, this corrects the bias of magnitude image
correction_factor = {0: 1, # No correction
# Precomputed factor from Koay 2006, this corrects the bias of magnitude
# image
correction_factor = {0: 1, # No correction
1: 0.42920367320510366,
4: 0.4834941393603609,
6: 0.4891759468548269,
Expand Down
1 change: 1 addition & 0 deletions dipy/denoise/tests/test_denoise.py
View file
Open in desktop
Expand Up @@ -5,6 +5,7 @@
import dipy.data as dpd
import nibabel as nib


def test_denoise():
"""

Expand Down
65 changes: 42 additions & 23 deletions dipy/denoise/tests/test_noise_estimate.py
View file
Open in desktop
Expand Up @@ -5,10 +5,13 @@

from numpy.testing import (assert_almost_equal, assert_equal, assert_,
assert_array_almost_equal)
from dipy.denoise.noise_estimate import _inv_nchi_cdf, piesno, estimate_sigma, _piesno_3D
from dipy.denoise.noise_estimate import _inv_nchi_cdf, piesno, estimate_sigma
from dipy.denoise.noise_estimate import _piesno_3D
import dipy.data

# See page 5 of the reference paper for tested values


def test_inv_nchi():
# Values taken from hispeed.MedianPIESNO.lambdaPlus
# and hispeed.MedianPIESNO.lambdaMinus
Expand All @@ -27,13 +30,15 @@ def test_piesno():
# Values taken from hispeed.OptimalPIESNO with the test data
# in the package computed in matlab
test_piesno_data = nib.load(dipy.data.get_data("test_piesno")).get_data()
sigma = piesno(test_piesno_data, N=8, alpha=0.01, l=1, eps=1e-10, return_mask=False)
sigma = piesno(test_piesno_data, N=8, alpha=0.01, l=1, eps=1e-10,
return_mask=False)
assert_almost_equal(sigma, 0.010749458025559)

noise1 = (np.random.randn(100, 100, 100) * 50) + 10
noise2 = (np.random.randn(100, 100, 100) * 50) + 10
rician_noise = np.sqrt(noise1**2 + noise2**2)
sigma, mask = piesno(rician_noise, N=1, alpha=0.01, l=1, eps=1e-10, return_mask=True)
sigma, mask = piesno(rician_noise, N=1, alpha=0.01, l=1, eps=1e-10,
return_mask=True)

# less than 3% of error?
assert_(np.abs(sigma - 50) / sigma < 0.03)
Expand All @@ -49,53 +54,63 @@ def test_piesno():
assert_(np.abs(sigma - 50) / sigma < 0.03)

sigma = _piesno_3D(rician_noise, N=1, alpha=0.01, l=1, eps=1e-10,
return_mask=False,
initial_estimation=initial_estimation)
return_mask=False,
initial_estimation=initial_estimation)
assert_(np.abs(sigma - 50) / sigma < 0.03)

sigma = _piesno_3D(np.zeros_like(rician_noise), N=1, alpha=0.01, l=1, eps=1e-10,
return_mask=False,
initial_estimation=initial_estimation)
sigma = _piesno_3D(np.zeros_like(rician_noise), N=1, alpha=0.01, l=1,
eps=1e-10, return_mask=False,
initial_estimation=initial_estimation)

assert_(np.all(sigma == 0))

sigma, mask = _piesno_3D(np.zeros_like(rician_noise), N=1, alpha=0.01, l=1, eps=1e-10,
return_mask=True,
initial_estimation=initial_estimation)
sigma, mask = _piesno_3D(np.zeros_like(rician_noise), N=1, alpha=0.01, l=1,
eps=1e-10, return_mask=True,
initial_estimation=initial_estimation)

assert_(np.all(sigma == 0))
assert_(np.all(mask == 0))

# Check if no noise points found in array it exits
sigma = _piesno_3D(1000*np.ones_like(rician_noise), N=1, alpha=0.01, l=1, eps=1e-10,
return_mask=False, initial_estimation=10)
sigma = _piesno_3D(1000*np.ones_like(rician_noise), N=1, alpha=0.01, l=1,
eps=1e-10, return_mask=False, initial_estimation=10)
assert_(np.all(sigma == 10))


def test_estimate_sigma():

sigma = estimate_sigma(np.ones((7, 7, 7)), disable_background_masking=True)
assert_equal(sigma, 0.)

sigma = estimate_sigma(np.ones((7, 7, 7, 3)), disable_background_masking=True)
sigma = estimate_sigma(np.ones((7, 7, 7, 3)),
disable_background_masking=True)
assert_equal(sigma, np.array([0., 0., 0.]))

sigma = estimate_sigma(5 * np.ones((7, 7, 7)), disable_background_masking=False)
sigma = estimate_sigma(5 * np.ones((7, 7, 7)),
disable_background_masking=False)
assert_equal(sigma, 0.)

sigma = estimate_sigma(5 * np.ones((7, 7, 7, 3)), disable_background_masking=False)
sigma = estimate_sigma(5 * np.ones((7, 7, 7, 3)),
disable_background_masking=False)
assert_equal(sigma, np.array([0., 0., 0.]))

arr = np.zeros((3, 3, 3))
arr[0, 0, 0] = 1
sigma = estimate_sigma(arr, disable_background_masking=False, N=1)
assert_array_almost_equal(sigma, 0.10286889997472792 / np.sqrt(0.42920367320510366))
assert_array_almost_equal(sigma,
(0.10286889997472792 /
np.sqrt(0.42920367320510366)))

arr = np.zeros((3, 3, 3, 3))
arr[0, 0, 0] = 1
sigma = estimate_sigma(arr, disable_background_masking=False, N=1)
assert_array_almost_equal(sigma, np.array([0.10286889997472792 / np.sqrt(0.42920367320510366),
0.10286889997472792 / np.sqrt(0.42920367320510366),
0.10286889997472792 / np.sqrt(0.42920367320510366)]))
assert_array_almost_equal(sigma,
np.array([0.10286889997472792 /
np.sqrt(0.42920367320510366),
0.10286889997472792 /
np.sqrt(0.42920367320510366),
0.10286889997472792 /
np.sqrt(0.42920367320510366)]))
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@samuelstjean Should I even change here?


arr = np.zeros((3, 3, 3))
arr[0, 0, 0] = 1
Expand All @@ -110,6 +125,10 @@ def test_estimate_sigma():

arr[0, 0, 0] = 1
sigma = estimate_sigma(arr, disable_background_masking=True, N=12)
assert_array_almost_equal(sigma, np.array([0.46291005 / np.sqrt(0.4946862482541263),
0.46291005 / np.sqrt(0.4946862482541263),
0.46291005 / np.sqrt(0.4946862482541263)]))
assert_array_almost_equal(sigma,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is hard to read, no need to break lines for the 80 characters rule everytime if it does not really help.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@samuelstjean Should I change it to

assert_array_almost_equal(sigma,
                          np.array([0.46291005 / np.sqrt(0.4946862482541263),
                                    0.46291005 / np.sqrt(0.4946862482541263),
                                    0.46291005 / np.sqrt(0.4946862482541263)]))

or

assert_array_almost_equal(sigma, np.array([0.46291005 / np.sqrt(0.4946862482541263), 0.46291005 /  np.sqrt(0.4946862482541263), 0.46291005 / np.sqrt(0.4946862482541263)]))

np.array([0.46291005 /
np.sqrt(0.4946862482541263),
0.46291005 /
np.sqrt(0.4946862482541263),
0.46291005 /
np.sqrt(0.4946862482541263)]))