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

DKI project: PR#1 Simulations to test DKI #582

Merged
merged 26 commits into from Jun 8, 2015

Conversation

Projects
None yet
4 participants
@RafaelNH
Contributor

RafaelNH commented Mar 1, 2015

I add some functions on dipy.sims.voxel module for DKI simulations (multi_tensor_dki, compute_Wijkl, single_diffkurt_tensors, dki_design_matrix). This PR was discussed to be the first milestone of the DKI project. To see how to use the functions please give a look to the added script on doc/examples/simulations_dki.py. Basically, the objective of the new functions is to produce ground truth DT and KT tensors to test diffusion fits (you can simulate tensors for any single voxel crossing fiber configuration ;) ). DW signals are also simulated from the computed tensors using the DKI model. Future steps on the DKI project can use these simulations to check if the DKI implementation is correct. DKI simulations modules are also relevant for DKI fit comparison studies.

Implementation was done according to PEP 8 standards and the nomenclature is consistent to previous dipy functions. However, I notice that the order of DT elements is different from previous dipy modules (as dipy.reconst.dti.py), I decided to keep this order since I notice that they are consistent to the in progress DKI implementations. Perhaps is better to transform all new functions to dipy assumed format - Dxx Dxy Dyy Dxz Dyz Dzz. What do you thing? Please let me know if you need help doing this on the submitted changes after its revision.

Enjoy,
R Neto Henriques

# kurtosis tensor - Impact on the development of robust tractography
# procedures and novel biomarkers", NeuroImage, 2015; 111:85-99.
#
diffusion_tensor = np.array([1638e-6, 444e-6, 444e-6, 0, 0, 0])

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

Would it be OK to use the same diffusion evals as the ones just above (on L16)? It's starting to get crowded in this namespace...

This comment has been minimized.

@RafaelNH

RafaelNH Mar 2, 2015

Contributor

Can we instead change the evals? The parameters diffusion_tensor and kurtosis_tensor are the defaults for a voxel of well aligned fibers and assuming the multi-compartmental model's parameters on Table 1 of Neto Henriques et al. (2015), so I want to keep then consistent to each other. Another option is to remove these global variables at all, and require users always giving the DT and KT as input of function single_diffkurt_tensors.

diffusion_tensor = np.array([1638e-6, 444e-6, 444e-6, 0, 0, 0])
kurtosis_tensor = np.array([1.7068, 0.8010, 0.8010, 0, 0, 0, 0, 0, 0, 0.3897,
0.3897, 0.2670, 0, 0, 0])

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

Just for my own understanding: are the kurtosis elements of similar units to the DT evals? Why are they 3 orders of magnitude different?

This comment has been minimized.

@RafaelNH

RafaelNH Mar 2, 2015

Contributor

Kurtosis don't have the same units of DT evals! Actually kurtosis is dimensionless - it is an statistical measure quantifying deviations of the water diffusion propagation function from gaussian. Typical values of kurtosis are on the range of [0 3].

def multi_tensor_dki(gtab, mevals, S0=100, angles=[(0., 0.), (90., 0.)],
fractions=[50, 50], snr=20):
r"""Simulate the elements of the diffusion and diffusion kurtosis tensor as

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

Possible typo, Should this be: "...the diffusion tensor and diffusion kurtosis tensor..."?

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

Would be nice to have a one-line summary at the top, and then a more expansive explanation.

This comment has been minimized.

@RafaelNH

RafaelNH Mar 2, 2015

Contributor

Yes, you are right it should be diffusion tensor

r"""Simulate the elements of the diffusion and diffusion kurtosis tensor as
well as the diffusion signal S based on the DKI model (i.e. taylor
accumulante larger than the fourth order are ignored). Simulations are

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

"Taylor expansion"?

This comment has been minimized.

@RafaelNH

RafaelNH Mar 2, 2015

Contributor

Perhaps is better to say 'Taylor expansion coefficients larger than ...'.

r"""Simulate the elements of the diffusion and diffusion kurtosis tensor as
well as the diffusion signal S based on the DKI model (i.e. taylor
accumulante larger than the fourth order are ignored). Simulations are
based on multicompartmental models which assumes that tissue are described

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

"tissue is described" or "tissue compartments are described"

This comment has been minimized.

@RafaelNH

RafaelNH Mar 2, 2015

Contributor

tissue is described

gtab : GradientTable
mevals : array (K, 3)
eigenvalues of the diffusion tensor for each individual compartment
S0 : float

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

Mark optional arguments as "(optional)" and tell us what the default values are.

This comment has been minimized.

@RafaelNH

RafaelNH Mar 2, 2015

Contributor

ok, will do

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

This is still needed

sticks = np.array(sticks)
# computing a 3D matrice containing the individual DT components
mD=np.zeros((len(fractions),3,3))

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

PEP8: spaces around operators and after comma:

mD = np.zeros((len(fractions, 3, 3))
mD[i] = dot(dot(R, np.diag(mevals[i])), R.T)
# compute DT
D=np.zeros((3,3))

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

PEP8

# compute DT
D=np.zeros((3,3))
for i in range(len(fractions)):
D=D+fractions[i]*mD[i]

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

PEP8

return dt, kt, S
def compute_Wijkl(Dc,frac,i,j,k,el,DT=None):
r""" # function that computes the elements of KT by especifying their

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

You can remove the '#' at the beginning. Also: "especifying" => "specifying"

def compute_Wijkl(Dc,frac,i,j,k,el,DT=None):
r""" # function that computes the elements of KT by especifying their
four indexes i, j, k and el. Each index can be equal to x, y or z which are
codified respectively by 0, 1 and 2. Simulations are based on

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

I am not sure I understand the sentence: "x, y, or z which are codified by 0, 1, and 2". Does this have to do with the principal orientation of the kurtosis tensor in 3D space?

This comment has been minimized.

@RafaelNH

RafaelNH Mar 2, 2015

Contributor

Unfortunately, defining a principal orientation of KT is not trivial (and what will this principal orientation mean on KT? - don't forget that KT is not an ellipsoid). KT have 15 elements: Wxxxx Wyyyy Wzzzz Wxxxy Wxxxz Wxyyy Wyyyz Wxzzz Wyzzz Wxxyy Wxxzz Wyyzz Wxxyz Wxyyz Wxyzz - what I wanted to say for codifying each element I use 0 for x, 1 for y and 2 for z.

wijkl=(wijkl-
DT[i][j]*DT[k][el]-DT[i][k]*DT[j][el]-DT[i][el]*DT[j][k])/(MD**2)
return wijkl

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

PEP8 in this entire function: space around operators and around =. Make sure to break lines into 80 characters when you add those spaces.

Returns
--------
S : (N,) ndarray
Simulated signal based on the DKI model S=S0*exp(-b*D+1/6*b**2*D**2*K)

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

You might want to render the equation with the latex engine, by adding $ signs around it, so:

Simulated signal based on the DKI model: $S=S0e^{-b D+\frac{1}{6} b^2 D^2 K}$
Elements of the diffusion kurtosis tensor.
snr : float
Signal to noise ratio, assuming Rician noise. None implies no noise.

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

Again: optional arguments and their default values.

A=dki_design_matrix(gtab)
# define vector of DKI parameters
MD=sum(dt[0:3])/3

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

Is this also MD when the last 3 elements are non-zero?

This comment has been minimized.

@RafaelNH

RafaelNH Mar 2, 2015

Contributor

Yup! This is a matrix trace propriety - trace is invariant to tensor rotations, thus MD=trace(DT)/3=(Dxx+Dyy+Dzz)/3=(e0 + e1 + e2)/3.

# define vector of DKI parameters
MD=sum(dt[0:3])/3
X=np.concatenate((dt,kt*MD*MD,np.array([np.log(S0)])),axis=0)

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

How is b**2 factored into this equation for the KT?

This comment has been minimized.

@RafaelNH

RafaelNH Mar 2, 2015

Contributor

b**2 is factorized on dki design matrix - see function dki_design_matrix

return S
def dki_design_matrix(gtab):

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

It's fine to put this here for now, but we'll eventually need to have just one implementation of this, preferably in the reconst/dki.py file (which doesn't exist yet on this branch...), so we need to remember to remove this from here, and import it from there when we start working on merging these things together.

This comment has been minimized.

@RafaelNH

RafaelNH Mar 2, 2015

Contributor

I agree

def dki_design_matrix(gtab):
r"""
Constructs B design matrix for DKI

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

What's B?

This comment has been minimized.

@RafaelNH

RafaelNH Mar 2, 2015

Contributor

As for DTI, in DKI framework you can define the b-values or the B matrix. This function is computing the B matrix for DKI.

-------
design_matrix : array (N,22)
Design matrix or B matrix for the DKI model
design_matrix[j, :] = (Bxx, Byy, Bzz, Bxy, Bxz, Byz,

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

Maybe write this in latex?

BlogS0)
"""
b=gtab.bvals
bvec=gtab.bvecs

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

PEP8: spaces around = in these two lines.

@arokem

This comment has been minimized.

Member

arokem commented Mar 1, 2015

Did you mean to push the Thumbs.db file into this PR?

In this example we show how someone can simulate the diffusion tensor,
diffusion kurtosis tensor and signals of a single voxel by taking into account
different cellular media using a multi-compartment simulatations. Simulations

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

typo =:> "using a multi-compartment simulation".

In this example we show how someone can simulate the diffusion tensor,
diffusion kurtosis tensor and signals of a single voxel by taking into account
different cellular media using a multi-compartment simulatations. Simulations
are preformed according to:

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

=> "performed"

Here we use the one we created in :ref:`example_gradients_spheres`.
"""
from gradients_spheres import gtab

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

This will run the gradients_spheres example. Is it very fast?

This comment has been minimized.

@RafaelNH

RafaelNH Mar 2, 2015

Contributor

Yes. I use this example also to keep consistency to the simulations shown in simulate.py.

"""
In ``mevals`` we save the eigenvalues of each tensor. To simulate crossing
fibers with two different media (representing intra and extra-cellular media)
4 tissue components are taken in to account (the first two compartments

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

"in to" => "into" (no space).

In ``mevals`` we save the eigenvalues of each tensor. To simulate crossing
fibers with two different media (representing intra and extra-cellular media)
4 tissue components are taken in to account (the first two compartments
correspond to the intra and extra cellular media for the first fiber population

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

=> "...intra- and extra-cellular medium..."

"""
In ``angles`` we save in polar coordinates (:math:`\theta, \phi`) the principal
axis of each compartment tensor. To simulate crossing fibres at 70 degrees

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

Let's go with the US spelling everywhere: "fibers".

This comment has been minimized.

@RafaelNH

RafaelNH Mar 2, 2015

Contributor

Ok

"""
In ``angles`` we save in polar coordinates (:math:`\theta, \phi`) the principal
axis of each compartment tensor. To simulate crossing fibres at 70 degrees
the compartments corresponding to the first fibre are aligned to the x-axis

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

Same here.

In ``angles`` we save in polar coordinates (:math:`\theta, \phi`) the principal
axis of each compartment tensor. To simulate crossing fibres at 70 degrees
the compartments corresponding to the first fibre are aligned to the x-axis
while the compartments corresponding to the second fibre are aligned to the

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

Likewise.

"""
In ``fractions`` we save the percentage of the contribution of each tensor,
which is computed by multipling the percentage of contribution of each fiber

This comment has been minimized.

@arokem

arokem Mar 1, 2015

Member

=> "multiplying"

RafaelNH added some commits May 13, 2015

STYLE: comments on functions compute_Wijkl and single_diffkurt_tensor…
…s were improved, formulas are now latex compatible
BW: remove some unnecessary defaut values on DKI simulates and remove…
… unnecessary inputs of function compute_Wijkl (dipy/sims/voxel.py)
BW: rename the inputs i,j,k, and el of function compute_Wijkl to avoi…
…d future bugs, note that the added for inside the function has index i as well (dipy/sims/voxel.py)
RF: adding an helpers function to check if ground truth directions ar…
…e in the write format. Before this, code was repeated in each type of simulations, thus now the codes are much more easy to read. This helper function will be tested on test modules
Test: testing helpers function to check if ground truth directions of…
… simulations are in the right format (dipy.sims.test.test_voxel.py). Starting dki simulation tests
TESTS: Tests for DKI simulations complete. Three test are added: test…
…_compute_Wijkl, test_DKI_simulations_aligned_fibers, test_DKI_crossing_fibers_simulations
DOC, BW: Minor changes on dipy.sims.voxels.py: 1) Insert new example …
…of the use of function multi_tensor_dki (DOC); 2) insure that dt and kt are arrays in function single_diffkurt_tensors (BW)
RF: Add two inputs in function compute_Wijkl (dipy/sims/voxel.py), no…
…w total DT and MD can given as input, this can be usefull for applications that need to use compute_Wijkl many times (e.g. if simulations need to be repeated several times)
DOC, STYLE: (dipy.sims.voxel.py)Rename some variables and functions t…
…o something more meanfull. Functions compute_Wijkl and single_diffkurt_tensor are know namedkurtosis_element and DKI_signal. Variables Dc and mD are named D_comps standingfor Diffusion tensors for individual COMPartmentS. Added documentation for check_directions.
BW, TEST: Rename dipy.sims.tests.test_voxel functions according to ch…
…anges at dipy.sims.voxel. Adding test for testing kurtosis_elements optional inputs. Remove test_dki which is useless
BW, DOC: Change the gradient table used on DKI simulation example (th…
…is update reduces computing time and clarifies DKI requirements). Added also a figure directive
BF: Fixed an error in computing individual DT on DKI multicompartment…
…al simulations (dipy/sims/voxel.py). A change on all_tensor_evecs outputon the master branch caused an bug on line 385 of voxel.py. Matrices are now multiplied in the right order
NF: Adding the signal simulated from only the dt components and compa…
…ring this to the signals obtain from the dki simulations in DKI simulation example
BW: Output of all_tensor_evecs returns eigenvectors column-wese. dipy…
….sims.voxel and respective tests were updated

@RafaelNH RafaelNH force-pushed the RafaelNH:DKI_PR1_test_simulations_RNH branch from fe9e0ef to 1e11d08 Jun 8, 2015

@RafaelNH

This comment has been minimized.

Contributor

RafaelNH commented Jun 8, 2015

I just add the example comparing DKI simulations with a simulation using only DTI components. I also add a small paragraph explaining the differences. Please let me know if everything is clear.

I fix the function all_tensor_evecs to return the eigenvectors in column-wise and made the necessary changes on the parts of voxel.py that were calling this function.

@RafaelNH

This comment has been minimized.

Contributor

RafaelNH commented Jun 8, 2015

PS: Having the eigenvectors column-wise is also consistent to numpy (np.linalg.eigh).

In this example we show how to simulate the diffusion kurtosis imaging (DKI)
data of a single voxel. DKI captures information about the non-Gaussian
properties of water diffusion which is a consequence of the existence of tissue
barriers and compartments. In this simulations compartmental heterogeneity is

This comment has been minimized.

@arokem

arokem Jun 8, 2015

Member

typo: "this simulations" => "these simulations"

This comment has been minimized.

@arokem

arokem Jun 8, 2015

Member

or "this simulation", if you want the singular

gaussian water diffusion. This feature can be shown from the figure above,
since signals simulated from the DKI models reveals larger DW signal
intensities than the signals obtained only from the diffusion tensor
components.

This comment has been minimized.

@arokem

arokem Jun 8, 2015

Member

Great!

@arokem

This comment has been minimized.

Member

arokem commented Jun 8, 2015

OK - this is ready to go! I'll fix the typo in the example in a subsequent commit, once I merge this one in.

Congratulations @RafaelNH on your first PR! On to DKI fitting!

arokem added a commit that referenced this pull request Jun 8, 2015

Merge pull request #582 from RafaelNH/DKI_PR1_test_simulations_RNH
DKI project: PR#1 Simulations to test DKI

@arokem arokem merged commit d9999dd into nipy:master Jun 8, 2015

1 check passed

continuous-integration/travis-ci/pr The Travis CI build passed
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment