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

NF: Free water tensor model #835

Merged
merged 111 commits into from Sep 12, 2016

Conversation

Projects
None yet
7 participants
@RafaelNH
Contributor

RafaelNH commented Jan 14, 2016

In this pull request the free water tensor model will be added to Dipy. The model's fit will be based on the method proposed in http://www.ncbi.nlm.nih.gov/pubmed/25271843.

As @arokem suggested this technique will be implemented in a separate reconst module which I named fwdti.py. Since this technique is an expansion of the simple DTI model, its classes will be defined from inheritance of the classes defined on Dipy's DTI module.

At the moment, I just add the structures for the fwdti class. My very next step will be editing class FreeWaterTensorModel.

@etpeterson

This comment has been minimized.

Contributor

etpeterson commented Jan 15, 2016

Looks like you're following the structure from DKI, right? That seems like a good way to start to me. I'll keep checking back to see if I can help out.

@etpeterson

This comment has been minimized.

Contributor

etpeterson commented Jan 20, 2016

Fortunately or unfortunately I'm getting some pressure to get working on this. @RafaelNH if you're otherwise occupied I'll forge on ahead and try to create my own fitting routine to merge into here.

@RafaelNH

This comment has been minimized.

Contributor

RafaelNH commented Jan 20, 2016

Hi @etpeterson! I have to confess that currently I cannot work on this full time, so feel free to carry with the work done here. Actually, github is a perfect platform for developers to work simultaneously in the same scripts. Just keep me updated with your progress to avoid duplicated work. We also have the advantage of being in different time zones, thus if the end of the day you upload your work, I will start my day reviewing your changes and carrying on from what you did.

@RafaelNH

This comment has been minimized.

Contributor

RafaelNH commented Jan 20, 2016

PS: I will be working on this in the following 2 hours.

@arokem

This comment has been minimized.

Member

arokem commented Jan 20, 2016

I'll just add for the benefit of both of you that you can also use our
gitter channel to coordinate: https://gitter.im/nipy/dipy

I've found it useful for coordinating work with others when working
simultaneously on the same thing. And I'd love to stay in the loop on this
work, so it will also allow me to see what's going on :-)

On Wed, Jan 20, 2016 at 2:34 AM, Rafael Neto Henriques <
notifications@github.com> wrote:

PS: I will be working on this in the following 2 hours.


Reply to this email directly or view it on GitHub
#835 (comment).

@RafaelNH

This comment has been minimized.

Contributor

RafaelNH commented Jan 21, 2016

I am getting close to a final version of the WLLS solution of the water elimination model :). It runs without problem, however its outputs are not what I would expect (the new added test is causing Travis to fail).

Anyway, please let me know if you have suggestions on what was done so far or if you detect the problem that is making the test fail. Also coding contributions are welcome. Above, are the parts of the module that still have to be done:

  • Find bug in WLS that is causing fwdti to fail
  • Define the properties of the class FreeWaterTensorFit (that is easy)
  • Add the model's prediction function (analogous to the dti and dki modules)
  • Prediction function have to be tested.
  • Implementing the NLS part of the articles procedure
  • Finding the Jacobian for the NLS fit
@etpeterson

This comment has been minimized.

Contributor

etpeterson commented Jan 22, 2016

Nice work. I copied yours and marked it up with changes and added the prediction function. Until I figure out how to contribute like a normal person I put it in my fork (https://github.com/etpeterson/dipy/blob/fwdti/dipy/reconst/fwdti.py).

I'm not sure I totally followed the steps so some comments may be off base, but I do think the piterations loop is zooming in to the wrong range of f-values.

@RafaelNH

This comment has been minimized.

Contributor

RafaelNH commented Jan 22, 2016

Hi @etpeterson, thanks for your input :). Give a look to Dipy's documentation on how you can contribute with code http://nipy.org/dipy/devel/gitwash/index.html (@arokem, @Garyfallidis - do you have some recommendations of further reading for Eric?). You even can create a pull request of my branch. In this way the revision of your changes will be much more easy. Also it gives me the possibility of commenting code lines using tools of the github website (the comments directly to python should be only the info to stay with the final version of the files for guiding future users). Learning the basics of github is also important because is the only way for formally adding your contributions to Dipy's developing history.

Anyway, I give a look on that you did. I didn't include the function apparent_diffusion_coef - it is a duplicate of function in dti (users can call it from there). The prediction function seems to be correct, however we save to test in using nosetest (see test_fwdti,py for what I tested so far). Regarding to your in script comments, I will answer them here using Github website tools.

# compute F2
S0r = np.exp(-np.array([all_new_params[6],]*nvol))
SIpred = (1-FS)*np.exp(np.dot(W, all_new_params)) + FS*S0r*SFW.T

This comment has been minimized.

@RafaelNH

RafaelNH Jan 22, 2016

Contributor

@etpeterson - this is already correct. First right side equation term already is multiplying by S0, because last element of all_new_params is log(S0), i.e:

np.exp(np.dot(W, all_new_params)) = S0*np.exp(np.dot(W[:, :6], all_new_params[:6]))

Second term FS_S0r_SFW.T needs S0 since the free water contribution computed in line 363 do not take into account S0.

# General free-water signal contribution
fwsig = np.exp(np.dot(design_matrix,
np.array([Diso, 0, Diso, 0, 0, Diso, -np.log(1.)])))

This comment has been minimized.

@RafaelNH

RafaelNH Jan 22, 2016

Contributor

@etpeterson regarding to your question - negative in np.log(1) since design matrix is defined as Bm = [Bxx, Bxy, Byy, Bxz, Byz, Bzz, -1] instead of Bm = [Bxx, Bxy, Byy, Bxz, Byz, Bzz, 1].

This comment has been minimized.

@samuelstjean

samuelstjean Jan 26, 2016

Contributor

log(1) = 0, probable not what you want.

# refining precision
flow = f - df
fhig = f + df
ns = 19

This comment has been minimized.

@RafaelNH

RafaelNH Jan 22, 2016

Contributor

@etpeterson - yes, you are right, this part is different to the article. However, I did this change intentionally to insure that we cover the full f space. For example, assume that the ground truth f is 0.54 however f=0.6 have a local F2 value lower value than 0.5. If you did what the article describes, the second step of the procedure will be searching in the range of 0.55 to 0.65 which doesn't include 0.54.

This comment has been minimized.

@etpeterson

etpeterson Jan 22, 2016

Contributor

Agreed, that is an issue with their method and your implementation is very nice. I'm wondering though if it's a practical issue to never evaluate the end points. According to the article, 18% of the measurements are <0.05 and >0.95 so it seems that it would be good to evaluate those to try to avoid getting stuck in a local minimum.

This comment has been minimized.

@RafaelNH

RafaelNH Jan 22, 2016

Contributor

Lets assume that the ground truth is f=0. In the first step f=0.1 will have the lower F2 value. In the next step the algorithm will sample f from 0.01 to 0.19 with lower F2 value in f=0.01. Following this, the algorithm in theory can converge to 0 with the defined precision. But do not worry, we can add these cases in the testing module. If in practice the algorithm doesn't converge we can add the extreme values, which is computing F2 for the DTI model (f=0) and for a model only having a free isotropic term.

@etpeterson

This comment has been minimized.

Contributor

etpeterson commented Jan 26, 2016

Looks great @RafaelNH! I just got it running on some data I have and at first glance the results look reasonable and it runs just fine. I'll test it some more later and let you know.

raise ValueError(e_s)
def fit(self, data, mask=None):
""" Fit method of the free water elinimation DTI model class

This comment has been minimized.

@samuelstjean

samuelstjean Jan 26, 2016

Contributor

elinimation -> elimination

This comment has been minimized.

@RafaelNH

RafaelNH Jan 29, 2016

Contributor

thanks!

# repeat fw contribution for all the samples
SFW = np.array([fwsig,]*ns)
FS, SI = np.meshgrid(fs, sig)
for r in range(riterations):

This comment has been minimized.

@samuelstjean

samuelstjean Jan 26, 2016

Contributor

you could add break statements when stuff converges, no point in doing iterations that bring little precision

This comment has been minimized.

@RafaelNH

RafaelNH Jan 29, 2016

Contributor

@samuelstjean - Well, for this to converge to a given precision it will require always the same amount of steps, that is why I didn't use break statements. Users can just decrease the number of steps it they want to decrease the precision - see the function's documentation.

This comment has been minimized.

@samuelstjean

samuelstjean Jan 29, 2016

Contributor

I don't really get it, but if you say so, it's only 2-3 iterations anyway after all.

to 3 corresponding to a precision of 0.01.
riterations : inter, optional
Number of iteration repetitions with adapted S0. To insure that S0 is
taken as a model free parameter Each precision iteration is repeated

This comment has been minimized.

@samuelstjean

samuelstjean Jan 29, 2016

Contributor

missing a .

@RafaelNH RafaelNH referenced this pull request Feb 8, 2016

Open

FWDTI nonlinear start #1

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Feb 9, 2016

Hi @RafaelNH, this looks like going towards the correct direction. For merging this soon you need to increase your coverage and report some timings and images with full datasets. After this is done you will need to go ahead with a tutorial. This can be in a second PR.

dipy/reconst/fwdti.py                                183     38     36      8    75%   123-127, 134-136, 153-156, 161, 171-172, 462-485, 544, 581-587, 120->128, 133->134, 148->153, 158->161, 167->171, 459->462, 539->542, 543->544
@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Feb 9, 2016

Also, you may want to ask from @etpeterson to write or extend the tutorial about this and you can mentor him a bit. In that way he can learn how to contribute to Dipy. Sounds like a good plan @etpeterson ?

@etpeterson

This comment has been minimized.

Contributor

etpeterson commented Feb 9, 2016

That works for me @Garyfallidis

@RafaelNH

This comment has been minimized.

Contributor

RafaelNH commented Feb 9, 2016

@Garyfallidis, @etpeterson - Ok! Added some more tests! The only lines that I am not covering is when lsq non-linear procedure does not converge. Anyone have any idea how to test this? Anyway, even if this lines are not working property, it should not be problematic. Since fw_params is already set to the initial guess parameters at the beginning of the function, these lines of code are redundant (perhaps we can even remove the non-covered lines, I was keeping them only to be more consistent to the restore function in dti.py).

I think the next step for me will be running this function in some real data and time its speed. After that I will be more than happy to help @etpeterson with a new PR for the tutorial =). How that sounds?

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Feb 9, 2016

Sounds awesome! Fly high! :)

@RafaelNH

This comment has been minimized.

Contributor

RafaelNH commented Feb 11, 2016

I used the fwdti model to process a slice of the CENIR multi b-values dataset. It took around 90 secs to be processed. Thus, to process the full dataset without a mask it will take around 1.75 hours. Not bad for a first non-optimized version of the procedures and considering that the free water elimination requires a non-linear convergence step.

The free water diffusion tensor derived measures are looking very good. As expected, near to the ventricles and parenchyma, FA values are larger and Diffusion is smaller (Fig1 upper panels) when compared to the values of standard DTI (Fig1 lower panels). There are some voxels with overestimated FA in the ventricles, however these are expected according to Hoy et al., 2014. Users can easily remove these overestimated values of FA by excluding voxels with high free water diffusion volume fractions (Fig.2, idea also mentioned by Hoy et al., 2014.

I added the lines of codes that I used to time the free water elimination model and plot the diffusion measures into doc/examples/reconst_fwdti.py which follows a similar structure of the DKI usage example. @etpeterson - I think you can create the example script in the next PR from this script.

@@ -0,0 +1,608 @@
#!/usr/bin/python
""" Classes and functions for fitting tensors """

This comment has been minimized.

@arokem

arokem Feb 11, 2016

Member

This top-level docstring seems to be a copy-paste residual

which is equal to $3 * 10^{-3} mm^{2}s^{-1} [1]_.
References
----------
.. [1] Pasternak, O., Sochen, N., Gur, Y., Intrator, N., Assaf, Y., 2009.

This comment has been minimized.

@arokem

arokem Feb 11, 2016

Member

I don't think this is the right reference here, considering this is not an implementation of the Pasternak algorithm.

This comment has been minimized.

@arokem

arokem Feb 11, 2016

Member

I see - is this supposed to be a reference for the fact that water diffusivity is set to 3e-3? I am sure there are other, earlier references for that.

@RafaelNH RafaelNH changed the title from (WIP) NF: Free water tensor model to NF: Free water tensor model Feb 11, 2016

assert_almost_equal(FAdti, FAfwe)
assert_almost_equal(Ffwe, gtf)
# Test non-linear fit, when no first quess is given

This comment has been minimized.

@arokem

arokem Feb 11, 2016

Member

Typo: "quess" => "guess"

def test_fwdti_precision():
# Simulation when water contamination is added
gtf = 0.63416 #ground truth volume fraction

This comment has been minimized.

@arokem

arokem Feb 11, 2016

Member

PEP8 (space between # and comment text)

def test_fwdti_predictions():
# single voxel case
# test funtion
gtf = 0.50 #ground truth volume fraction

This comment has been minimized.

@arokem

arokem Feb 11, 2016

Member

Same here (PEP8 comment)

assert_array_almost_equal(fwdtiF.f, GTF)
# 4th error - if a sigma is selected by no value of sigma is given for
# in the non-linear approach to performe restore

This comment has been minimized.

@arokem

arokem Feb 11, 2016

Member

Typo: "performe" => "perform"

In general, maybe: "... if weighting is set to 'sigma', but no value of sigma is given (in case RESTORE is used)."

RafaelNH added some commits Sep 1, 2016

RF: Using the multi_voxel_fit decorator
for this some additional changes were required: 1. S0 estimations is not
required as input, instead it gets the average of b-value=0 signals. 2) S0
is not given as output. 3) NLS always uses WLS as initial guess (the
option of giving a different initial guess had to be removed). 4)
min_signal cannot be defined as min signal of the all dataset, now it is
defined as an optional variable with default value of 1e-6
@coveralls

This comment has been minimized.

coveralls commented Sep 1, 2016

Coverage Status

Coverage decreased (-0.03%) to 82.997% when pulling 90bd462 on RafaelNH:free_water_DTI into 31bc6fe on nipy:master.

1 similar comment
@coveralls

This comment has been minimized.

coveralls commented Sep 1, 2016

Coverage Status

Coverage decreased (-0.03%) to 82.997% when pulling 90bd462 on RafaelNH:free_water_DTI into 31bc6fe on nipy:master.

pred_sig = np.zeros(f.shape + (gtab.bvals.shape[0],))
mask = _positive_evals(evals[..., 0], evals[..., 1], evals[..., 2])
index = ndindex(f.shape)
for v in index:

This comment has been minimized.

@arokem

arokem Sep 2, 2016

Member

Would it be faster to loop over the indices in np.where(mask) instead?

RafaelNH added some commits Sep 3, 2016

@coveralls

This comment has been minimized.

coveralls commented Sep 9, 2016

Coverage Status

Coverage increased (+0.06%) to 83.081% when pulling 8a40caf on RafaelNH:free_water_DTI into 31bc6fe on nipy:master.

1 similar comment
@coveralls

This comment has been minimized.

coveralls commented Sep 9, 2016

Coverage Status

Coverage increased (+0.06%) to 83.081% when pulling 8a40caf on RafaelNH:free_water_DTI into 31bc6fe on nipy:master.

"""
Without spatial constrains the free water elimination model cannot be solved
in data acquired from one non-zero b-value _[Hoy2014]. Therefore, here we
download a dataset that was required from multi b-values.

This comment has been minimized.

@Garyfallidis

Garyfallidis Sep 11, 2016

Member

acquired with multiple b-values

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Sep 11, 2016

This looks very close to be merged. I want to try it with a different dataset too and then will merge (see also comments above). Any other final comments from the rest?

@arokem

This comment has been minimized.

Member

arokem commented Sep 11, 2016

Yup: +1 from me - ready to go as soon as you've had a chance to check things out. @RafaelNH : is there anything we are forgetting, or is this ready to go in from your point-of-view?

@RafaelNH

This comment has been minimized.

Contributor

RafaelNH commented Sep 11, 2016

Hey guys! Yes for my point of view this is ready to be merge! Please let me know if there are any comments that I forgot to address!

Please note that I did some refactoring to use the multi voxel decorator that force me to change some of the features. For example, now when using the non-linear approach this will always use the wls as the first parameters guess. Also S0 is always computed from the mean of the b0 - before I was allowing users to give their own b0 estimate or take it as a model parameter. Because of the multi voxel decorator this is not possible anymore, however I think this is not an issue because the S0 from the mean of the b0 images showed to provide the more robust fits. In general the code is more simplified and easier to read.

In addition, I add the f value correction step in the example (the aspect that I discussed above with @samuelstjean).

@arokem and @Garyfallidis after testing this please feel free to merge this =).

@coveralls

This comment has been minimized.

coveralls commented Sep 12, 2016

Coverage Status

Coverage increased (+0.06%) to 83.085% when pulling 31661d1 on RafaelNH:free_water_DTI into 31bc6fe on nipy:master.

1 similar comment
@coveralls

This comment has been minimized.

coveralls commented Sep 12, 2016

Coverage Status

Coverage increased (+0.06%) to 83.085% when pulling 31661d1 on RafaelNH:free_water_DTI into 31bc6fe on nipy:master.

@coveralls

This comment has been minimized.

coveralls commented Sep 12, 2016

Coverage Status

Coverage increased (+0.06%) to 83.085% when pulling 31661d1 on RafaelNH:free_water_DTI into 31bc6fe on nipy:master.

1 similar comment
@coveralls

This comment has been minimized.

coveralls commented Sep 12, 2016

Coverage Status

Coverage increased (+0.06%) to 83.085% when pulling 31661d1 on RafaelNH:free_water_DTI into 31bc6fe on nipy:master.

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Sep 12, 2016

Good job!

@Garyfallidis Garyfallidis merged commit 630e952 into nipy:master Sep 12, 2016

4 checks passed

codecov/patch 97.03% of diff hit (target 80.96%)
Details
codecov/project 80.99% (+0.02%) compared to 31bc6fe
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
coverage/coveralls Coverage increased (+0.06%) to 83.085%
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment