# WIP: Dki #233

Closed
wants to merge 19 commits into
from
+341 −4

## Conversation

Projects
None yet
4 participants
Member

### arokem commented Aug 22, 2013

 This is still WIP, but I think that I am about ready for others to also start taking a look at this. I am relying mostly on these two papers for the implementation details: [1] Lu, H, Jensen, JH, Ramani, A, Helpern, JA (2006). Three-dimensional characterization of non-gaussian water diffusion in humans using diffusion kurtosis imaging. NMR in Biomedicine 19: 236-247 [2] Jensen, JH and Helpern JA (2010). MRI quantification of non-Gaussian water diffusion by kurtosis analysis. NMR in Biomedicine 23: 698-710. However, I am trying to generalize this to cases in which measurements were conducted in different directions in the different b values. Not sure if that's legit.
Member

### arokem commented Aug 22, 2013

 Also - as a comment - this will eventually depend to some degree on stuff that is in this pull request: #230 So that one has to be first merged for this to be finalized!

### MrBago reviewed Aug 22, 2013

dipy/reconst/dki.py Outdated
 X = np.matrix(A.copy()) return np.dot(linalg.pinv(np.dot(X.T, X)), X.T)

#### MrBago Aug 22, 2013

Contributor

I think ols_matrix is a re-implementation of np.linalg.pinv. When the rows of A are linearly independent pinv(A) is (A' x A)^{-1} A' (but I think it's implemented slightly differently to be more numerically stable).

http://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_pseudoinverse#Linearly_independent_rows

#### arokem Aug 22, 2013

Member

You are correct. I wanted a way to perform the operation of linalg.solve for a matrix with dimensionality (x,y,z, b). If we use linalg.solve we would have to loop over voxels, which is what I am trying to avoid here. Does that even make sense though?

Here - the rows will sometime be linearly dependent - in cases where measurements in two b values are done along the same direction. I haven't tested that yet, but that might be a problem (?).

#### MrBago Aug 22, 2013

Contributor

Sorry, the linearly independent thing is a bit of a red herring. You're fine as long as you have at least 6 unique gradient directions per shell, but repeats are fine on top of that. For the 1d case, the following code is equivalent to using np.linalg.solve if X is (b, b) and np.linalg.lstsq if X is (b, m), but it also works when data is (x, y, z, b).

invX = np.linalg.pinv(X)
result1 = np.dot(data, invX.T)

# Or

result2 = np.tensordot(invX, data, axes=(1, -1))

# notice that result2.T will make it (m, z, y, x)

result2 = np.rollaxis(result2, 0, result2.ndim)

On Thu, Aug 22, 2013 at 9:20 AM, Ariel Rokem notifications@github.comwrote:

In dipy/reconst/dki.py:

• is given by:
• ..math ::
•    \hat{\beta} = (A' x A)^{-1} A' y

• A = np.asarray(A)
• X = np.matrix(A.copy())
• return np.dot(linalg.pinv(np.dot(X.T, X)), X.T)

You are correct. I wanted a way to perform the operation of linalg.solve
for a matrix with dimensionality (x,y,z, b). If we use linalg.solve we
would have to loop over voxels, which is what I am trying to avoid here.
Does that even make sense though?

Here - the rows will sometime be linearly dependent - in cases where
measurements in two b values are done along the same direction. I haven't
tested that yet, but that might be a problem (?).

Reply to this email directly or view it on GitHubhttps://github.com//pull/233/files#r5928786
.

#### MrBago Aug 22, 2013

Contributor

Sorry there was a bunch of typos in my email, I've updated the post but I'm
not sure edits get emailed out. To read the cogent version try this link
https://github.com/nipy/dipy/pull/233/files#r5928786.

Bago

On Thu, Aug 22, 2013 at 10:12 AM, Bago mrbago@gmail.com wrote:

Sorry, the linearly independent thing is a bit of a read herring. The
following code is equivalent to using np.linalg.solve if X is (b, b) and
np.linalg.lstsq if X is (b, m) for the 1d case, but it also works when data
is (x, y, z, b).

invX = np.linalg.pinv(X)
result1 = np.dot(data, invX.T)

# Or

result2 = np.tensordot(invX, data, axes=(1, -1))

# if data is (x, y, z, b) result 2 will be (b, x, y, z) to get it to be

(x, y, z, b) we can use rollaxis

# notice that result2.T will make it (b, z, y, x)

result2 = np.rollaxis(result2, 0, result2.ndim)

On Thu, Aug 22, 2013 at 9:20 AM, Ariel Rokem notifications@github.comwrote:

In dipy/reconst/dki.py:

• is given by:
• ..math ::
•    \hat{\beta} = (A' x A)^{-1} A' y

• A = np.asarray(A)
• X = np.matrix(A.copy())
• return np.dot(linalg.pinv(np.dot(X.T, X)), X.T)

You are correct. I wanted a way to perform the operation of linalg.solve
for a matrix with dimensionality (x,y,z, b). If we use linalg.solve we
would have to loop over voxels, which is what I am trying to avoid here.
Does that even make sense though?

Here - the rows will sometime be linearly dependent - in cases where
measurements in two b values are done along the same direction. I haven't
tested that yet, but that might be a problem (?).

Reply to this email directly or view it on GitHubhttps://github.com//pull/233/files#r5928786
.

### MrBago reviewed Aug 22, 2013

 # We will have a separate tensor model for each shell: self.tensors = [] self.sh_idx = [] self.shells = np.unique(gtab.bvals[~gtab.b0s_mask])

#### MrBago Aug 22, 2013

Contributor

This is probably ok, but we should be on the lookout for floating point precision issues with this. IE bvals=[0 1000, 999.9999].

### MrBago reviewed Aug 22, 2013

dipy/reconst/dki.py Outdated
 """ # This is the mean of the diagonal terms and the squared off-diagonals return np.mean(np.concatenate([self.model_params[..., :3],

#### MrBago Aug 22, 2013

Contributor

Should we consider re-ordering the elements of model_params to avoid having to do this concatenate here?

### MrBago reviewed Aug 22, 2013

dipy/reconst/dki.py Outdated
 """ # Extract the ADC for the diffusion-weighted directions : sphere = dps.Sphere(xyz=self.gtab.bvecs[~self.gtab.b0s_mask]) self.tensor_fits = []

#### MrBago Aug 22, 2013

Contributor

Do you want the tensor fit to be an attribute of the model? If anything it should probably be an attribute of the fit. Could it just be a variable here? Same for adc.

Contributor

### MrBago commented Aug 22, 2013

 This looks like a good start! Thanks.
Member

### Garyfallidis commented Sep 6, 2013

 A general comment about this PR is that it does not provide an example. Please add an example in the docs. This will help with the reviewing process. I think you can use the sherbrooke_3shell dataset for this type of model.

### Garyfallidis reviewed Sep 6, 2013

 class DiffusionKurtosisModel(object): """ The diffusion kurtosis model:

#### Garyfallidis Sep 6, 2013

Member

Why here you are not using the multivoxel decorator? I am not saying that you should. But maybe it would reduce some of the code.

Member

### arokem commented Sep 6, 2013

 Yes. Good point. I will do that. On Fri, Sep 6, 2013 at 8:27 AM, Eleftherios Garyfallidis < notifications@github.com> wrote: A general comment about this PR is that it does not provide an example. Please add an example in the docs. This will help with the reviewing process. I think you can use the sherbrooke_3shell dataset for this type of model. — Reply to this email directly or view it on GitHubhttps://github.com//pull/233#issuecomment-23948062 .

### Garyfallidis reviewed Sep 6, 2013

 class DiffusionKurtosisModel(object): """

#### Garyfallidis Sep 6, 2013

Member

I think that this model could inherit from the OdfModel because in DKI you can as well create ODFs. Something to have in mind.

Merged

Member

### Garyfallidis commented Jan 29, 2015

 @arokem Would it be possible to rebase this PR? I think @chantal and the DKI team will try to build on top of it or use some parts. So, it would be nice if it could run as it is now with the latest version. If you think this is a bad idea and they should make their own PR from scratch please say so. Thx.
Member

### arokem commented Jan 29, 2015

 I'll give it a try and report back. On Thu, Jan 29, 2015 at 7:55 AM, Eleftherios Garyfallidis < notifications@github.com> wrote: @arokem https://github.com/arokem Would it be possible to rebase this PR? I think @chantal https://github.com/Chantal and the DKI team will try to build on top of it or use some parts. So, it would be nice if it could run as it is now with the latest version. If you think this is a bad idea and they should make their own PR from scratch please say so. Thx. — Reply to this email directly or view it on GitHub #233 (comment).

Member

### arokem commented Jan 29, 2015

 I've rebased on master. I ended up merging in dti and dki in from master, because this ended up behind master for some reason: I used --theirs to resolve the conflicts on these files, but it seems that this got interpreted as resolving towards the dki branch. Any ideas on how to make the last two commits nicer in the history?

### arokem added some commits Jul 21, 2013

ENH: Added some DKI-type data. Refactored the calculation to rely on …
…eq 38-39.

This is in Jensen and Helpern's paper. The idea being that you compare the
estimate of ADC in different b value to calculate a 'real' ADC and then use
that to bootstrap the calculation of ADK.
RF: moved a lot of the action into the __init__ of the Fit class.
Also : started to implement prediction for DKI.
WIP: Progress on predicting from DKI.
Still need to figure out a way to calculate diffusivity from the model itself.

### arokem added some commits Aug 25, 2013

RF + BF: A few things
1. Rename some stuff
2. Use the MD estimated from tensors, using the same formula as for the ADC
3. In predict: use np.rollaxis to reshape the result of mutliplying the
parameters with the design-matrix, so that the output matches the data shape.

Merged

### hassemlal commented Jul 27, 2015

 Hi, It looks like the pull request #664 supersedes this request, maybe it should be closed?
Member

### arokem commented Jul 27, 2015

 Good call. Been meaning to do this.
Member

### arokem commented Jul 27, 2015

 Closing, because this was superseded by #664