Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

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

Mean Signal DKI #1230

Merged
merged 38 commits into from Mar 12, 2019

Conversation

@RafaelNH
Copy link
Contributor

commented Apr 22, 2017

In this PR you can find the implementation of the mean spherical DKI model. This model provides non-Gaussian indexes that do not depend on the orientation distribution fiber and that are less sensitive to artefacts when compared to standard mean kurtosis.

The missing parts to complete this PR are the following:

  • Add test to cover 100% of the lines of code (97% is the current test coverage)
  • Documentation proofreading
  • Example of usage (this should include the comparison between the novel index and the standard mk, and simulations showing the biological relevance of the mean spherical kurtosis model).

Please let me know if you have any additional comments that you want me to address at this point!

Cheers,
Rafael

@RafaelNH RafaelNH changed the title (WIP) Mean Spherical DKI (WIP) Mean Signal DKI Aug 31, 2017

@skoudoro

This comment has been minimized.

Copy link
Member

commented Apr 19, 2018

Hi @RafaelNH, What is the status of this PR? It would be good if you can complete the missing part :-)
you seem really close to finishing this PR.

@nilgoyette
Copy link
Contributor

left a comment

I didn't test anything, it's simply a code review.

if self.min_signal is None:
min_signal = MIN_POSITIVE_SIGNAL
else:
min_signal = self.min_signal

This comment has been minimized.

Copy link
@nilgoyette

nilgoyette May 1, 2018

Contributor

This condition should be set in the constructor, where you write

self.min_signal = self.kwargs.pop('min_signal', None)

Can't you replace None by MIN_POSITIVE_SIGNAL and do a little trick for the ValueError after?

index = index + (slice(None),) * (N - len(index))
if model_S0 is not None:
model_S0 = model_S0[index[:-1]]
return type(self)(self.model, model_params[index], model_S0=model_S0)

This comment has been minimized.

Copy link
@nilgoyette

nilgoyette May 1, 2018

Contributor

I don't see this class subclassed anywhere. If it's true and there's no future plan to subclass it, you should replace type(self) by MeanDiffusionKurtosisFit to avoid useless complexity.

if not mask[v]:
continue
# Skip if no signal is present
if np.mean(msignal[v]) <= min_signal:

This comment has been minimized.

Copy link
@nilgoyette

nilgoyette May 1, 2018

Contributor

Is there no way to zip over mask and msignal at the same time instead of using the v index?

This comment has been minimized.

Copy link
@RafaelNH

RafaelNH Jan 10, 2019

Author Contributor

Hi @nilgoyette! I am not finding a trivial way to zip over mask and msignal at the same time! Note that these arrays have different dimensions.

This comment has been minimized.

Copy link
@nilgoyette

nilgoyette Jan 10, 2019

Contributor

What do you mean, mask is in 3D and msignal in 4D?

This comment has been minimized.

Copy link
@RafaelNH

RafaelNH Jan 10, 2019

Author Contributor

Indeed - any concrete suggestion on how to zip over mask and msignal simultaneously?

This comment has been minimized.

Copy link
@nilgoyette

nilgoyette Jan 11, 2019

Contributor

I never did that in Python. If you say that you searched and found nothing, I'll believe you. Thank you for checking.

fractions=frac_sph,
snr=None)
# Compute ground truth values
MDgt = f*Di + (1-f)*De

This comment has been minimized.

Copy link
@nilgoyette

nilgoyette May 1, 2018

Contributor

Please stick to the standard that you use everywhere else: spaces between operators. Here and elsewhere.

bmag : int
The order of magnitude that the bvalues have to differ to be
considered an unique b-value. B-values are also rounded to relative

This comment has been minimized.

Copy link
@nilgoyette

nilgoyette May 1, 2018

Contributor

to relative to


def mdki_prediction(mdki_params, gtab, S0=1.0):
"""
Predict the mean signal given the parameters of the mean spherical DKI, an

This comment has been minimized.

Copy link
@nilgoyette

nilgoyette May 1, 2018

Contributor

a GradientTable object

Parameters
----------
params : ndarray ([X, Y, Z, ...], 2)
Array containing in the last axis the mean spherical diffusivity and

This comment has been minimized.

Copy link
@nilgoyette

nilgoyette May 1, 2018

Contributor

in?

return_S0_hat : bool
Boolean to return (True) or not (False) the S0 values for the fit.
args, kwargs : arguments and key-word arguments passed to the

This comment has been minimized.

Copy link
@nilgoyette

nilgoyette May 1, 2018

Contributor

keyword is always written as a single word. If you change it, change it everywhere else in this review.

coefficients of the mean spherical diffusion kurtosis model. Note that
nub is the number of unique b-values
msignal : ndarray ([X, Y, Z, ..., nub])
Mean signal along all gradient direction for each unique b-value

This comment has been minimized.

Copy link
@nilgoyette

nilgoyette May 1, 2018

Contributor

directionS

Voxel with mean signal intensities lower than the min positive signal
are not processed. Default: 0.0001
return_S0_hat : bool
Boolean to return (True) or not (False) the S0 values for the fit.

This comment has been minimized.

Copy link
@nilgoyette

nilgoyette May 1, 2018

Contributor

I find this sentence (and others like it) unclear. Why not a simpler "If True, also return S0 values for the fit"?

----------
.. [1] Henriques, R.N., Correia, M.M., 2017. Interpreting age-related
changes based on the mean signal diffusion kurtosis. 25th Annual
Meeting of the ISMRM; Honolulu. April 22-28

This comment has been minimized.

Copy link
@arokem

arokem Jun 21, 2018

Member

Year?

raise ValueError(mes)

def fit(self, data, mask=None):
""" Fit method of the DTI model class

This comment has been minimized.

Copy link
@arokem

arokem Jun 21, 2018

Member

=> "Fit method of the MDKI model class"

# Skip if no signal is present
if np.mean(msignal[v]) <= min_signal:
continue
# Define weights as diag(sqrt(ng) * yn**2)

This comment has been minimized.

Copy link
@arokem

arokem Jun 21, 2018

Member

I think the "sqrt" in this comment is incorrect?

@RafaelNH RafaelNH force-pushed the RafaelNH:Mean_DKI branch from 9d942d6 to 5d1d3fb Jan 9, 2019

@codecov-io

This comment has been minimized.

Copy link

commented Jan 11, 2019

Codecov Report

❗️ No coverage uploaded for pull request base (master@7b76d74). Click here to learn what that means.
The diff coverage is 98.47%.

Impacted file tree graph

@@            Coverage Diff            @@
##             master    #1230   +/-   ##
=========================================
  Coverage          ?   83.84%           
=========================================
  Files             ?      120           
  Lines             ?    14561           
  Branches          ?     2294           
=========================================
  Hits              ?    12209           
  Misses            ?     1827           
  Partials          ?      525
Impacted Files Coverage Δ
dipy/reconst/dki.py 96.82% <ø> (ø)
dipy/reconst/fwdti.py 94.28% <100%> (ø)
dipy/core/gradients.py 90% <100%> (ø)
dipy/reconst/msdki.py 98.27% <98.27%> (ø)
@RafaelNH

This comment has been minimized.

Copy link
Contributor Author

commented Jan 18, 2019

Hello! Apologies - progress on this was much slower than I've expected. Basically, I've decided to push back on this to double check if recent advances on microstructural modelling were making MSDKI an obsolete and useless technique. However, according to my recent study (Henriques et al., 2019), I can concluded that this type of signal representation approaches are still a valuable tool for the scientific community. Actually, I feel that this is one of the most exciting techniques that I am sharing in an open source environment. Hope you find it interesting as well.

As I've promised you, I've created an example where I've described be relevance of this technique and showed how one can reconstruct diffusion-weighted data using MSDKI. Hope you enjoy.

In my side this PR is ready to go, but let me know if you have further comments.

The tests might not be passing because I am having an intermittent error in dipy.workflows.tests.test_tracking.test_det_track (not related to this PR) in the environment with python 3.5. Should I create a issue on this?

@RafaelNH RafaelNH changed the title (WIP) Mean Signal DKI Mean Signal DKI Jan 18, 2019

@arokem

This comment has been minimized.

Copy link
Member

commented Jan 18, 2019

Rafael -- nice work! I will take a look later today.

Regarding the error that you hit: did you rebase this on top of master?

@RafaelNH

This comment has been minimized.

Copy link
Contributor Author

commented Jan 18, 2019

Hi Ariel! Yes, I did. Is this issue not happening in other PRs?

@arokem

This comment has been minimized.

Copy link
Member

commented Jan 18, 2019

Doesn't seem to be (for example: https://ci.appveyor.com/project/skoudoro/dipy-5dd9b/builds/21721304 from a PR I have in review), but maybe put up an issue about it, and we can try to resolve that together. It might be something intermittent. I agree that it seems unrelated to your changes, so it's a bit puzzling.

@RafaelNH RafaelNH force-pushed the RafaelNH:Mean_DKI branch from b93e42b to 32f8edd Jan 18, 2019

@RafaelNH

This comment has been minimized.

Copy link
Contributor Author

commented Jan 18, 2019

Well - I did it 9 days ago with no conflict! I've just did it again - still no conflict! Let's see if it passes!

@skoudoro
Copy link
Member

left a comment

Hi @RafaelNH, Nice work and tutorial! Thank you for this! After a quick look, I added a couple of comments below. I will go deeper later.

import numpy as np
import matplotlib.pyplot as plt

# Reconstriuction modules

This comment has been minimized.

Copy link
@skoudoro

skoudoro Jan 18, 2019

Member

typo: reconstruction instead of Reconstriuction

three different b-values).
"""

# Sample the spherical cordinates of 60 random diffusion-weighted directions.

This comment has been minimized.

Copy link
@skoudoro

skoudoro Jan 18, 2019

Member

typo: coordinates

@@ -0,0 +1,350 @@
"""

This comment has been minimized.

Copy link
@skoudoro

skoudoro Jan 18, 2019

Member

Can you add # -*- coding: utf-8 -*- at the top of the file? Otherwise, I get the following error on python27 - windows:

$ python doc/examples/reconst_mdki.py
  File "doc/examples/reconst_mdki.py", line 324
SyntaxError: Non-ASCII character '\xe2' in file doc/examples/reconst_mdki.py on line 325, but no encoding declared; see http://python.org/dev/peps/pep-0263/ for details
from dipy.sims.voxel import multi_tensor_dki
from dipy.io.gradients import read_bvals_bvecs
from dipy.core.gradients import (gradient_table, unique_bvals, round_bvals)
from dipy.data import get_data

This comment has been minimized.

Copy link
@skoudoro

skoudoro Jan 18, 2019

Member

get_data is deprecated, can you use get_fnames ?

I do not remember if this deprecation is on master or the last release, I need to check


from dipy.core.gradients import (check_multi_b, unique_bvals, round_bvals)
from dipy.reconst.base import ReconstModel
from dipy.reconst.dti import (MIN_POSITIVE_SIGNAL)

This comment has been minimized.

Copy link
@skoudoro

skoudoro Jan 18, 2019

Member

I think you do not need the parenthesis

args, kwargs : arguments and keyword arguments passed to the
fit_method. See mdti.wls_fit_mdki for details
min_signal : float

This comment has been minimized.

Copy link
@skoudoro

skoudoro Jan 18, 2019

Member

This docstring does not appear as __init__ function parameter.

self.bmag = bmag
self.args = args
self.kwargs = kwargs
self.min_signal = self.kwargs.pop('min_signal', MIN_POSITIVE_SIGNAL)

This comment has been minimized.

Copy link
@skoudoro

skoudoro Jan 18, 2019

Member

ok, I see, I wonder whether it will be better to be explicit and just add this parameter in the __init__ function like you did on line 329

@RafaelNH

This comment has been minimized.

Copy link
Contributor Author

commented Jan 20, 2019

Hi @skoudoro! Many thanks for your comments - I've address them all! Let me know if you have further ones!

@arokem
Copy link
Member

left a comment

Great. I think this is really close to ready to merge. I had just a couple of small comments.

In particular, I wonder whether something could be done to speed up the prediction from this model. Right now, you are looping through voxels (here:

dipy/dipy/reconst/msdki.py

Lines 89 to 106 in 9a7acfc

# Define MSDKI design matrix given gtab.bvals
A = design_matrix(round_bvals(gtab.bvals))
# Initialize pred_sig
pred_sig = np.zeros(msdki_params.shape[:-1] + gtab.bvals.shape)
# looping for all voxels
index = ndindex(msdki_params.shape[:-1])
for v in index:
params = msdki_params[v].copy()
params[1] = params[1] * params[0] ** 2
if isinstance(S0, np.ndarray):
S0v = S0[v]
else:
S0v = S0
pred_sig[v] = S0v * np.exp(np.dot(A[:, :2], params))
return pred_sig
). Would it be possible to do the dot product over multiple voxels at a time?

At the very least, could we use a mask for this function, to avoid predicting for voxels with no meaningful signal?

@@ -1224,7 +1224,7 @@ def dki_prediction(dki_params, gtab, S0=1.):
The gradient table for this prediction
S0 : float or ndarray (optional)
The non diffusion-weighted signal in every voxel, or across all
voxels. Default: 150
voxels. Default: 1

This comment has been minimized.

Copy link
@arokem

arokem Jan 21, 2019

Member

Good catch!

@@ -1286,7 +1286,13 @@ def wrapped_fit_tensor(design_matrix, data, return_S0_hat=False,
return_S0_hat=return_S0_hat,
*args, **kwargs)
data = data.reshape(-1, data.shape[-1])
dtiparams = np.empty((size, 12), dtype=np.float64)
params0 = fit_tensor(design_matrix, data[0:2],

This comment has been minimized.

Copy link
@arokem

arokem Jan 21, 2019

Member

Is this something to do with reusing this function for other models? Not the standard 12 parameter tensor? Why are you taking two rows from the reshaped data, rather than just the first one?

This comment has been minimized.

Copy link
@RafaelNH

RafaelNH Mar 10, 2019

Author Contributor

Indeed, in one of my earliest implementations, I was reusing this function for msdki. I've just realised that I am not using this anymore. Therefore, I've restored dti to the master's branch version.

index = ndindex(msdki_params.shape[:-1])
for v in index:
params = msdki_params[v].copy()
params[1] = params[1] * params[0] ** 2

This comment has been minimized.

Copy link
@arokem

arokem Jan 21, 2019

Member

I think that at least this computation can be taken out of the loop:

msdki_params[:, 1] = msdki_params[:, 1] * msdki_params[: 0] ** 2

for ...
raise ValueError(e_s)

# Check if at least three b-values are given
enough_b = check_multi_b(self.gtab, 3, non_zero=False)

This comment has been minimized.

Copy link
@arokem

arokem Jan 21, 2019

Member

Shouldn't you pass bmag in here?

This comment has been minimized.

Copy link
@RafaelNH

RafaelNH Mar 10, 2019

Author Contributor

Yes! That is correct 👍

(Doctoral thesis). Downing College, University of Cambridge.
https://doi.org/10.17863/CAM.29356
"""
return msdki_prediction(msdki_params, self.gtab, S0)

This comment has been minimized.

Copy link
@arokem

arokem Jan 21, 2019

Member

Things might go faster if msdki_prediction can use the mask.

simulations (for more information on this type of simulations see
:ref:`example_simulate_multi_tensor`). For this example, simulations are
produced based on the sum of four diffusion tensors to represent the intra-
and extra-cellular spaces of two fiber populations. The parameters of theses

This comment has been minimized.

Copy link
@arokem

arokem Jan 21, 2019

Member

"theses" => "these"


"""
For the acquisition parameters of the synthetic data, we use 60 gradient
directions for two non-zero-b-values (1000 and 2000 $s/mm^{2}$) and two

This comment has been minimized.

Copy link
@arokem

arokem Jan 21, 2019

Member

"non-zero-b-values" => "non-zero b-values" (remove the dash between the two parts).

"""
For the acquisition parameters of the synthetic data, we use 60 gradient
directions for two non-zero-b-values (1000 and 2000 $s/mm^{2}$) and two
zero-bvalue (note that, such as the standard DKI, MSDKI requires at least

This comment has been minimized.

Copy link
@arokem

arokem Jan 21, 2019

Member

"zero-bvalue" => "zero b-values"

gtab = gradient_table(bvals, bvecs)


""" Simulations are lopped for different intra- and extra-cellular water

This comment has been minimized.

Copy link
@arokem

arokem Jan 21, 2019

Member

"lopped" => "looped"

MD = dki_fit.md
MK = dki_fit.mk(0, 3)

""" Now we plot the results as a function of the ground truth interstection

This comment has been minimized.

Copy link
@arokem

arokem Jan 21, 2019

Member

"interstection" => "intersection"

@ShreyasFadnavis

This comment has been minimized.

Copy link
Member

commented Feb 7, 2019

Hi @RafaelNH ! I am working on merging your PR! Apart from the minor changes that @arokem has suggested, LGTM! The code, tests and examples look super!

Can you please make the suggested changes by @arokem ?

@arokem

This comment has been minimized.

Copy link
Member

commented Feb 15, 2019

Hey @RafaelNH : any chance you might be able to make these small fixes? This is almost ready to go!

@skoudoro skoudoro referenced this pull request Mar 5, 2019
14 of 15 tasks complete
RafaelNH added 10 commits Apr 20, 2017
RF, TEST: add function that rounds b-values given the order of magnit…
…ude that values should be round

function check_multi_b also uses this function to check the number of unique b-values
added more tests to check both check_multi_b and round_bvals
NF: Add new script that will contain the implementation of the mean s…
…pherical DKI model

File has a priliminary version of the names of the file's functions and object s
TESTS: Add a file that will contain the tests for mdki
this file already have two simulations of spherical tensors for single and multi-voxels
RF, NF, TESTS: Add function that gives the unique b-values in dipy.co…
…re.gradients:

function check_multi_b is now using this function

@RafaelNH RafaelNH force-pushed the RafaelNH:Mean_DKI branch from 1003bde to 792cad4 Mar 10, 2019

@RafaelNH

This comment has been minimized.

Copy link
Contributor Author

commented Mar 10, 2019

@arokem @ShreyasFadnavis ! Minor changes made! Sorry for the delay on this! I've improved the prediction function to avoid voxels looping, removed unnecessary changes in dti.py and added bmag optional input in maximum b-value requirement check!

@skoudoro

This comment has been minimized.

Copy link
Member

commented Mar 11, 2019

Wow, Thanks for this compatibility update on commit ff68f5a!

I will wait until this Travis matrix succeed.

@skoudoro

This comment has been minimized.

Copy link
Member

commented Mar 12, 2019

Thanks @RafaelNH for this nice work!

@skoudoro skoudoro merged commit ae74815 into nipy:master Mar 12, 2019

4 checks passed

codecov/patch No report found to compare against
Details
codecov/project No report found to compare against
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
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
Projects
None yet
6 participants
You can’t perform that action at this time.