# DKI fitting (DKI PR2) #664

Merged
merged 60 commits into from Jul 13, 2015

## Conversation

Projects
None yet
5 participants
Contributor

### RafaelNH commented Jun 11, 2015

 I am starting this work from the code that Prof. Maurizio Marrale kindly supplied. He did a nice work on DKI fitting and estimation of standard DKI rotational invariant measures. Particularly, mean and radial kurtosis are implemented according to the analytical solution proposed by Tabesh et al., 2011 (Magn Reson Med. 65(3):823-36. doi: 10.1002/mrm.22655). Nevertheless, there is still plenty of work to be done on this pull request. I already remove some redundancy of code relative to DTI's module. Now I will add some documentation on the poorly documented functions, rewrite functions in PEP8 style, create some missing parts of the script for better consistency with the module for DTI reconstruction (for example signal predicted from DKI is missing). As suggested by @Garyfallidis on PR #233, I am also creating an example document for DKI usage (doc/example/reconst_dki). Finally, using testing modules, I will check if implementations are correct by comparing the output of the DKI reconstructions with the DKI simulations and comparing the MK/RK estimations to more conservative implementations of DKI invariant measures (Eqs, 11, 12, 13 and 14 of Jensen and Helpern, 2010 - NMR Biomed. 23(7):698-710. doi: 10.1002/nbm.1518)

### RafaelNH reviewed Jun 11, 2015

 # [1] Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011. # Estimation of tensors and tensor-derived measures in diffusional # kurtosis imaging. Magn Reson Med. 65(3), 823-836 #

#### RafaelNH Jun 11, 2015

Contributor

After proper documentation on each function, this will be removed

Member

### arokem commented Jun 12, 2015

 Cool. Looks like you have your work cut out for you here. Don't forget to eliminate the duplication of the DKI design matrix function (which also appears here: https://github.com/nipy/dipy/blob/master/dipy/sims/voxel.py#L536). I think that it makes more sense to put the design matrix definition in this module, and then to import it from here to the simulations.

### omarocegueda reviewed Jun 13, 2015

 import scipy.optimize as opt from dipy.reconst.dti import (fractional_anisotropy, geoderic_anisotropy,

#### omarocegueda Jun 13, 2015

Contributor

Typo: geoderic_anisotropy-->geodesic_anisotropy, this is causing the Travis build failure

### MrBago reviewed Jun 26, 2015

 for kl in range(3): for ll in range(3): multiplyB = B[il][indi]*B[jl][indj]*B[kl][indk]*B[ll][indl] Wre = Wre + multiplyB * W4D[il][jl][kl][ll]

#### MrBago Jun 26, 2015

Contributor

numpy indexing is done by, W4D[il, jl, kl, ll]. The other way works sometimes, but can cause bugs.

### MrBago reviewed Jun 26, 2015

 return dki_params def _wls_iter(design_matrix, inv_design, sig, min_signal, min_diffusivity):

#### MrBago Jun 26, 2015

Contributor

This function seems to be very close to the similar function it dipy.reconst.dti, could this function call the other function to avoid code duplication?

#### RafaelNH Jun 29, 2015

Contributor

The function dipy.recont.dti is implemented specifically for the DTI, so we cannot call this function here. This will be only possible if we generalize the DTI fitting algorithms. For this, we have to remove the decompose_tensor step and change the variable output format. Tensor decomposition have to be done instead in the main. This changes are feasible, however I will recommend doing this in another pull request, this one is getting long.

### MrBago reviewed Jun 26, 2015

 # loop over all voxels for vox in range(len(kt)): R = evecs[vox] dt = lower_triangular(np.dot(np.dot(R, np.diag(evals[vox])), R.T))

#### MrBago Jun 26, 2015

Contributor

I believe there is a function in reconst.dti that computes the lower triangular elements from the evals/evecs that might be more efficient than this.

#### RafaelNH Jun 29, 2015

Contributor

Indeed, I am importing the function lower_triangular from reconst.dti

### MrBago reviewed Jun 26, 2015

 return DKIFit(self, dki_params) def predict(self, dki_params, S0=1):

#### MrBago Jun 26, 2015

Contributor

The predict method should take an optional gtab argument so that predictions can be made on arbitrary gradient sets.

#### RafaelNH Jun 29, 2015

Contributor

Does this really make sense for the predict function defined in the class KurtosisModel? Why should you define a class model with a gtab just to make a prediction with a given model parameters and a different gradient set, if you can do it directly for function dki_prediction without defining any class?

If you think this should be done, this change is trivial. However, I will recommend doing this also in reconst.dti to keep consistency.

Contributor

### MrBago commented Jun 26, 2015

 I like this PR but it's pretty long, I wasn't able to review it in detail. One concern I have is that the math looks very dense. Has anyone looked to see if the long expressions can be vectorized or better implemented using linear algebra operations? When I see V[:, 0] * V[:, 1] * dt[1] + ... I think we should be either using broadcasting or linear algebra but I'm not sure in this case.

### Garyfallidis reviewed Jun 27, 2015

 from .base import ReconstModel def rdpython(x,y,z):

Member

### Garyfallidis reviewed Jun 27, 2015

 def rdpython(x,y,z): r"""

#### Garyfallidis Jun 27, 2015

Member

Hi @RafaelNH and @arokem,

Many of the functions need correction for PEP8 and the math equations will not render as they are in docstrings.
Please you save some good time to correct these.

Also in your functions you usually write WIP in the documentation. This is not a recommended practice. Try to write the documentation of your functions as you go. Otherwise you will end up with bad documentation (quickly written). As you write a function and you still have good memory of what it does it is easier to write the documentation. It makes it easier also for us to give you good reviews.

#### arokem Jun 28, 2015

Member

Hi @Garyfallidis. Thanks for all the comments!

The "WIP" markers here are to tell you that these functions are still work
in progress, and will change substantially before they are really ready for
review. For now, it would be good to focus on the code below line 435 or so
of this file (def mean_kurtosis). In particular, we are still working
through the math of the first few functions, and when we sort it all out,
there will also be clearer documentation here.

On Sat, Jun 27, 2015 at 1:31 PM, Eleftherios Garyfallidis <

In dipy/reconst/dki.py
#664 (comment):

•                          apparent_diffusion_coef, from_lower_triangular,

•                          lower_triangular, decompose_tensor)

+from dipy.sims.voxel import DKI_signal
+from dipy.utils.six.moves import range
+from dipy.data import get_sphere
+from ..core.geometry import vector_norm
+from ..core.sphere import Sphere
+from .vec_val_sum import vec_val_vect
+from ..core.onetime import auto_attr
+from .base import ReconstModel
+
+
+def rdpython(x,y,z):
• r"""

Many of the functions need correction for PEP8 and the math equations will
not render as they are in docstrings.
Please you save some good time to correct these.

Also in your functions you usually write WIP in the documentation. This is
not a recommended practice. Try to write the documentation of your
functions as you go. Otherwise you will end up with bad documentation
(quickly written). As you write a function and you still have good memory
of what it does it is easier to write the documentation. It makes it easier
also for us to give you good reviews.

Reply to this email directly or view it on GitHub
https://github.com/nipy/dipy/pull/664/files#r33415901.

#### arokem Jun 28, 2015

Member

More specifically to your point: these functions are taken from Maurizio's
implementation, so they were copy pasted verbatim in here. As mentioned, we
need to still fully understand these implementations (and verify them
through tests), but once we have that under control, we can also document
them. In principle, I absolutely agree with you that documenting as you
implement is much better.

On Sat, Jun 27, 2015 at 1:31 PM, Eleftherios Garyfallidis <

In dipy/reconst/dki.py
#664 (comment):

•                          apparent_diffusion_coef, from_lower_triangular,

•                          lower_triangular, decompose_tensor)

+from dipy.sims.voxel import DKI_signal
+from dipy.utils.six.moves import range
+from dipy.data import get_sphere
+from ..core.geometry import vector_norm
+from ..core.sphere import Sphere
+from .vec_val_sum import vec_val_vect
+from ..core.onetime import auto_attr
+from .base import ReconstModel
+
+
+def rdpython(x,y,z):
• r"""

Many of the functions need correction for PEP8 and the math equations will
not render as they are in docstrings.
Please you save some good time to correct these.

Also in your functions you usually write WIP in the documentation. This is
not a recommended practice. Try to write the documentation of your
functions as you go. Otherwise you will end up with bad documentation
(quickly written). As you write a function and you still have good memory
of what it does it is easier to write the documentation. It makes it easier
also for us to give you good reviews.

Reply to this email directly or view it on GitHub
https://github.com/nipy/dipy/pull/664/files#r33415901.

#### RafaelNH Jun 29, 2015

Contributor

Hi @Garyfallidis!
As @arokem mentioned the functions marked in WIP are still in the version taken from Mauruzio's. Normally, I write the documentation at the same time that I am writing the functions. In this case, I am writing the documentation at the same time that I am revising the articles and the same time that I am correcting the code :).

### Garyfallidis reviewed Jun 27, 2015

 WIP """ d1mach=np.zeros(5) d1mach[0]=1.1*10**(-308)

#### Garyfallidis Jun 27, 2015

Member

pep8: space between operators. Please correct everywhere.

### Garyfallidis reviewed Jun 27, 2015

 return drd def rfpython(x,y,z):

#### Garyfallidis Jun 27, 2015

Member

This function name is not clear. Is "python" really useful here? And what rf means?

### Garyfallidis reviewed Jun 27, 2015

 def F1m(a,b,c): """ Helper function that computes function $F_1$ which is required to compute

#### Garyfallidis Jun 27, 2015

Member

It is recommended to start the documentation from the first line of the comments.
""" Helper function

### Garyfallidis reviewed Jun 27, 2015

 Parameters ---------- a : ndarray (...)

#### Garyfallidis Jun 27, 2015

Member

Here the parentheses are not informative. Look at other dipy functions for inspirations of how to explain your variables.

#### Garyfallidis Jun 27, 2015

Member

So, what do you mean here that your variables can have any dimensions? Please clarify.

### Garyfallidis reviewed Jun 27, 2015

 -------- Function F_1 is defined as [1]_: \begin{multline}

#### Garyfallidis Jun 27, 2015

Member

I bet that this equation here does not render correctly. Use the math symbol as you do later below.

### Garyfallidis reviewed Jun 27, 2015

 The MK analytical solution is calculated using the following equation [2]_: .. math::

#### Garyfallidis Jun 27, 2015

Member

There is a lot of math here. Make sure that this documentation renders correctly.

### Garyfallidis reviewed Jun 27, 2015

 kurtosis imaging. Magn Reson Med. 65(3), 823-836 """ # Float error used to compare two floats, abs(l1 - l2) < er for l1 = l2 er = 1e-10

#### Garyfallidis Jun 27, 2015

Member

An easier and more robust way to get the machine epsilon across different architectures for a
given float type is to use np.finfo():

print(np.finfo(float).eps)
2.22044604925e-16

print(np.finfo(np.float32).eps)
1.19209e-07

So, 1e-10 that you use here could be not enough for 64bit and too low for 32bit.

Contributor

Cool!!!!

### Garyfallidis reviewed Jun 27, 2015

 return MK def _MK_analytical_solution(dki_params):

#### Garyfallidis Jun 27, 2015

Member

In python we try to not use capital letters in function names. But we use capital letters in classes.

### Garyfallidis reviewed Jun 27, 2015

 return (MD/ADC) ** 2 * AKC def DKI_prediction(dki_params, gtab, S0=150, snr=None):

#### Garyfallidis Jun 27, 2015

Member

Again here you use capital letters in function names. Please correct everywhere.

### Garyfallidis reviewed Jun 27, 2015

 return pred_sig class DKIModel(ReconstModel):

#### Garyfallidis Jun 27, 2015

Member

Non standard naming for class. Please write
class DiffusionKurtosisModel

### Garyfallidis reviewed Jun 27, 2015

 return DKI_prediction(dki_params, self.gtab, S0) class DKIFit(TensorFit):

#### Garyfallidis Jun 27, 2015

Member

DiffusionKurtosisFit

Using longer and more expressive names rather than shortcuts make the code much easier to read.
Having as much readable code as we can is one of our missions. Keep it up!

### Garyfallidis reviewed Jun 27, 2015

#### Garyfallidis Jun 27, 2015

Member

Reduce unneeded blank lines in docstrings.

### Garyfallidis reviewed Jun 27, 2015

 return evals, evecs, kt def dki_design_matrix(gtab):

#### Garyfallidis Jun 27, 2015

Member

Because you are already at dki.py module there is no reason repeating dki in front of function names. So, you can use design_matrix here rather than dki_design_matrix.
Now someone would say that that could be a problem if you want to call design_matrix from two different models e.g. dti and dki. In that case you can write

from dipy.reconst.dti import design_matrix as dti_design_matrix
from dipy.reconst.dki import design_matrix as dki_design_matrix


and resolve any naming issues.

### Garyfallidis reviewed Jun 27, 2015

 assert_) from dipy.sims.voxel import multi_tensor_dki

#### Garyfallidis Jun 27, 2015

Member

remove blank line

### Garyfallidis reviewed Jun 27, 2015

 from dipy.sims.voxel import multi_tensor_dki import dipy.reconst.dti as dti

Member

as above

### Garyfallidis reviewed Jun 27, 2015

 def test_dki_fits(): """DKI fits are tested on noise free simulates"""

Member

simulations

### Garyfallidis reviewed Jun 27, 2015

 non-Gaussian diffusion are of interest since they can be used to charaterize tissue microstructural heterogeneity [Jensen2010]_ and to derive concrete biophysical parameters as the density of axonal fibres and diffusion tortuosity [Fierem2011]_. Moreover, DKI can be used to resolve crossing fibers on

#### Garyfallidis Jun 27, 2015

Member

Correct english please. I would stop at crossing fibers. On tractography doesn't read well.

### Garyfallidis reviewed Jun 27, 2015

 In the following example we show how to reconstruct your diffusion multi-shell datasets using the kurtosis tensor model. First import the necessary modules:

#### Garyfallidis Jun 27, 2015

Member

There is no need to explain the imports. We did that already in some other tutorials. There is no need to repeat the same information.

#### arokem Jun 28, 2015

Member

What if this is the first tutorial a person reads? I don't know that we can
assume that users will read other tutorials before this one, if all they
want is to do DKI.

On Sat, Jun 27, 2015 at 2:41 PM, Eleftherios Garyfallidis <

In doc/examples/reconst_dki.py
#664 (comment):

+fourth-order KT tensors, respectively, and $MD$ is the mean diffusivity.
+As the DT, KT has antipodal symmetry and thus only 15 Wijkl elemments are
+needed to fully characterize the KT:
+
+.. math::

• \begin{matrix} ( & W_{xxxx} & W_{yyyy} & W_{zzzz} & W_{xxxy} & W_{xxxz}
•                & ... \

•                & W_{xyyy} & W_{yyyz} & W_{xzzz} & W_{yzzz} & W_{xxyy}

•                & ... \

•                & W_{xxzz} & W_{yyzz} & W_{xxyz} & W_{xyyz} & W_{xyzz}

•                & & )\end{matrix}

+In the following example we show how to reconstruct your diffusion multi-shell
+datasets using the kurtosis tensor model.
+
+First import the necessary modules:

There is no need to explain the imports. We did that already in some other
tutorials. There is no need to repeat the same information.

Reply to this email directly or view it on GitHub
https://github.com/nipy/dipy/pull/664/files#r33416682.

Contributor

### RafaelNH commented Jun 30, 2015

 @MrBago and @Garyfallidis - thanks for the comments! Indeed this PR was to long! As suggested by @arokem, I removed for now all the functions regarding the estimation of standard kurtosis measures. This PR contains now only the methods for fitting the diffusion and kurtosis tensor from the DKI model (as the PR title suggests). Standard rotational invariant kurtosis measures will be submitted on a third PR after the revision of this one. Please let me know what further changes should be done to complete this PR. Further comments are more than welcome =).

### RafaelNH added some commits Jun 30, 2015

 TESTS, DOC: Removing tests that are not directly relevant for the kur… 
…tosis tensor fit procedures and for the calculation of the apparent kurtosis coefficient values.
 e6558ab 
 DOC: Improve the example of usage of the DKI module. Remove the parts… 
… regarding to the estimation of kurtosis standard measures
 58336cb 
 STYLE: Format dki.py, test_dki.py, and reconst_dki.py according to PEP8 
 5d62abd 
 DOC: Improve documentations according to the comments in Github 
 d157bbd 
 RF: Now DKI module deals with minimun signal in an analogous way than… 
… DTI module
 2e637e7 
 RF, DOC: Remove duplicate function dki_design_matrix on dipy.sims.vox… 
…el. Continue working in documentation of dipy.reconst.dki
 18d9975 
 DOC: Change the DKI fitting example according to the comments done by @… 
…arokem in Github
 20136b7 
 BF: Fixing bug in dki module. Problem happens when importing dki modu… 
…le in dipy.sims.voxel. For now bug was sorted by defining again the dki desing matrix function in dipy.sims.voxel.
 4c25c35 
 DOC: Small changes in the example documentation according to the last… 
… comments done by @arokem.
 3e09d0e 
 TEST, BF: Adding test to evaluate dki_prediction. This was not done b… 
…efore. With this test some bugs where found and fix.
 979c8b1 
 PEP8 
 c3c1fa7 

### RafaelNH referenced this pull request Jul 9, 2015

Merged

#### NF: Function to sample perpendicular directions relative to a given vector #674

 DOC: Minor changes multi shell -> multi-shell 
 2666383 
Contributor

### RafaelNH commented Jul 9, 2015

 Does anyone have suggestions on how to sort the problem that I mentioned on my previous comment? Anyone have any comments before merging this pull request?
Member

### arokem commented Jul 9, 2015

 Oh - sorry - I missed that last comment. We definitely need to resolve this first. Might be a circular import, or some-such. Will take a look.
Member

### arokem commented Jul 9, 2015

 Definitely a circular import: dki => dti => data => sims => dki
Member

### arokem commented Jul 9, 2015

 Still thinking of a solution
 RF: To avoid duplication of code, and circular imports, move dki_desi… 
…gn_matrix

Created a new utils.py module under the reconst namespace to put things like
this.
 1a178b6 

Merged

Member

### arokem commented Jul 9, 2015

 This is my proposal: RafaelNH#1 The idea I have is to create a reconst.utils module that will contain the implementation of dki_design_matrix, and then import from here both in dki and in sims.voxel

### RafaelNH added some commits Jul 10, 2015

 Merge pull request #1 from arokem/dki-move-design 
RF: To avoid duplication of code, and circular imports, move dki_desi…
 29cf6d3 
 BW: Dealing with the min positive signal in the similar way dealed in… 
… module DTI. For this function _min_positive_signal was moved out of TensorModel class. DiffusionKurtosisModel class function fit was also updated with the lastest lines of code on TensorModel class function fit
 f038880 
Contributor

### RafaelNH commented Jul 10, 2015

 I think I answer all your comments regarding this Pull request. Are there any other aspect that I should address at this point? Or is this ready to be merged?
Member

### arokem commented Jul 10, 2015

 This is good to go from my end. Let's give others a few days to take a look. If no one objects, I will merge this on Monday.
Member

### arokem commented Jul 13, 2015

 I'm going ahead and merging this. It will be quickly followed by the next PR in the series, so feel free to comment on any of this for the next PR as well.

### arokem added a commit that referenced this pull request Jul 13, 2015

 Merge pull request #664 from RafaelNH/DKI_PR2_fitting_RNH 
DKI fitting (DKI PR2)
 247450d 

### arokem merged commit 247450d into nipy:master Jul 13, 2015 1 check passed

#### 1 check passed

continuous-integration/travis-ci/pr The Travis CI build passed
Details

Merged

Closed