Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP

Loading…

Noise #26

Merged
merged 24 commits into from

7 participants

@arokem
Owner

This pull request refactors the noise-generating functions in the simulation module. I have written down precisely how I define SNR for the purpose of this. The one thing that I am feeling unsure about is the definition of the Rician noise from two Gaussians. My idea there follows up on Bago's suggestion on the mailing list.

dipy/sims/phantom.py
((111 lines not shown))
+ # all the voxels and with the right shape:
+ noise = np.random.randn(*vol.shape) * (np.sqrt(p_noise))
+ elif noise_type == 'rician':
+ # To generate rician noise, we add two IID Gaussian noise sources in
+ # the complex domain and combine them together:
+ noise1 = np.random.randn(*vol.shape)
+ noise2 = np.random.randn(*vol.shape)
+ noise_initial = np.sqrt(noise1**2 + noise2**2)
+ # Now let's get control of the variance, to make sure that we have the
+ # right power:
+ var_initial = np.var(noise_initial, -1)
+
+ # We will scale a demeaned version of the noise
+ mean_initial = np.mean(noise_initial,-1)[...,np.newaxis]
+ demeaned = noise_initial - mean_initial
+ # By our goal for the variance:
@MrBago Owner
MrBago added a note

If I'm reading this code correctly you're assuming that abs(signal + noise) is the same as signal + abs(noise), which I believe isn't the case. To get a rician distribution I believe you need abs(signal + noise) (otherwise I believe it's a Rayleigh distribution).

@arokem Owner
arokem added a note

Which line does this refer to? Both the signal and the noise here (in the Rician case) are supposed to be strictly positive, so, as I understand it abs(signal + noise) is the same as signal + abs(noise).

@MrBago Owner
MrBago added a note

Sorry, this was meant to go on line 199, I'm not sure why it ended up here. signal + abs(complex) is not the same as abs(signal + complex), the mr noise is complex, that is why we end up with a rician distribution.

@arokem Owner
arokem added a note

So, if I understand you correctly, you suggest that we corrupt the signal with something like:

signal = signal + np.random.randn(shape) + np.random.randn(shape) * 1j

And then take the absolute value of that as our final signal?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
dipy/sims/phantom.py
((105 lines not shown))
+ # SNR = var(signal)/var(noise) => var(noise) = var(signal)/SNR:
+ p_noise = np.mean(p_signal/snr)
+
+ if noise_type == 'gaussian':
+ # Generate the noise with the correct standard deviation, averaged over
+ # all the voxels and with the right shape:
+ noise = np.random.randn(*vol.shape) * (np.sqrt(p_noise))
+ elif noise_type == 'rician':
+ # To generate rician noise, we add two IID Gaussian noise sources in
+ # the complex domain and combine them together:
+ noise1 = np.random.randn(*vol.shape)
@MrBago Owner
MrBago added a note

Is there a reason we're getting the variance from the noise sample and not the distribution? This could be a real issue if we're adding noise to a small number of voxels, (will this work for one voxel? will the variance be 0?)

@arokem Owner
arokem added a note

Do you mean we should get the distribution from scipy.stats and sample from that? They way I wrote this, we are getting a sample (or samples, for the Rician case) of Gaussian noise and then all the rest is just to get it to have the right variance, so that:

SNR == var(signal)/var(noise)

@MrBago Owner
MrBago added a note

I think we need to distinguish between var(noise-sample) and var(noise-distribution). These two will be very close for a large noise-sample, but at the end of the day, var(noise-sample) is only an estimate of var(noise-distribution) and when noise-sample is small, it's a pretty bad estimate. I believe is usually defined as SNR = mean(signal) / sqrt(var(noise-distribution)). We can either hard code the analytical form for var(rice-distribution) for use scipy.stats either should work. Also the variance of a rice distribution is very close to the variance of a normal distribution with the same sigma (sigma**2), maybe that estimate is close enough for what we're doing.

Bago

@arokem Owner
arokem added a note
@arokem Owner
arokem added a note
@MrBago Owner
MrBago added a note

We should set up a skype or something to talk about this, It's not clear to me why you're using var(signal) / var(noise). Also it looks like you're implementation of the Gaussian noise scales the noise differently for each voxel. Even though there might be applications where this is useful, it's not physical for MRI. MRI noise generally has the same power over the whole volume.

Also, in the MRI literature there are two ways that noise is usually measured, either the mean signal to noise of diffusion weighted images or the mean signal to noise of b0(no diffusion weighting) images. I believe that the first is a more useful number, but because the SNR of dwi images depends on b-value some people report the second. I believe we should keep that in mind.

@arokem Owner
arokem added a note
@MrBago Owner
MrBago added a note
@arokem Owner
arokem added a note
@MrBago Owner
MrBago added a note

Ya, the std(rician) isn't exactly sigma, it's pretty close to sigma, especially as sigma gets small (snr gets big). The exact form is on Wikipedia.

@arokem Owner
arokem added a note
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@arokem
Owner

Or maybe a factor sqrt(7/3) ?!

@arokem
Owner

Any more comments here? Otherwise, I will go ahead and merge this some time during the weekend

@MrBago
Owner

Thanks for the hard work ariel, I wanted to give this a little more thought. I think we should be able to solve for the scaling analytically instead of approximating it iteratively, I'm also wondering if we should be using rms(noise) instead of std(noise). I just need to set aside a little time to think about it. If you want, we can merge it and I can propose changes when I have time.

@arokem
Owner
@arokem
Owner

I have thought a bit more about the calculation of SNR for DWI: it doesn't really make sense to calculate what is the equivalent of the tSNR that people calculate for fMRI (which is what we have here right now), because part of the variance in the signal is actually what we are interested in. An interesting idea that I got from @rfdougherty is to use the b0 measurements to calculate the noise and to use the Stejskal/Tanner equation to then predict what the signal in each voxel should be like based on the b0 measurement, the b value used in the measurement of the diffusion-weighted images and the mean-diffusivity in each voxel, calculated using DTI from the data.

Essentially, the calculating would like something like this:

sig = np.mean(data[..., b0_idx], -1)
noise = np.var(data[..., bo_idx], -1)
snr = (sig/noise) * np.exp( -b * mean_diffusivity)

Does that make sense to you all?

@arokem
Owner

Refactor plan:

  1. Specify sigma instead of the SNR for the noise-generating function itself.

  2. Wrap that function in another API that can get a gradient table as its input and use the SNR relative to the b0.

    For the SNR: use the RMS of the noise, which is not mean-subtracted (in contrast to the variance).

    So: SNR = mean(true signal)/RMS(true added noise).

    For Gaussian, this will be equivalent to taking the std, but for Rician it might be different.

    As an alternative, consider setting SNR to be: mean/sigma

  3. To apply SNR estimates from actual data, we might want to be able to estimate the parameters of the Rician distribution from lower order moments of the distribution of b0, or (probably better), use an ROI where the Rician is approximately equal to the Gaussian, so that you can use std.

@travisbot

This pull request fails (merged ce322ef6 into 59eab5c).

@travisbot

This pull request fails (merged 560508d2 into 898a401).

@travisbot

This pull request fails (merged 65419341 into 4ff1bd1).

@travisbot

This pull request fails (merged 969efcb5 into 38c1ec1).

@travisbot

This pull request fails (merged 8f3048a6 into 38c1ec1).

@travisbot

This pull request fails (merged 6fa78d2e into 38c1ec1).

@travisbot

This pull request fails (merged 6fa78d2e into 38c1ec1).

@arokem
Owner

OK - I think that this should fulfill what we talked about a few weeks ago. @MrBago - what do you think?

@travisbot

This pull request fails (merged ca3f0ef1 into 38c1ec1).

@travisbot

This pull request fails (merged ca3f0ef1 into 38c1ec1).

@travisbot

This pull request passes (merged 80fd4bf3 into dc8b7fd).

dipy/sims/phantom.py
@@ -62,11 +64,12 @@ def orbital_phantom(bvals=None,
angles and radii define the total thickness options
S0 : double, simulated signal without diffusion gradients applied
Default 100.
- snr : signal to noise ratio
- Used for applying rician noise to the data.
- Default 200. Common is 20.
- background_noise : boolean, Default False
-
+ snr : float, optional
+ The signal to noise ratio sed to apply Rician noise to the data.
@stefanv Owner
stefanv added a note

sed?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
dipy/sims/phantom.py
((9 lines not shown))
+ Gudbjartsson, H and Patz, S (2008). The Rician Distribution of Noisy MRI
@stefanv Owner
stefanv added a note
References
---------------
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
dipy/sims/phantom.py
@@ -74,7 +77,16 @@ def orbital_phantom(bvals=None,
Notes
--------
Crossings can be created by adding multiple orbitual_phantom outputs.
+
+ In these simulations, we can ask for Rician noise to be added. In
+ that case, the definition of SNR is as follows:
+
+ SNR = mean(true signal)/RMS(true added noise).
@stefanv Owner
stefanv added a note

Std = RMS only under assumption of zero-mean noise. Is this the case?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
dipy/sims/phantom.py
@@ -74,7 +77,16 @@ def orbital_phantom(bvals=None,
Notes
--------
Crossings can be created by adding multiple orbitual_phantom outputs.
@stefanv Owner
stefanv added a note

orbitual -> orbital

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
dipy/sims/phantom.py
@@ -84,7 +96,7 @@ def f(t):
z=np.linspace(-1,1,len(x))
return x,y,z
- data=orbitual_phantom(func=f)
+ data=orbital_phantom(func=f)
@stefanv Owner
stefanv added a note

data = orbital_phantom(func=f)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
dipy/sims/phantom.py
((66 lines not shown))
- return np.sqrt(K/np.float(snr))*np.random.randn(*vol.shape)
- noise1=gaussian_noise(vol)
- noise2=gaussian_noise(vol)
- return np.sqrt((vol+noise1)**2+noise2**2)
-
-def add_gaussian_noise(vol,snr=20):
- """ add gaussian noise in a 4D array with a specific snr
-
- Parameters
- -----------
- vol : array, shape (X,Y,Z,W)
- snr : float,
- signal to noise ratio
-
+
+ sigma: float
@stefanv Owner
stefanv added a note

sigma : float, and elsewhere

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
dipy/sims/phantom.py
@@ -131,76 +143,127 @@ def f(t):
#FA=ten.fa()
#FA[np.isnan(FA)]=0
#vol[np.isnan(vol)]=0
+
+ if snr is not None:
+ # We start by guessing that sigma should be approximately such that the
+ # snr is right and using: snr = mean_sig/rms => rms = mean_sig/snr
+ mean_sig = np.mean(vol)
+ sigma = mean_sig/snr
+ vol_w_noise = add_noise(vol, sigma, noise_type='rician')
+ noise = vol - vol_w_noise
@stefanv Owner
stefanv added a note

No gripes with this PR, so I'll just comment in general on the dangers of floating point subtraction and losing precision:

s = np.ones((10))*1e10; n = np.random.random((10)) * 1e-5; sn = s+n
np.mean(sn - s), np.mean(n)

However, since this situation would only occur if you have very small noise levels, it really doesn't apply here :)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
dipy/sims/phantom.py
@@ -131,76 +143,127 @@ def f(t):
#FA=ten.fa()
#FA[np.isnan(FA)]=0
#vol[np.isnan(vol)]=0
+
+ if snr is not None:
+ # We start by guessing that sigma should be approximately such that the
+ # snr is right and using: snr = mean_sig/rms => rms = mean_sig/snr
+ mean_sig = np.mean(vol)
+ sigma = mean_sig/snr
+ vol_w_noise = add_noise(vol, sigma, noise_type='rician')
@MrBago Owner
MrBago added a note

We can replace this with vol_w_noise = add_noise(vol, sigma, noise_type='rician', scale=1./scipy.stats.rice(snr).std())

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
dipy/sims/phantom.py
@@ -131,76 +143,127 @@ def f(t):
#FA=ten.fa()
#FA[np.isnan(FA)]=0
#vol[np.isnan(vol)]=0
+
+ if snr is not None:
+ # We start by guessing that sigma should be approximately such that the
+ # snr is right and using: snr = mean_sig/rms => rms = mean_sig/snr
+ mean_sig = np.mean(vol)
+ sigma = mean_sig/snr
@MrBago Owner
MrBago added a note

Sorry, I don't know what I'm smoking.

I think in order to avoid the while loop we can replace sigma = mean_sig/snr with sigma = mean_sig / snr / scipy.stats.rice(snr).std(). The only problem is that scipy.stats.rice wiggs out if snr == 0 or > 37.5 so we need to do something like:

if snr > 37.5:
    scale = 1
elif snr == 0:
    scale = np.sqrt((4 - np.pi)/2)
else:
    scale = scipy.stats.rice(snr).std()

sigma = mean_sig / snr / scale

I'm kinda tiered so someone else should have a second look at that and make sure it makes sense.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@travisbot

This pull request fails (merged 552a13c3 into dc8b7fd).

@arokem
Owner

@MrBago: is this what you meant? As travisbot was quick to point out, test_snr now fails with this. First, do you think the test is OK? Second, I debugged on this and for the cases used in the test, the scale variable is very close to 1, so this block of code doesn't really do much. That is, you end up with a sigma very similar to the one you would have gotten from mean_sig/snr. Unfortunately, for the Rician case, this is pretty far off from the sigma you actually want. Note that sigma is not the standard deviation of the Rician, but is a bit off of that (because of what the code in add_noise does with 'sigma').

@arokem
Owner

@MrBago - does this make sense to you? We had to add a 'fudge factor' of sqrt(2) to make the SNR numbers come out right here (that last commit). But why sqrt(2)? Any ideas?

@coryahrens

sounds like the length of a unit complex number...

@MrBago
Owner

You have sqrt(x / 2) * sqrt(2) which is sqrt(x) so I'm not sure what you're doing here.

@MrBago
Owner

I'm apprehensive about the way this function is set up. MRI noise is poorly described by "std" or "variance", and will have a different variance in different parts of the image and in different images, (ie different b-values in a dwi acquisition). The sigma is a much better description and should be relatively consistent, at least within one image (and hopefully even at different b-value acquisitions).

@arokem
Owner
@arokem
Owner
@MrBago
Owner

The problem is that you're thinking of noise as (signal + noise). In MRI measured signal is abs(signal + complex_noise). This means that var(measured - signal) does not have a consistent distribution independent of signal. The fudge factor comes from assuming (signal == 0), which you can do outside the brain.

@MrBago
Owner

Calculating a different sigma in each voxel is unphysical, it's an interesting exercise, but MRI scanners will generate noise with a consistent sigma, and an inconsistent variance or "SNR". If you don't believe me we should set up a phantom simulation.

@arokem
Owner
@MrBago
Owner

From Maxime Descoteaux's 2007 qball paper

The noise is
generated with a complex Gaussian noise with a standard
deviation of σ, producing a S0 signal with SNR = 1/σ , that
is we define SNR as the ratio of maximum signal intensity
of S0 to the standard deviation σ of the complex Gaussian
noise.

ie noise is defined as b0-signal / sigma.

@MrBago
Owner

Sorry I mean SNR is defined as b0-signal / sigma.

@arokem
Owner
@arokem
Owner
@arokem
Owner
@MrBago
Owner
@arokem
Owner
@MrBago
Owner
@arokem
Owner
@arokem arokem RF: Moved noise-addition functions to the voxel sub-module.
Also: adopted the convention used by Descoteaux et al. to define SNR.
e7200c9
@coryahrens

Could you please explain why MRI images are not well characterized by mean and variance. My understanding of Rician noise was that it covers MRI pretty well. An unusual feature of this noise is that it is non-stationary, i.e., the moments (mean, variance, etc) are not constant; they change with the value of the signal.

Isn't the standard deviation of the Gaussian noise for each channel fairly constant over the whole image? In this case, the sigma in the Rician distribution would be constant, but the signal value would be changing and hence the standard deviation of the whole Rician distribution changes over the volume.

@MrBago
Owner

That's a good way of describing it. I just meant that because the Rician distribution is non-stationary, we expect variance(signal) to vary over a volume, even though variance(complex_noise) is pretty consistent over the volume.

@arokem
Owner

OK @Garyfallidis, this is ready for you to take a look. In particular, there's still one test failing for the orbital phantom. We're not getting the FA you would expect for the single tensor voxel. Any ideas?

@arokem
Owner

Maybe this is because there is no single-tensor voxel anywhere in the volume?

@Garyfallidis

Okay Travis reports more than one errors here.

1077FAIL: Doctest: dipy.sims.phantom.orbital_phantom
1078----------------------------------------------------------------------
1079Traceback (most recent call last):
1080 File "/usr/lib/python2.7/doctest.py", line 2201, in runTest
1081 raise self.failureException(self.format_failure(new.getvalue()))
1082AssertionError: Failed doctest test for dipy.sims.phantom.orbital_phantom
1083 File "/usr/local/lib/python2.7/dist-packages/dipy/sims/phantom.py", line 92, in orbital_phantom
1084
1085----------------------------------------------------------------------
1086File "/usr/local/lib/python2.7/dist-packages/dipy/sims/phantom.py", line 142, in dipy.sims.phantom.orbital_phantom
1087Failed example:
1088 data = orbital_phantom(func=f)
1089Exception raised:
1090 Traceback (most recent call last):
1091 File "/usr/lib/python2.7/doctest.py", line 1289, in __run
1092 compileflags, 1) in test.globs
1093 File "", line 1, in
1094 data = orbital_phantom(func=f)
1095 File "/usr/local/lib/python2.7/dist-packages/dipy/sims/phantom.py", line 146, in orbital_phantom
1096 if gtab.bvals is None:
1097 AttributeError: 'NoneType' object has no attribute 'bvals'

In phantom.py line 146. if gtab.bvals is None should be if gtab is None.

1035ERROR: dipy.reconst.tests.test_dti.test_TensorModel
1036----------------------------------------------------------------------
1037Traceback (most recent call last):
1038 File "/usr/lib/python2.7/dist-packages/nose/case.py", line 197, in runTest
1039 self.test(*self.arg)
1040 File "/usr/local/lib/python2.7/dist-packages/dipy/reconst/tests/test_dti.py", line 22, in test_TensorModel
1041 assert_equal(dtifit.fa < 0.5, True)
1042 File "/usr/local/lib/python2.7/dist-packages/dipy/core/onetime.py", line 175, in get
1043 val = self.getter(obj)
1044 File "/usr/local/lib/python2.7/dist-packages/dipy/reconst/dti.py", line 143, in fa
1045 all_zero = (np.isclose(ev1, 0) & np.isclose(ev2, 0) & np.isclose(ev3, 0))
1046AttributeError: 'module' object has no attribute 'isclose'

Which version of numpy introduces isclose? It is not supported in 1.6.1 which at least I use.

1100======================================================================
1101FAIL: Doctest: dipy.sims.voxel.add_noise
1102----------------------------------------------------------------------
1103Traceback (most recent call last):
1104 File "/usr/lib/python2.7/doctest.py", line 2201, in runTest
1105 raise self.failureException(self.format_failure(new.getvalue()))
1106AssertionError: Failed doctest test for dipy.sims.voxel.add_noise
1107 File "/usr/local/lib/python2.7/dist-packages/dipy/sims/voxel.py", line 48, in add_noise
1108
1109----------------------------------------------------------------------
1110File "/usr/local/lib/python2.7/dist-packages/dipy/sims/voxel.py", line 86, in dipy.sims.voxel.add_noise
1111Failed example:
1112 signal_w_noise = add_noise(signal, snr=10, noise_type='rician')
1113Exception raised:
1114 Traceback (most recent call last):
1115 File "/usr/lib/python2.7/doctest.py", line 1289, in __run
1116 compileflags, 1) in test.globs
1117 File "", line 1, in
1118 signal_w_noise = add_noise(signal, snr=10, noise_type='rician')
1119 TypeError: add_noise() takes at least 3 arguments (3 given)
1120
1121
1122----------------------------------------------------------------------

You forgot to set S0. I think?

Finally, concerning the error FA which is different. Is the FA value very different?

@matthew-brett
@arokem
Owner

OK - we are down to just the one I meant. Strangely, while travis is getting:

ACTUAL: 1000.0
DESIRED: 686

I am getting on my machine:

ACTUAL: 979.0
DESIRED: 686

In both cases, the actual resulting FA is larger than expected from the calculation, based on the input axial and radial diffusivities.

@stefanv
Owner

Sorry, isclose was only added in NumPy 1.7.0.

@MrBago
Owner

@arokem, can you check to see if the signal being passed to TensorModel.fit if < 1. I believe that TensorModel sets all signals < min_signal to min_signal (so that we can take the log). You can try setting min_signal to something very small, like 1e-8 and see if that gives the right answer.

@arokem
Owner
@arokem
Owner
@Garyfallidis
Owner

@arokem Do you mean S0 to 100 rather than 0?

I am getting on my machine:

ACTUAL: 812.0
DESIRED: 686

Indeed the problem seems to take place in regions outside the fiber in the phantom. What I suppose should be the thresholded or masked region?

How do you set the mask using TensorModel? In multi_voxel.py the mask is a parameter at the model.fit(data, mask).
I don't see how this is available in TensorModel but I see it available in the old Tensor class.

Apart from that we need also to understand why the max is different from one computer to the other.
GGT!

@arokem arokem BF: Small eigenvalues => FA == 0
Make sure that when all eigenvalues of a tensor are small, they are effectively
set to 0, so that the FA doesn't misbehave with really small values of evals.

This also sets straight the tests for the orbital_phantom.
5d773c2
@arokem
Owner

Can we just merge this already?

@Garyfallidis
Owner

Yes.

@Garyfallidis Garyfallidis merged commit 6a69476 into from
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Commits on Nov 9, 2012
  1. @arokem

    TST: fvtk not needed here.

    arokem authored
  2. @arokem
  3. @arokem
  4. @arokem

    DOC: Small one.

    arokem authored
  5. @arokem
  6. @arokem

    DOC: Added reference.

    arokem authored
  7. @MrBago @arokem

    Using mean(signal)/std(noise) as definition of SNR

    MrBago authored arokem committed
  8. @arokem
  9. @arokem

    RF: This seems to work.

    arokem authored
    For some reason.
  10. @arokem
  11. @arokem
  12. @arokem
  13. @arokem
  14. @arokem
  15. @arokem

    DOC: Typos.

    arokem authored
  16. @arokem
  17. @arokem
Commits on Nov 10, 2012
  1. @arokem
  2. @arokem

    RF: Moved noise-addition functions to the voxel sub-module.

    arokem authored
    Also: adopted the convention used by Descoteaux et al. to define SNR.
Commits on Nov 19, 2012
  1. @arokem

    RF: More noise-related stuff.

    arokem authored
Commits on Nov 29, 2012
  1. @stefanv
Commits on Nov 30, 2012
  1. @arokem
  2. @arokem

    BF: Fixing failing tests.

    arokem authored
Commits on Dec 5, 2012
  1. @arokem

    BF: Small eigenvalues => FA == 0

    arokem authored
    Make sure that when all eigenvalues of a tensor are small, they are effectively
    set to 0, so that the FA doesn't misbehave with really small values of evals.
    
    This also sets straight the tests for the orbital_phantom.
This page is out of date. Refresh to see the latest.
View
1  dipy/core/gradients.py
@@ -72,6 +72,7 @@ def info(self):
print(' min %f ' % self.bvecs.min())
print(' max %f ' % self.bvecs.max())
+
def gradient_table_from_bvals_bvecs(bvals, bvecs, b0_threshold=0, atol=1e-2,
**kwargs):
"""Creates a GradientTable from a bvals array and a bvecs array
View
15 dipy/reconst/dti.py
@@ -65,6 +65,7 @@ def fit(self, data):
"""
dti_params = self.fit_method(self.design_matrix, data,
*self.args, **self.kwargs)
+
return TensorFit(self, dti_params)
@@ -140,7 +141,7 @@ def fa(self):
ev3 = evals[..., 2]
# Make sure not to get nans:
- all_zero = (ev1 == 0) & (ev2 == 0) & (ev3 == 0)
+ all_zero = np.allclose([ev1, ev2, ev3], 0)
fa = np.sqrt(0.5 * ((ev1 - ev2)**2 + (ev2 - ev3)**2 + (ev3 - ev1)**2)
/ (ev1*ev1 + ev2*ev2 + ev3*ev3 + all_zero))
@@ -273,9 +274,10 @@ def _wls_iter(ols_fit, design_matrix, sig, min_signal=1):
sig = np.maximum(sig, min_signal) #throw out zero signals
log_s = np.log(sig)
w = np.exp(np.dot(ols_fit, log_s))
- D = np.dot(np.linalg.pinv(design_matrix*w[:,None]), w*log_s)
+ D = np.dot(np.linalg.pinv(design_matrix * w[:,None]), w*log_s)
+ # D, _, _, _ = np.linalg.lstsq(design_matrix * w[:, None], log_s)
tensor = from_lower_triangular(D)
- return decompose_tensor(tensor)
+ return decompose_tensor(tensor, minimum_eval=np.finfo(float).eps)
def _ols_iter(inv_design, sig, min_signal=1):
@@ -286,7 +288,7 @@ def _ols_iter(inv_design, sig, min_signal=1):
log_s = np.log(sig)
D = np.dot(inv_design, log_s)
tensor = from_lower_triangular(D)
- return decompose_tensor(tensor)
+ return decompose_tensor(tensor, minimum_eval=np.finfo(float).eps)
def ols_fit_tensor(design_matrix, data, min_signal=1):
@@ -473,7 +475,7 @@ def tensor_eig_from_lo_tri(data):
return dti_params
-def decompose_tensor(tensor):
+def decompose_tensor(tensor, minimum_eval=0):
"""
Returns eigenvalues and eigenvectors given a diffusion tensor
@@ -499,6 +501,7 @@ def decompose_tensor(tensor):
See Also
--------
numpy.linalg.eig
+
"""
#outputs multiplicity as well so need to unique
@@ -510,7 +513,7 @@ def decompose_tensor(tensor):
eigenvals = eigenvals[order]
#Forcing negative eigenvalues to 0
- eigenvals = np.maximum(eigenvals, 0)
+ eigenvals = np.maximum(eigenvals, minimum_eval)
# b ~ 10^3 s/mm^2 and D ~ 10^-4 mm^2/s
# eigenvecs: each vector is columnar
View
17 dipy/reconst/tests/test_dsi.py
@@ -18,12 +18,11 @@ def test_dsi():
#load icosahedron sphere
sphere2 = create_unit_sphere(5)
btable = np.loadtxt(get_data('dsi515btable'))
- bvals = btable[:,0]
- bvecs = btable[:,1:]
- data, golden_directions = SticksAndBall(bvals, bvecs, d=0.0015,
+ gtab = gradient_table(btable[:,0], btable[:,1:])
+ data, golden_directions = SticksAndBall(gtab, d=0.0015,
S0=100, angles=[(0, 0), (90, 0)],
fractions=[50, 50], snr=None)
- gtab = gradient_table(bvals, bvecs)
+
ds = DiffusionSpectrumModel(gtab)
#symmetric724
ds.direction_finder.config(sphere=sphere, min_separation_angle=25,
@@ -60,22 +59,20 @@ def test_dsi():
def sticks_and_ball_dummies(gtab):
- bvals=gtab.bvals
- bvecs=gtab.bvecs
sb_dummies={}
- S, sticks = SticksAndBall(bvals, bvecs, d=0.0015, S0=100,
+ S, sticks = SticksAndBall(gtab, d=0.0015, S0=100,
angles=[(0, 0)],
fractions=[100], snr=None)
sb_dummies['1fiber'] = (S, sticks)
- S, sticks = SticksAndBall(bvals, bvecs, d=0.0015, S0=100,
+ S, sticks = SticksAndBall(gtab, d=0.0015, S0=100,
angles=[(0, 0), (90, 0)],
fractions=[50, 50], snr=None)
sb_dummies['2fiber'] = (S, sticks)
- S, sticks = SticksAndBall(bvals, bvecs, d=0.0015, S0=100,
+ S, sticks = SticksAndBall(gtab, d=0.0015, S0=100,
angles=[(0, 0), (90, 0), (90, 90)],
fractions=[33, 33, 33], snr=None)
sb_dummies['3fiber'] = (S, sticks)
- S, sticks = SticksAndBall(bvals, bvecs, d=0.0015, S0=100,
+ S, sticks = SticksAndBall(gtab, d=0.0015, S0=100,
angles=[(0, 0), (90, 0), (90, 90)],
fractions=[0, 0, 0], snr=None)
sb_dummies['isotropic'] = (S, sticks)
View
12 dipy/reconst/tests/test_dti.py
@@ -289,3 +289,15 @@ def test_from_lower_triangular():
D = D * np.ones((5, 4, 1))
tensor = from_lower_triangular(D)
assert_array_equal(tensor, result)
+
+def test_all_constant():
+ """
+
+ """
+ bvecs, bvals = read_bvec_file(get_data('55dir_grad.bvec'))
+ gtab = grad.gradient_table_from_bvals_bvecs(bvals, bvecs.T)
+ fit_methods = ['LS', 'OLS', 'NNLS']
+ for fit_method in fit_methods:
+ dm = dti.TensorModel(gtab, )
+ assert_almost_equal(dm.fit(np.zeros(bvals.shape[0])).fa, 0)
+ assert_almost_equal(dm.fit(100 * np.ones(bvals.shape[0])).fa, 0)
View
4 dipy/reconst/tests/test_gqi.py
@@ -22,10 +22,10 @@ def test_gqi():
btable = np.loadtxt(get_data('dsi515btable'))
bvals = btable[:,0]
bvecs = btable[:,1:]
- data, golden_directions = SticksAndBall(bvals, bvecs, d=0.0015,
+ gtab = gradient_table(bvals, bvecs)
+ data, golden_directions = SticksAndBall(gtab, d=0.0015,
S0=100, angles=[(0, 0), (90, 0)],
fractions=[50, 50], snr=None)
- gtab = gradient_table(bvals, bvecs)
gq = GeneralizedQSamplingModel(gtab, method='gqi2', sampling_length=1.4)
#symmetric724
gq.direction_finder.config(sphere=sphere, min_separation_angle=25,
View
4 dipy/reconst/tests/test_shm.py
@@ -142,8 +142,8 @@ def make_fake_signal():
evecs1 = evecs0
a = evecs0[0]
b = evecs1[1]
- S1 = single_tensor(gtab.bvals, gtab.bvecs, .55, evals[0], evecs0)
- S2 = single_tensor(gtab.bvals, gtab.bvecs, .45, evals[1], evecs1)
+ S1 = single_tensor(gtab, .55, evals[0], evecs0)
+ S2 = single_tensor(gtab, .45, evals[1], evecs1)
return S1 + S2, gtab, np.vstack([a, b])
View
343 dipy/sims/phantom.py
@@ -1,12 +1,69 @@
import numpy as np
-from dipy.sims.voxel import SingleTensor
+import scipy.stats as stats
+
+from dipy.sims.voxel import SingleTensor, diffusion_evals
+import dipy.sims.voxel as vox
from dipy.core.geometry import vec2vec_rotmat
-from dipy.data import get_data
-#from dipy.viz import fvtk
-#from dipy.reconst.dti import Tensor
+from dipy.data import get_data
+from dipy.core.gradients import gradient_table
+
+
+def add_noise(vol, snr=1.0, S0=None, noise_type='rician'):
+ """ Add noise of specified distribution to a 4D array.
+
+ Parameters
+ -----------
+ vol : array, shape (X,Y,Z,W)
+ Diffusion measurements in `W` directions at each ``(X, Y, Z)`` voxel
+ position.
+ snr : float, optional
+ The desired signal-to-noise ratio. (See notes below.)
+ S0 : float, optional
+ Reference signal for specifying `snr` (defaults to 1).
+ noise_type : string, optional
+ The distribution of noise added. Can be either 'gaussian' for Gaussian
+ distributed noise, 'rician' for Rice-distributed noise (default) or
+ 'rayleigh' for a Rayleigh distribution.
+
+ Returns
+ --------
+ vol : array, same shape as vol
+ Volume with added noise.
+
+ Notes
+ -----
+ SNR is defined here, following [1]_, as ``S0 / sigma``, where ``sigma`` is
+ the standard deviation of the two Gaussian distributions forming the real
+ and imaginary components of the Rician noise distribution (see [2]_).
+
+ References
+ ----------
+ .. [1] Descoteaux, Angelino, Fitzgibbons and Deriche (2007) Regularized,
+ fast and robust q-ball imaging. MRM, 58: 497-510
+ .. [2] Gudbjartson and Patz (2008). The Rician distribution of noisy MRI
+ data. MRM 34: 910-914.
+
+ Examples
+ --------
+ >>> signal = np.arange(800).reshape(2, 2, 2, 100)
+ >>> signal_w_noise = add_noise(signal, snr=10, noise_type='rician')
+
+ """
+ orig_shape = vol.shape
+ vol_flat = np.reshape(vol.copy(), (-1, vol.shape[-1]))
+
+ if S0 is None:
+ S0 = np.max(vol)
+
+ for vox_idx, signal in enumerate(vol_flat):
+ vol_flat[vox_idx] = vox.add_noise(signal, snr=snr, S0=S0,
+ noise_type=noise_type)
+
+ return np.reshape(vol_flat, orig_shape)
+
def diff2eigenvectors(dx,dy,dz):
- """ numerical derivatives 2 eigenvectors
+ """ numerical derivatives 2 eigenvectors
"""
basis=np.eye(3)
u=np.array([dx,dy,dz])
@@ -18,213 +75,151 @@ def diff2eigenvectors(dx,dy,dz):
eigs=np.zeros((3,3))
eigs[:,0]=eig0
eigs[:,1]=eig1
- eigs[:,2]=eig2
+ eigs[:,2]=eig2
return eigs, R
-def orbital_phantom(bvals=None,
- bvecs=None,
- evals=np.array([1.4,.35,.35])*10**(-3),
- func=None,
- t=np.linspace(0,2*np.pi,1000),
- datashape=(64,64,64,65),
- origin=(32,32,32),
- scale=(25,25,25),
- angles=np.linspace(0,2*np.pi,32),
- radii=np.linspace(0.2,2,6),
- S0=100.):
- """ Create a phantom based on a 3d orbit f(t)->(x,y,z)
-
+def orbital_phantom(gtab=None,
+ evals=diffusion_evals,
+ func=None,
+ t=np.linspace(0, 2 * np.pi, 1000),
+ datashape=(64, 64, 64, 65),
+ origin=(32, 32, 32),
+ scale=(25, 25, 25),
+ angles=np.linspace(0, 2 * np.pi, 32),
+ radii=np.linspace(0.2, 2, 6),
+ S0=100.,
+ snr=None):
+ """Create a phantom based on a 3-D orbit ``f(t) -> (x,y,z)``.
+
Parameters
-----------
- bvals : array, shape (N,)
- bvecs : array, shape (N,3)
+ gtab : GradientTable
+ Gradient table of measurement directions.
evals : array, shape (3,)
- tensor eigenvalues
- func : user defined function f(t)->(x,y,z)
- It could be desirable for -1=<x,y,z <=1
- if None creates a circular orbit
+ Tensor eigenvalues.
+ func : user defined function f(t)->(x,y,z)
+ It could be desirable for ``-1=<x,y,z <=1``.
+ If None creates a circular orbit.
t : array, shape (K,)
- represents time for the orbit
- Default is np.linspace(0,2*np.pi,1000)
+ Represents time for the orbit. Default is
+ ``np.linspace(0, 2 * np.pi, 1000)``.
datashape : array, shape (X,Y,Z,W)
- size of the output simulated data
+ Size of the output simulated data
origin : tuple, shape (3,)
- define the center for the volume
+ Define the center for the volume
scale : tuple, shape (3,)
- scale the function before applying to the grid
+ Scale the function before applying to the grid
angles : array, shape (L,)
- density angle points, always perpendicular to the first eigen vector
- Default np.linspace(0,2*np.pi,32),
+ Density angle points, always perpendicular to the first eigen vector
+ Default np.linspace(0, 2 * np.pi, 32).
radii : array, shape (M,)
- thickness radii
- Default np.linspace(0.2,2,6)
- angles and radii define the total thickness options
- S0 : double, simulated signal without diffusion gradients applied
- Default 100.
- snr : signal to noise ratio
- Used for applying rician noise to the data.
- Default 200. Common is 20.
- background_noise : boolean, Default False
-
+ Thickness radii. Default ``np.linspace(0.2, 2, 6)``.
+ angles and radii define the total thickness options
+ S0 : double, optional
+ Maximum simulated signal. Default 100.
+ snr : float, optional
+ The signal to noise ratio set to apply Rician noise to the data.
+ Default is to not add noise at all.
+
Returns
- ---------
+ -------
data : array, shape (datashape)
-
- Notes
+
+ See Also
--------
- Crossings can be created by adding multiple orbitual_phantom outputs.
-
+ add_noise
+
Examples
---------
-
- def f(t):
- x=np.sin(t)
- y=np.cos(t)
- z=np.linspace(-1,1,len(x))
- return x,y,z
-
- data=orbitual_phantom(func=f)
-
+
+ >>> def f(t):
+ ... x = np.sin(t)
+ ... y = np.cos(t)
+ ... z = np.linspace(-1, 1, len(x))
+ ... return x, y, z
+
+ >>> data = orbital_phantom(func=f)
+
"""
-
- if bvals==None:
- fimg,fbvals,fbvecs=get_data('small_64D')
- bvals=np.load(fbvals)
- bvecs=np.load(fbvecs)
- bvecs[np.isnan(bvecs)]=0
-
- if func==None:
- x=np.sin(t)
- y=np.cos(t)
- z=np.zeros(t.shape)
+
+ if gtab is None:
+ fimg, fbvals, fbvecs = get_data('small_64D')
+ gtab = gradient_table(fbvals, fbvecs)
+
+ if func is None:
+ x = np.sin(t)
+ y = np.cos(t)
+ z = np.zeros(t.shape)
else:
- x,y,z=func(t)
-
- #stop
-
- dx=np.diff(x)
- dy=np.diff(y)
- dz=np.diff(z)
-
- x=scale[0]*x+origin[0]
- y=scale[1]*y+origin[1]
- z=scale[2]*z+origin[2]
-
- bx=np.zeros(len(angles))
- by=np.sin(angles)
- bz=np.cos(angles)
-
- vol=np.zeros(datashape)
-
+ x, y, z = func(t)
+
+ dx = np.diff(x)
+ dy = np.diff(y)
+ dz = np.diff(z)
+
+ x = scale[0] * x + origin[0]
+ y = scale[1] * y + origin[1]
+ z = scale[2] * z + origin[2]
+
+ bx = np.zeros(len(angles))
+ by = np.sin(angles)
+ bz = np.cos(angles)
+
+ # The entire volume is considered to be inside the brain.
+ # Voxels without a fiber crossing through them are taken
+ # to be isotropic with signal = S0.
+ vol = np.zeros(datashape) + S0
+
for i in range(len(dx)):
- evecs,R=diff2eigenvectors(dx[i],dy[i],dz[i])
- S=SingleTensor(bvals,bvecs,S0,evals,evecs,snr=None)
- #print sigma, S0/snr, S0, snr
- vol[x[i],y[i],z[i],:]+=S
+ evecs, R = diff2eigenvectors(dx[i], dy[i], dz[i])
+ S = SingleTensor(gtab, S0, evals, evecs, snr=None)
+
+ vol[x[i], y[i], z[i], :] += S
+
for r in radii:
for j in range(len(angles)):
- rb=np.dot(R,np.array([bx[j],by[j],bz[j]]))
- vol[x[i]+r*rb[0],y[i]+r*rb[1],z[i]+r*rb[2]]+=S
-
- #ten=Tensor(vol,bvals,bvecs)
- #FA=ten.fa()
- #FA[np.isnan(FA)]=0
- #vol[np.isnan(vol)]=0
- return vol
+ rb = np.dot(R,np.array([bx[j], by[j], bz[j]]))
-def add_rician_noise(vol,snr=20):
- """ add rician noise in 4D diffusion data
-
- Parameters
- -----------
- vol : array, shape (X,Y,Z,W)
- snr : float,
- signal to noise ratio
-
- Returns
- --------
- voln : array, same shape as vol
- vol with additional rician noise
-
- Reference
- ----------
- http://my.fit.edu/~kozaitis/Matlab/MatNoiseIm.html
- """
-
- if len(vol.shape)==4:
- def gaussian_noise(vol):
- voln=np.random.randn(*vol.shape[:3])
- pvol=np.sum(vol[...,0]**2) #power of initial volume
- pnoise=np.sum(np.random.randn(*voln.shape[:3])**2) #power of noise volume
- K=pvol/pnoise
- #print pvol,pnoise,K
- return np.sqrt(K/np.float(snr))*np.random.randn(*vol.shape)
- noise1=gaussian_noise(vol)
- noise2=gaussian_noise(vol)
- return np.sqrt((vol+noise1)**2+noise2**2)
-
-def add_gaussian_noise(vol,snr=20):
- """ add gaussian noise in a 4D array with a specific snr
-
- Parameters
- -----------
- vol : array, shape (X,Y,Z,W)
- snr : float,
- signal to noise ratio
-
- Returns
- --------
- voln : array, same shape as vol
- vol with additional rician noise
-
- Reference
- ----------
- http://my.fit.edu/~kozaitis/Matlab/MatNoiseIm.html
- """
-
- if len(vol.shape)==4:
- voln=np.random.randn(*vol.shape[:3])
- pvol=np.sum(vol[...,0]**2) #power of initial volume
- pnoise=np.sum((np.random.randn(*voln.shape)-.5)**2) #power of noise volume
- K=pvol/pnoise
- print pvol,pnoise,K
- noise=np.sqrt(K/np.float(snr))*(np.random.randn(*vol.shape))
- return vol+noise
-
-
-
-
-
+ vol[x[i] + r * rb[0],
+ y[i] + r * rb[1],
+ z[i] + r * rb[2]] += S
+ vol = vol / np.max(vol, axis=-1)[..., np.newaxis]
+ vol *= S0
+
+ if snr is not None:
+ vol = add_noise(vol, snr, S0=S0, noise_type='rician')
+
+ return vol
if __name__ == "__main__":
-
- ##TODO: this can become a nice tutorial for generating phantoms
-
+
+ ## TODO: this can become a nice tutorial for generating phantoms
+
def f(t):
x=np.sin(t)
y=np.cos(t)
#z=np.zeros(t.shape)
z=np.linspace(-1,1,len(x))
return x,y,z
-
+
#helix
- vol=orbital_phantom(func=f)
-
+ vol=orbital_phantom(func=f)
+
def f2(t):
x=np.linspace(-1,1,len(t))
- y=np.linspace(-1,1,len(t))
+ y=np.linspace(-1,1,len(t))
z=np.zeros(x.shape)
return x,y,z
#first direction
vol2=orbital_phantom(func=f2)
-
+
def f3(t):
x=np.linspace(-1,1,len(t))
- y=-np.linspace(-1,1,len(t))
+ y=-np.linspace(-1,1,len(t))
z=np.zeros(x.shape)
return x,y,z
@@ -232,22 +227,22 @@ def f3(t):
vol3=orbital_phantom(func=f3)
#double crossing
vol23=vol2+vol3
-
+
#"""
def f4(t):
x=np.zeros(t.shape)
y=np.zeros(t.shape)
z=np.linspace(-1,1,len(t))
return x,y,z
-
+
#triple crossing
vol4=orbital_phantom(func=f4)
vol234=vol23+vol4
-
+
voln=add_rician_noise(vol234)
-
+
#"""
-
+
#r=fvtk.ren()
#fvtk.add(r,fvtk.volume(vol234[...,0]))
#fvtk.show(r)
View
134 dipy/sims/tests/test_phantom.py
@@ -1,56 +1,90 @@
+from __future__ import division
+
import numpy as np
import nose
import nibabel as nib
-from nose.tools import assert_true, assert_false, assert_equal, assert_almost_equal
-from numpy.testing import assert_array_equal, assert_array_almost_equal
+import numpy.testing.decorators as dec
+
+from numpy.testing import (assert_, assert_equal, assert_array_equal,
+ assert_array_almost_equal, assert_almost_equal,
+ run_module_suite)
+
from dipy.core.geometry import vec2vec_rotmat
-from dipy.data import get_data
-from dipy.viz import fvtk
-from dipy.reconst.dti import Tensor
-from dipy.sims.phantom import orbital_phantom
+from dipy.data import get_data
+from dipy.reconst.dti import TensorModel
+from dipy.sims.phantom import orbital_phantom, add_noise
+from dipy.sims.voxel import single_tensor
+from dipy.core.gradients import gradient_table
+
+
+fimg,fbvals,fbvecs=get_data('small_64D')
+bvals=np.load(fbvals)
+bvecs=np.load(fbvecs)
+bvecs[np.isnan(bvecs)]=0
+
+gtab = gradient_table(bvals, bvecs)
+
+
+def f(t):
+ """
+ Helper function used to define a mapping time => xyz
+ """
+ x = np.linspace(-1,1,len(t))
+ y = np.linspace(-1,1,len(t))
+ z = np.linspace(-1,1,len(t))
+ return x,y,z
+
def test_phantom():
-
- def f(t):
- x=np.sin(t)
- y=np.cos(t)
- z=np.linspace(-1,1,len(x))
- return x,y,z
-
- fimg,fbvals,fbvecs=get_data('small_64D')
- bvals=np.load(fbvals)
- bvecs=np.load(fbvecs)
- bvecs[np.isnan(bvecs)]=0
-
- N=50 #timepoints
-
- vol=orbital_phantom(bvals=bvals,
- bvecs=bvecs,
- func=f,
- t=np.linspace(0,2*np.pi,N),
- datashape=(10,10,10,len(bvals)),
- origin=(5,5,5),
- scale=(3,3,3),
- angles=np.linspace(0,2*np.pi,16),
- radii=np.linspace(0.2,2,6))
-
- ten=Tensor(vol,bvals,bvecs)
- FA=ten.fa()
- FA[np.isnan(FA)]=0
-
- assert_equal(np.round(FA.max()*1000),707)
-
-
-
-
-
-if __name__ == "__main__":
- test_phantom()
-
-
-
-
-
-
-
-
+ N = 50
+
+ vol = orbital_phantom(gtab,
+ func=f,
+ t=np.linspace(0, 2 * np.pi, N),
+ datashape=(10, 10, 10, len(bvals)),
+ origin=(5, 5, 5),
+ scale=(3, 3, 3),
+ angles=np.linspace(0, 2 * np.pi, 16),
+ radii=np.linspace(0.2, 2, 6),
+ S0=100)
+
+ m = TensorModel(gtab)
+ t = m.fit(vol)
+ FA = t.fa
+ # print vol
+ FA[np.isnan(FA)] = 0
+ # 686 -> expected FA given diffusivities of [1500, 400, 400]
+ l1, l2, l3 = 1500e-6, 400e-6, 400e-6
+ expected_fa = (np.sqrt(0.5) * np.sqrt((l1 - l2)**2 + (l2-l3)**2 + (l3-l1)**2 )/np.sqrt(l1**2 + l2**2 + l3**2))
+
+ assert_array_almost_equal(FA.max(), expected_fa, decimal=2)
+
+
+def test_add_noise():
+ np.random.seed(1980)
+
+ N = 50
+ S0 = 100
+
+ options = dict(func=f,
+ t=np.linspace(0, 2 * np.pi, N),
+ datashape=(10, 10, 10, len(bvals)),
+ origin=(5, 5, 5),
+ scale=(3, 3, 3),
+ angles=np.linspace(0, 2 * np.pi, 16),
+ radii=np.linspace(0.2, 2, 6),
+ S0=S0)
+
+ vol = orbital_phantom(gtab, **options)
+
+ for snr in [10, 20, 30, 50]:
+ vol_noise = orbital_phantom(gtab, snr=snr, **options)
+
+ sigma = S0 / snr
+
+ assert_(np.abs(np.var(vol_noise - vol) - sigma ** 2) < 1)
+
+
+
+if __name__ == "__main__":
+ run_module_suite()
View
83 dipy/sims/tests/test_voxel.py
@@ -1,12 +1,22 @@
import numpy as np
-import nose
-import nibabel as nib
-from nose.tools import assert_true, assert_false, assert_equal, assert_almost_equal
-from numpy.testing import assert_array_equal, assert_array_almost_equal
-from dipy.sims.voxel import SingleTensor, multi_tensor_odf, all_tensor_evecs
+
+from nose.tools import (assert_true, assert_false, assert_equal,
+ assert_almost_equal)
+from numpy.testing import (assert_array_equal, assert_array_almost_equal,
+ assert_)
+
+from dipy.sims.voxel import (SingleTensor, multi_tensor_odf, all_tensor_evecs,
+ add_noise, single_tensor, sticks_and_ball)
from dipy.core.geometry import vec2vec_rotmat
from dipy.data import get_data, get_sphere
-from dipy.viz import fvtk
+from dipy.core.gradients import gradient_table
+
+
+fimg, fbvals, fbvecs = get_data('small_64D')
+bvals = np.load(fbvals)
+bvecs = np.load(fbvecs)
+gtab = gradient_table(bvals, bvecs)
+
def diff2eigenvectors(dx,dy,dz):
""" numerical derivatives 2 eigenvectors
@@ -24,39 +34,58 @@ def diff2eigenvectors(dx,dy,dz):
return eigs, R
+def test_sticks_and_ball():
+ d = 0.0015
+ S, sticks = sticks_and_ball(gtab, d=d, S0=1, angles=[(0,0),],
+ fractions=[100], snr=None)
+ assert_array_equal(sticks, [[0, 0, 1]])
+ S_st = SingleTensor(gtab, 1, evals=[d, 0, 0], evecs=[[0, 0, 0],
+ [0, 0, 0],
+ [1, 0, 0]])
+ assert_array_almost_equal(S, S_st)
+
+
def test_single_tensor():
+ evals = np.array([1.4, .35, .35]) * 10**(-3)
+ evecs = np.eye(3)
+ S = SingleTensor(gtab, 100, evals, evecs, snr=None)
+ assert_array_almost_equal(S[gtab.b0s_mask], 100)
+ assert_(np.mean(S[~gtab.b0s_mask]) < 100)
- fimg,fbvals,fbvecs=get_data('small_64D')
- bvals=np.load(fbvals)
- bvecs=np.load(fbvecs)
- #bvals=np.loadtxt(fbvals)
- #bvecs=np.loadtxt(fbvecs).T
- img=nib.load(fimg)
- data=img.get_data()
-
- evals=np.array([1.4,.35,.35])*10**(-3)
- evecs=np.eye(3)
- S=SingleTensor(bvals,bvecs,100,evals,evecs,snr=None)
- """
- colours=fvtk.colors(S,'jet')
- r=fvtk.ren()
- fvtk.add(r,fvtk.point(bvecs,colours))
- fvtk.show(r)
- """
+ from dipy.reconst.dti import TensorModel
+ m = TensorModel(gtab)
+ t = m.fit(S)
+
+ assert_array_almost_equal(t.fa, 0.707, decimal=3)
def test_multi_tensor():
sphere = get_sphere('symmetric724')
- vertices, faces = sphere.vertices, sphere.faces
+ vertices = sphere.vertices
mevals=np.array(([0.0015, 0.0003, 0.0003],
[0.0015, 0.0003, 0.0003]))
e0 = np.array([1, 0, 0.])
e1 = np.array([0., 1, 0])
mevecs=[all_tensor_evecs(e0), all_tensor_evecs(e1)]
- odf = multi_tensor_odf(vertices, [0.5,0.5], mevals, mevecs)
+ odf = multi_tensor_odf(vertices, [0.5, 0.5], mevals, mevecs)
+
+ assert_(odf.shape == (len(vertices),))
+ assert_(np.all(odf <= 1) & np.all(odf >= 0))
+
+
+def test_snr():
+ np.random.seed(1978)
+
+ s = single_tensor(gtab)
+
+ # For reasonably large SNR, var(signal) ~= sigma**2, where sigma = 1/SNR
+ for snr in [5, 10, 20]:
+ sigma = 1.0 / snr
+ for j in range(1000):
+ s_noise = add_noise(s, snr, 1, noise_type='rician')
+
+ assert_array_almost_equal(np.var(s_noise - s), sigma**2, decimal=2)
- assert odf.shape == (len(vertices),)
- assert np.all(odf <= 1) & np.all(odf >= 0)
if __name__ == "__main__":
View
148 dipy/sims/voxel.py
@@ -15,16 +15,104 @@
diffusion_evals = np.array([1500e-6, 400e-6, 400e-6])
-def sticks_and_ball(bvals, gradients, d=0.0015, S0=100, angles=[(0,0), (90,0)],
+def _add_gaussian(sig, noise1, noise2):
+ """
+ Helper function to add_noise
+
+ This one simply adds one of the Gaussians to the sig and ignores the other
+ one.
+ """
+ return sig + noise1
+
+
+def _add_rician(sig, noise1, noise2):
+ """
+ Helper function to add_noise.
+
+ This does the same as abs(sig + complex(noise1, noise2))
+
+ """
+ return np.sqrt((sig + noise1)**2 + noise2**2)
+
+
+def _add_rayleigh(sig, noise1, noise2):
+ """
+ Helper function to add_noise
+
+ The Rayleigh distribution is $\sqrt\{Gauss_1^2 + Gauss_2^2}$.
+
+ """
+ return sig + np.sqrt(noise1**2 + noise2**2)
+
+
+def add_noise(signal, snr, S0, noise_type='rician'):
+ r""" Add noise of specified distribution to the signal from a single voxel.
+
+ Parameters
+ -----------
+ signal : 1-d ndarray
+ The signal in the voxel.
+ snr : float
+ The desired signal-to-noise ratio. (See notes below.)
+ If `snr` is None, return the signal as-is.
+ S0 : float
+ Reference signal for specifying `snr`.
+ noise_type : string, optional
+ The distribution of noise added. Can be either 'gaussian' for Gaussian
+ distributed noise, 'rician' for Rice-distributed noise (default) or
+ 'rayleigh' for a Rayleigh distribution.
+
+ Returns
+ --------
+ signal : array, same shape as the input
+ Signal with added noise.
+
+ Notes
+ -----
+ SNR is defined here, following [1]_, as ``S0 / sigma``, where ``sigma`` is
+ the standard deviation of the two Gaussian distributions forming the real
+ and imaginary components of the Rician noise distribution (see [2]_).
+
+ References
+ ----------
+ .. [1] Descoteaux, Angelino, Fitzgibbons and Deriche (2007) Regularized,
+ fast and robust q-ball imaging. MRM, 58: 497-510
+ .. [2] Gudbjartson and Patz (2008). The Rician distribution of noisy MRI
+ data. MRM 34: 910-914.
+
+ Examples
+ --------
+ >>> signal = np.arange(800).reshape(2, 2, 2, 100)
+ >>> signal_w_noise = add_noise(signal, 10., 100., noise_type='rician')
+
+ """
+ if snr is None:
+ return signal
+
+ sigma = S0 / snr
+
+ noise_adder = {'gaussian': _add_gaussian,
+ 'rician': _add_rician,
+ 'rayleigh': _add_rayleigh}
+
+ noise1 = np.random.normal(0, sigma, size=signal.shape)
+
+ if noise_type == 'gaussian':
+ noise2 = None
+ else:
+ noise2 = np.random.normal(0, sigma, size=signal.shape)
+
+ return noise_adder[noise_type](signal, noise1, noise2)
+
+
+def sticks_and_ball(gtab, d=0.0015, S0=100, angles=[(0,0), (90,0)],
fractions=[35,35], snr=20):
""" Simulate the signal for a Sticks & Ball model.
Parameters
-----------
- bvals : (N,) ndarray
- B-values for measurements.
- gradients : (N,3) ndarray
- Also known as b-vectors.
+ gtab : GradientTable
+ Signal measurement directions.
d : float
Diffusivity value.
S0 : float
@@ -33,9 +121,10 @@ def sticks_and_ball(bvals, gradients, d=0.0015, S0=100, angles=[(0,0), (90,0)],
List of K polar angles (in degrees) for the sticks or array of M
sticks as Cartesian unit vectors.
fractions : float
- Percentage of each stick.
+ Percentage of each stick. Remainder to 100 specifies isotropic
+ component.
snr : float
- Signal to noise ratio, assuming gaussian noise. If set to None, no
+ Signal to noise ratio, assuming Rician noise. If set to None, no
noise is added.
Returns
@@ -52,10 +141,9 @@ def sticks_and_ball(bvals, gradients, d=0.0015, S0=100, angles=[(0,0), (90,0)],
Neuroimage, 2007.
"""
-
fractions = [f / 100. for f in fractions]
f0 = 1 - np.sum(fractions)
- S = np.zeros(len(gradients))
+ S = np.zeros(len(gtab.bvals))
angles=np.array(angles)
if angles.shape[-1] == 3:
@@ -65,37 +153,28 @@ def sticks_and_ball(bvals, gradients, d=0.0015, S0=100, angles=[(0,0), (90,0)],
for pair in angles]
sticks = np.array(sticks)
- for (i, g) in enumerate(gradients[1:]):
- S[i + 1] = f0 * np.exp(-bvals[i+1] * d) + \
+ for (i, g) in enumerate(gtab.bvecs[1:]):
+ S[i + 1] = f0 * np.exp(-gtab.bvals[i+1] * d) + \
np.sum([
- fractions[j] * np.exp(-bvals[i + 1] * d * np.dot(s, g)**2)
+ fractions[j] * np.exp(-gtab.bvals[i + 1] * d * np.dot(s, g)**2)
for (j,s) in enumerate(sticks)
])
S[i + 1] = S0 * S[i + 1]
- S[0] = S0
- if snr is not None:
- std = S0 / snr
- S = S + np.random.randn(len(S)) * std
+ S[gtab.b0s_mask] = S0
+ S = add_noise(S, snr, S0)
return S, sticks
-def single_tensor(bvals, gradients, S0=1, evals=None, evecs=None, snr=None):
+def single_tensor(gtab, S0=1, evals=None, evecs=None, snr=None):
""" Simulated Q-space signal with a single tensor.
Parameters
-----------
- bvals : (N,) array
- B-values for measurements. The b-value is also ``b = \tau |q|^2``,
- where ``\tau`` is the time allowed for attenuation and ``q`` is the
- measurement position vector in Q-space (signal-space or Fourier-space).
- If b is too low, there is not enough attenuation to measure. With b
- too high, the signal to noise ratio increases.
- gradients : (N, 3) or (M, N, 3) ndarray
- Measurement gradients / directions, also known as b-vectors, as 3D unit
- vectors (either in a list or on a grid).
+ gtab : GradientTable
+ Measurement directions.
S0 : double,
Strength of signal in the presence of no diffusion gradient (also
called the ``b=0`` value).
@@ -106,7 +185,7 @@ def single_tensor(bvals, gradients, S0=1, evals=None, evecs=None, snr=None):
Eigenvectors of the tensor. You can also think of this as a rotation
matrix that transforms the direction of the tensor.
snr : float
- Signal to noise ratio, assuming gaussian noise. None implies no noise.
+ Signal to noise ratio, assuming Rician noise. None implies no noise.
Returns
--------
@@ -129,24 +208,17 @@ def single_tensor(bvals, gradients, S0=1, evals=None, evecs=None, snr=None):
if evecs is None:
evecs = np.eye(3)
- out_shape = gradients.shape[:gradients.ndim - 1]
+ out_shape = gtab.bvecs.shape[:gtab.bvecs.ndim - 1]
+ gradients = gtab.bvecs.reshape(-1, 3)
- gradients = gradients.reshape(-1, 3)
R = np.asarray(evecs)
S = np.zeros(len(gradients))
D = R.dot(np.diag(evals)).dot(R.T)
for (i, g) in enumerate(gradients):
- S[i] = S0 * np.exp(-bvals[i] * g.T.dot(D).dot(g))
-
- """ Alternative suggestion which works with multiple b0s
- design = design_matrix(bval, gradients.T)
- S = np.exp(np.dot(design, lower_triangular(D)))
- """
+ S[i] = S0 * np.exp(-gtab.bvals[i] * g.T.dot(D).dot(g))
- if snr is not None:
- std = S0 / snr
- S = S + np.random.randn(len(S)) * std
+ S = add_noise(S, snr, S0)
return S.reshape(out_shape)
Something went wrong with that request. Please try again.