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
Mean Signal DKI #1230
Conversation
Hi @RafaelNH, What is the status of this PR? It would be good if you can complete the missing part :-) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I didn't test anything, it's simply a code review.
dipy/reconst/mdki.py
Outdated
if self.min_signal is None: | ||
min_signal = MIN_POSITIVE_SIGNAL | ||
else: | ||
min_signal = self.min_signal |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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?
dipy/reconst/mdki.py
Outdated
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) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
dipy/reconst/mdki.py
Outdated
if not mask[v]: | ||
continue | ||
# Skip if no signal is present | ||
if np.mean(msignal[v]) <= min_signal: |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is there no way to zip over mask
and msignal
at the same time instead of using the v
index?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What do you mean, mask
is in 3D and msignal
in 4D?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Indeed - any concrete suggestion on how to zip over mask and msignal simultaneously?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I never did that in Python. If you say that you searched and found nothing, I'll believe you. Thank you for checking.
dipy/reconst/tests/test_mdki.py
Outdated
fractions=frac_sph, | ||
snr=None) | ||
# Compute ground truth values | ||
MDgt = f*Di + (1-f)*De |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Please stick to the standard that you use everywhere else: spaces between operators. Here and elsewhere.
dipy/core/gradients.py
Outdated
|
||
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 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
to relative to
dipy/reconst/mdki.py
Outdated
|
||
def mdki_prediction(mdki_params, gtab, S0=1.0): | ||
""" | ||
Predict the mean signal given the parameters of the mean spherical DKI, an |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
a GradientTable object
dipy/reconst/mdki.py
Outdated
Parameters | ||
---------- | ||
params : ndarray ([X, Y, Z, ...], 2) | ||
Array containing in the last axis the mean spherical diffusivity and |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
in?
dipy/reconst/mdki.py
Outdated
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 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
keyword is always written as a single word. If you change it, change it everywhere else in this review.
dipy/reconst/mdki.py
Outdated
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 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
directionS
dipy/reconst/mdki.py
Outdated
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. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I find this sentence (and others like it) unclear. Why not a simpler "If True, also return S0 values for the fit"?
dipy/reconst/mdki.py
Outdated
---------- | ||
.. [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 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Year?
dipy/reconst/mdki.py
Outdated
raise ValueError(mes) | ||
|
||
def fit(self, data, mask=None): | ||
""" Fit method of the DTI model class |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
=> "Fit method of the MDKI model class"
dipy/reconst/mdki.py
Outdated
# Skip if no signal is present | ||
if np.mean(msignal[v]) <= min_signal: | ||
continue | ||
# Define weights as diag(sqrt(ng) * yn**2) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think the "sqrt" in this comment is incorrect?
Codecov Report
@@ Coverage Diff @@
## master #1230 +/- ##
=========================================
Coverage ? 83.84%
=========================================
Files ? 120
Lines ? 14561
Branches ? 2294
=========================================
Hits ? 12209
Misses ? 1827
Partials ? 525
|
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? |
Rafael -- nice work! I will take a look later today. Regarding the error that you hit: did you rebase this on top of master? |
Hi Ariel! Yes, I did. Is this issue not happening in other PRs? |
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. |
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! |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
- Can you add your example in dipy/doc/examples/valid_exemples.txt
- Can you add your example on the reconstruction section on dipy/doc/examples_index.rst
doc/examples/reconst_mdki.py
Outdated
import numpy as np | ||
import matplotlib.pyplot as plt | ||
|
||
# Reconstriuction modules |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
typo: reconstruction
instead of Reconstriuction
doc/examples/reconst_mdki.py
Outdated
three different b-values). | ||
""" | ||
|
||
# Sample the spherical cordinates of 60 random diffusion-weighted directions. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
typo: coordinates
doc/examples/reconst_mdki.py
Outdated
@@ -0,0 +1,350 @@ | |||
""" |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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
dipy/reconst/tests/test_mdki.py
Outdated
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 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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
dipy/reconst/mdki.py
Outdated
|
||
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) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think you do not need the parenthesis
dipy/reconst/mdki.py
Outdated
args, kwargs : arguments and keyword arguments passed to the | ||
fit_method. See mdti.wls_fit_mdki for details | ||
|
||
min_signal : float |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This docstring does not appear as __init__
function parameter.
dipy/reconst/mdki.py
Outdated
self.bmag = bmag | ||
self.args = args | ||
self.kwargs = kwargs | ||
self.min_signal = self.kwargs.pop('min_signal', MIN_POSITIVE_SIGNAL) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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
Hi @skoudoro! Many thanks for your comments - I've address them all! Let me know if you have further ones! |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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: https://github.com/nipy/dipy/blob/9a7acfcec8085d0cd9a58156fd1db4e5d7b27a4c/dipy/reconst/msdki.py#L89-L106). 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 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Good catch!
dipy/reconst/dti.py
Outdated
@@ -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], |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
dipy/reconst/msdki.py
Outdated
index = ndindex(msdki_params.shape[:-1]) | ||
for v in index: | ||
params = msdki_params[v].copy() | ||
params[1] = params[1] * params[0] ** 2 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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 ...
dipy/reconst/msdki.py
Outdated
raise ValueError(e_s) | ||
|
||
# Check if at least three b-values are given | ||
enough_b = check_multi_b(self.gtab, 3, non_zero=False) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Shouldn't you pass bmag
in here?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Things might go faster if msdki_prediction
can use the mask.
doc/examples/reconst_msdki.py
Outdated
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 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
"theses" => "these"
doc/examples/reconst_msdki.py
Outdated
|
||
""" | ||
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 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
"non-zero-b-values" => "non-zero b-values" (remove the dash between the two parts).
doc/examples/reconst_msdki.py
Outdated
""" | ||
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 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
"zero-bvalue" => "zero b-values"
doc/examples/reconst_msdki.py
Outdated
gtab = gradient_table(bvals, bvecs) | ||
|
||
|
||
""" Simulations are lopped for different intra- and extra-cellular water |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
"lopped" => "looped"
doc/examples/reconst_msdki.py
Outdated
MD = dki_fit.md | ||
MK = dki_fit.mk(0, 3) | ||
|
||
""" Now we plot the results as a function of the ground truth interstection |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
"interstection" => "intersection"
Hey @RafaelNH : any chance you might be able to make these small fixes? This is almost ready to go! |
…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
…pherical DKI model File has a priliminary version of the names of the file's functions and object s
this file already have two simulations of spherical tensors for single and multi-voxels
function will has to be tested
… fit thanks Ariel Rokem
…re.gradients: function check_multi_b is now using this function
in this way number of shells do not match number of parameters
@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! |
Wow, Thanks for this compatibility update on commit ff68f5a! I will wait until this Travis matrix succeed. |
Thanks @RafaelNH for this nice work! |
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:
Please let me know if you have any additional comments that you want me to address at this point!
Cheers,
Rafael