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

IVIM #1058

Merged
merged 96 commits into from Sep 27, 2016

Conversation

Projects
None yet
6 participants
@sahmed95
Contributor

sahmed95 commented May 25, 2016

Initial simulations of IVIM

Notebook with examples
https://github.com/sahmed95/notebooks/blob/ivim/ivim_dev.ipynb

Shahnawaz Ahmed
@sahmed95

This comment has been minimized.

Contributor

sahmed95 commented May 25, 2016

@RafaelNH @arokem This is a very basic fitting using scipy's nonlinear curve_fit(). I haven't tried two stage fitting (higher and lower bvalues) and tried to keep this PR short. This notebook gives the result for different sets of simulated data and also for one of dipy's dataset. The generated data fits well but dipy's dataset gives poor fitting (perhaps due to very hight bvalues). I will look into Eric's dataset and try again. Till then, comments and suggestions are welcome. :-)

Shahnawaz Ahmed

@sahmed95 sahmed95 force-pushed the sahmed95:ivim_dev branch from 593b5f7 to 951c19f May 25, 2016

mevals = np.array(([D_star, D_star, D_star], [D, D, D]))
# This gives an isotropic signal
signal = multi_tensor(gtab, mevals, snr=None, fractions=[f, 100 - f])

This comment has been minimized.

@etpeterson

etpeterson May 25, 2016

Contributor

I think you want f to be between 0 and 1, at least that's how you defined it in ivim_function. That may be why you're getting overflow errors in the fitting. Not that 0 to 100 is wrong, it just needs to be consistent.

Also, f (perfusion fraction) should be around 10%, not 90%. It should fit it fine but it's just not the typical brain composition to have 90% of the signal flowing and 10% stationary.

This comment has been minimized.

@sahmed95

sahmed95 May 26, 2016

Contributor

Ah, yes. I am changing it to stay consistent. In multi_tensor it is taken as 0 - 100. Should I follow that or use 0 - 1 for "f" and while passing to multi_tensor do a fractions = [100*f, 100*(1-f)] @arokem. What is dipy's take on this ?

@etpeterson Thanks for pointing it out. But even after changing it, the overflow errors persist. Can it be because the data and bvalues not what it is supposed to be for this model. (I am using 'small_101D', I 'll take a look at the data Eric posted)

I think we might need tests to determine the range of f, D_star and D such that results are valid while using this model. Incorrect parameters give a very bad fit.

This comment has been minimized.

@etpeterson

etpeterson May 26, 2016

Contributor

Looks like dipy may be a little conflicted also. Single_tensor uses S0=1 and multi_tensor uses S0=100, though others seem to prefer the 100 area.

This comment has been minimized.

@arokem

arokem May 26, 2016

Member

Fixes to make things more consistent would be most welcome!

On Thu, May 26, 2016 at 9:18 AM, Eric Peterson notifications@github.com
wrote:

In dipy/reconst/tests/test_ivim.py
#1058 (comment):

+def test_nlls_fit():

  • """
  • Test the implementation of NLLS
  • """
  • fimg, fbvals, fbvecs = get_data('small_101D')
  • bvals, bvecs = read_bvals_bvecs(fbvals, fbvecs)
  • gtab = gradient_table(bvals, bvecs)
  • f, D_star, D = 90, 0.0015, 0.0008
  • mevals = np.array(([D_star, D_star, D_star], [D, D, D]))
  • This gives an isotropic signal

  • signal = multi_tensor(gtab, mevals, snr=None, fractions=[f, 100 - f])

Looks like dipy may be a little conflicted also. Single_tensor uses S0=1
and multi_tensor uses S0=100, though others seem to prefer the 100 area.


You are receiving this because you were mentioned.
Reply to this email directly or view it on GitHub
https://github.com/nipy/dipy/pull/1058/files/951c19fe208c5d54250cf7088c16b962fad29d51#r64775504

def ivim_function(bvals, S0, f, D_star, D):
S = np.vectorize(lambda b, S0, f, D_star, D: S0 *

This comment has been minimized.

@etpeterson

etpeterson May 25, 2016

Contributor

Just curious from an implementation standpoint why you are using vectorize and lambda rather than straight up evaluation here.

This comment has been minimized.

@sahmed95

sahmed95 May 26, 2016

Contributor

Ah, I was thinking of vectorizing the code and tried this initially. @arokem pointed out that right now I don't need to worry about that.

@etpeterson

This comment has been minimized.

Contributor

etpeterson commented May 25, 2016

Looks like a good start. Just a few comments beyond the notes above.

The single stage fitting should be fine for noiseless data (I believe) so don't worry about not having the 2-stage yet.

The noiseless tests should fit exactly (or very close) so especially in the 2nd noiseless test the fit isn't right. On a side note you may want to print out or add to the plots the true and fitted values so you can see just how close all the numbers are.

And finally, try using the b-value and b-vector set the Federau, LeBihan, or I used with more b-values below 300 and maxing out at around 1000, it's more suited to these D and D_star values you're using.

flat_data = data.reshape((-1, data.shape[-1]))
# Flatten for the iteration over voxels:
bvals = gtab.bvals
ivim_params = np.empty((flat_data.shape[0], 4))

This comment has been minimized.

@etpeterson

etpeterson May 26, 2016

Contributor

It looks like you're letting curve_fit use ones for the initial guess. The initial guess should be reasonable values, like approximate values from the a paper for example. It also would be nice for the user to send these in if they choose. This could also be why you're getting overflows.

Shahnawaz Ahmed

@sahmed95 sahmed95 force-pushed the sahmed95:ivim_dev branch from c549cef to b8570fb May 26, 2016

@sahmed95

This comment has been minimized.

Contributor

sahmed95 commented May 26, 2016

@etpeterson Updated Ipython notebook with the plots and various values mentioned.
https://github.com/sahmed95/dipy/blob/ipythonNB/dipy/reconst/ivim_dev.ipynb

Plots and values look much better and give very less error. Funny thing happening though.
Estimated f is becoming (1 - f(actual)) and similarly the values of D and D* are getting interchanged after estimation.

Trying to figure that out (perhaps something with curve fit or my function) but otherwise, this looks okay and is working as expected. I guess the next step would be to use the data you gave and test some more or try a two-stage fitting.

@etpeterson

This comment has been minimized.

Contributor

etpeterson commented May 27, 2016

Those noiseless fits look good. That error is on the scale of the accuracy of the nonlinear fit so I'd say it's working. The swapping is probably just because the initial guess is [1,1,1,1] so the algorithm has no way of knowing which should be D and which should be D* so it just picks one. A final check could be if D*>D and if not swap them and recalculate f just to make sure you return physically reasonable values. Also, the 2 stage fitting would probably fix this as well because in a way it finds close to accurate guesses.

@@ -147,7 +146,8 @@ def nlls_fit(data, gtab, jac=False):
bvals = gtab.bvals
ivim_params = np.empty((flat_data.shape[0], 4))
for vox in range(flat_data.shape[0]):
popt, pcov = curve_fit(ivim_function, bvals, flat_data[vox])
popt, pcov = curve_fit(ivim_function, bvals, flat_data[vox],

This comment has been minimized.

@etpeterson

etpeterson May 27, 2016

Contributor

Oh I see you added a guess but your D and D* guesses are still the same here which may be why they're getting swapped.

@etpeterson

This comment has been minimized.

Contributor

etpeterson commented May 27, 2016

Looks like you are setting the initial value, sorry. See my inline comment regarding that.

One other thing, I'd report some kind of norm (L1 or L2 probably) as the error (numpy.linalg.norm), just easier and currently you could get 0 error from positives and negatives cancelling.

ivim_params = np.empty((flat_data.shape[0], 4))
for vox in range(flat_data.shape[0]):
popt, pcov = curve_fit(ivim_function, bvals, flat_data[vox],
p0=[1.0, 0.10, 0.001, 0.001])

This comment has been minimized.

@arokem

arokem May 27, 2016

Member

Can we please pass p0 as an input to this function as a key-word argument (and allow users to initialize this through the IvimModel __init__, again as a kwarg with these values as defaults)

This comment has been minimized.

@sahmed95

sahmed95 May 30, 2016

Contributor

@arokem The dimension of the initial guess vector [S, f, D_star, D] for the parameters should be the same as the number of voxels ?

But, I guess most of the times the user will just give one set of guesses for all the voxels. The code can check for this and pass the p0 array as p0[voxel] if a set of initial guesses are provided. Does it make sense to keep that provision ? @etpeterson

For now, I am having a keyword argument in the fit method and using the same guess for all voxels.

This comment has been minimized.

@etpeterson

etpeterson May 30, 2016

Contributor

I think allowing for a single voxel and a whole array (image) would be good here. You're right that for the one stage fit a single set of guesses is what you'd expect but when you implement the two part fitting you'll probably be sending in an array of initial guesses so it would be good to have that capability.

Data or response variables holding the data. Note that the last
dimension should contain the data. It makes no copies of data.
jac : bool

This comment has been minimized.

@arokem

arokem May 27, 2016

Member

We should probably drop the jac key-word argument. I don't think that curve_fit allows passing a Jacobian matrix anyway

This comment has been minimized.

@arokem

arokem May 27, 2016

Member

If we replace curve_fit with another function (e.g., scipy.optimize.leastsq), it might be worth figuring out what the Jacobian function is. It might speed up the calculations quite a bit.

This comment has been minimized.

@etpeterson

etpeterson May 27, 2016

Contributor

Agreed. I think it would speed it up quite a bit. I'm not sure which optimization routine is best. I think scipy.optimize.leastsq makes it easy to try a few though.

This comment has been minimized.

@sahmed95

sahmed95 May 27, 2016

Contributor

Should I replace the current fitting routine with leastsq or have a different fitting function and let the user chose the fit_method ? (for example "nlls" or "least_sq")

This comment has been minimized.

@etpeterson

etpeterson May 27, 2016

Contributor

I think leastsq is the most flexible one. I believe one of the inputs is the algorithm so maybe the user should just input that (and have a reasonable default)?

This comment has been minimized.

@sahmed95

sahmed95 May 30, 2016

Contributor

@etpeterson
From https://python4mpia.github.io/fitting_data/least-squares-fitting.html

curve_fit is part of scipy.optimize and a wrapper for scipy.optimize.leastsq that overcomes its poor usability. Like leastsq, curve_fit internally uses a Levenburg-Marquardt gradient method (greedy algorithm) to minimise the objective function.

Similarly from https://lmfit.github.io/lmfit-py/model.html

With scipy, such problems are commonly solved with scipy.optimize.curve_fit, which is a wrapper around scipy.optimize.leastsq.

Curve_fit uses leastsq and allows for more flexibility. So having two fitting routines (curve_fit and leastsq) might be redundant.

Also, the type of noise in IVIM data is Gaussian I presume.

This comment has been minimized.

@etpeterson

etpeterson May 30, 2016

Contributor

I'd welcome input from others on this but looking at the documentation here http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html I'd say that minimize looks to be the most mature and flexible routine. I think the big question is which algorithm (lm, dogbox, BFGS, etc) is best (fast and robust) for fitting these data. I think something with bounds and the possibility of a jacobian are the major requirements. I think for now pick one that allows for that and make it work. Once it's working we may want to explore other solvers to see if they're faster or better.

Actually IVIM (and all MRI data) are rician. Just see the add_noise function, it defaults to rician because of that. https://github.com/nipy/dipy/blob/master/dipy/sims/voxel.py#L76. But don't worry about that, assuming the noise is gaussian is fairly safe in most cases and definitely a good place to start. If you're curious try running your tests with gaussian vs rician noise and see how it changes the fit.

"""Ivim model
"""
def __init__(self, gtab, fit_method="NLLS", *args, **kwargs):

This comment has been minimized.

@arokem

arokem May 27, 2016

Member

This seems to "inherit" a lot of its structure from the DTI model. I don't think that we need a the fit_method or the *args and **kwargs inputs here. In particular, you can move the code that is currently in nlls_fit into the fit method of this class.

This comment has been minimized.

@sahmed95

sahmed95 May 30, 2016

Contributor

Yes, I have followed the dti model and I thought that this was a good way to keep the fitting routines separate if we decide to have multiple fitting routines. If we decide to use only one routine [curve_fit or leastsq (which seems to be at the center of all fitting routines)] then I can insert the nlls_fit code here itself. I m now working on "leastsq" as a separate function similar to nlls fit. Once I complete the implementation and write the tests I can move the code in the fit method if we decide to use just one fitting routine.
What do you suggest ?

This comment has been minimized.

@arokem

arokem May 31, 2016

Member

Add a key-word argument for the split_b_val that we use to split the data for the two-stage fitting.

data_in_mask = np.maximum(data_in_mask, min_signal)
# Lookup design matrix for vectorization ?
params_in_mask = self.fit_method(data_in_mask, self.gtab)

This comment has been minimized.

@arokem

arokem May 27, 2016

Member

Insert the nlls_fit code here.

@@ -0,0 +1,71 @@
""" Testing DTI

This comment has been minimized.

@arokem

arokem May 27, 2016

Member

Whoops. Please change the docstring.

gradient_table_from_bvals_bvecs,
reorient_bvecs)
from dipy.sims.voxel import multi_tensor
from dipy.core.sphere import disperse_charges, Sphere, HemiSphere

This comment has been minimized.

@arokem

arokem May 27, 2016

Member

Could you run pyflakes on this? It seems that there are many imports that don't get used in these tests.

bvecs = get_bvecs(N)
gtab = gradient_table(bvals, bvecs.T)
f, D_star, D = 0.06, 0.0072, 0.00097

This comment has been minimized.

@arokem

arokem May 27, 2016

Member

Where are these values taken from? A comment stating this would be useful.

assert_equal(est_signal.shape, data.shape)
assert_array_almost_equal(est_signal, data)
# assert_array_almost_equal(ivim_fit.model_params, [S0, f, D_star, D])

This comment has been minimized.

@arokem

arokem May 27, 2016

Member

This is not working at the moment? Would be good to reinstate this test, and then make sure that it does pass, before this PR is merged.

This comment has been minimized.

@sahmed95

sahmed95 May 27, 2016

Contributor

I had removed it as the values of D and D* were getting swapped. I will correct it and then reinstate this test.
Also, I was of the opinion that this was just for testing and checking out simulations. Do you plan to merge it to the dipy core package ? In that case, there needs to be several other tests and improvements to the code I am sure.

# assert_array_almost_equal(ivim_fit.model_params, [S0, f, D_star, D])
def get_bvecs(N):

This comment has been minimized.

@arokem

arokem May 27, 2016

Member

This function might be useful as a utility function somewhere. @Garyfallidis : do you have any opinion about where something like this might be useful? In sims or somewhere in core?

@sahmed95 sahmed95 referenced this pull request May 30, 2016

Closed

WIP: Ivim implementation #1055

@RafaelNH RafaelNH added the gsoc2016 label May 30, 2016

@RafaelNH

This comment has been minimized.

Contributor

RafaelNH commented May 30, 2016

I just gave a look to the code! Did you try to use a procedure similar to the one suggested by Bihan et al., (1988)? -i.e., fitting first the data with b-value > 200 using the monoexponential fit, then estimate the f-value. This will be a nice initial estimates for your biexponential fit, and might improve the robustness of the estimates with noise. Federau et al., 2012 followed a similar approach.

@coveralls

This comment has been minimized.

coveralls commented Aug 3, 2016

Coverage Status

Coverage increased (+0.1%) to 83.042% when pulling 6913b8e on sahmed95:ivim_dev into 5edee6b on nipy:master.

@sahmed95

This comment has been minimized.

Contributor

sahmed95 commented Aug 3, 2016

@arokem @rafael @etpeterson I have updated and implemented all the changes we discussed. The test coverage is coming out to be 99%.

@coveralls

This comment has been minimized.

coveralls commented Aug 3, 2016

Coverage Status

Coverage increased (+0.1%) to 83.042% when pulling e07a989 on sahmed95:ivim_dev into 5edee6b on nipy:master.

@sahmed95

This comment has been minimized.

Contributor

sahmed95 commented Aug 3, 2016

Updated blog with the example
https://sahmed95.blogspot.in/

@etpeterson

This comment has been minimized.

Contributor

etpeterson commented Aug 9, 2016

I ran the IVIM code and it's looking like it's failing a lot. The pic is from the f parameter and is scaled with black at 0 and white at 0.5. In my experience f should be fairly clean. I'm running it again now with a mask (which I forgot to use last time). The rough code is below. Am I missing something here?

ivimmodel = IVIMModel(gtab)
ivimfit = ivimmdel.fit(img,mask)
ivimparams = ivimfit.model_params

ivim_dipy

@arokem

This comment has been minimized.

Member

arokem commented Aug 10, 2016

Just to be sure that I understand the failure mode: those values of 'f'
should be much lower?

On Wed, Aug 10, 2016 at 1:57 AM, Eric Peterson notifications@github.com
wrote:

I ran the IVIM code and it's looking like it's failing a lot. The pic is
from the f parameter and is scaled with black at 0 and white at 0.5. In my
experience f should be fairly clean. I'm running it again now with a mask
(which I forgot to use last time). The rough code is below. Am I missing
something here?

ivimmodel = IVIMModel(gtab)
ivimfit = ivimmdel.fit(img,mask)
ivimparams = ivimfit.model_params

[image: ivim_dipy]
https://cloud.githubusercontent.com/assets/1086934/17536305/106343d2-5e49-11e6-88e3-cd12a9fc0e8e.png


You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
#1058 (comment), or mute
the thread
https://github.com/notifications/unsubscribe-auth/AAHPNtrQ8EIsewJZQTBy9MLnksfWX7S-ks5qeQXVgaJpZM4Imo2R
.

@RafaelNH

This comment has been minimized.

Contributor

RafaelNH commented Aug 10, 2016

Yup! I was getting the same f profiles with the example data! @sahmed95 can you show us the f maps before the non-linear convergence?

@etpeterson

This comment has been minimized.

Contributor

etpeterson commented Aug 10, 2016

Mostly just much smoother. Sorry I didn't give much info there. Here's a side by side of what my code is generating (right) vs this PR (left). Not that mine is the gold standard by any means but just that the CSF regions have higher f values relative to the tissue and the values seem to vary more naturally. But on a side note, this PR is about 5x faster than mine to run, this whole brain fitting was 87 minutes, so that's good.

I can probably re-run it and save out the linear fitting as well, that's a good diagnosis. Also perhaps the nonlinear fitting status.

@sahmed95 what voxel coordinates were you using for testing? I want to see if I'm getting the same numbers there that you were.

screenshot

@sahmed95

This comment has been minimized.

Contributor

sahmed95 commented Aug 11, 2016

S0, f, D_star, D = 1000.0, 0.132, 0.00885, 0.000921
# params for a single voxel
params = np.array([S0, f, D_star, D])

mevals = np.array(([D_star, D_star, D_star], [D, D, D]))
# This gives an isotropic signal.
signal = multi_tensor(gtab, mevals, snr=None, S0=S0,
                      fractions=[f * 100, 100 * (1 - f)])
# Single voxel data
data_single = signal[0]

data_multi = np.zeros((2, 2, 1, len(gtab.bvals)))
data_multi[0, 0, 0] = data_multi[0, 1, 0] = data_multi[
   1, 0, 0] = data_multi[1, 1, 0] = data_single

This is what I am using for the tests. Can you tell me the slice for which you have posted the images ? (z value). Masking should not have affected the results anyway.

@sahmed95

This comment has been minimized.

Contributor

sahmed95 commented Aug 11, 2016

If Scipy 17 is not used, many times the value for D is becoming -ve and that is driving the f values to be incorrect I guess.

Also, check these two plots for which f value is almost same but the D* value is 10 times more. Is it supposed to mean anything ?

z = 25
plt.scatter(gtab.bvals, data[125, 125, z, :])
plt.scatter(gtab.bvals, data[105, 105, z, :])

plot

plot2

@etpeterson In your code you find 'f' by taking 1 - S0/S0_intercept in
https://github.com/etpeterson/IVIM_fitting/blob/master/IVIM_fit.py#L270
For the second plot, if we do that, then the 'f' value should be very less since the intercept and the value of S0 are almost the same. Can this be where things are going wrong.
Similarly, in your code https://github.com/etpeterson/IVIM_fitting/blob/master/IVIM_fit.py#L268 D* is obtained only by the linear fit I think.

I am not able to think of why the current PR should give incorrect results since it seems to fit data points well.

@sahmed95

This comment has been minimized.

Contributor

sahmed95 commented Aug 11, 2016

Notebook which shows the CSF section and the fitting.

https://gist.github.com/sahmed95/1af7d4550ee3862c25849946427e05d6

@etpeterson

This comment has been minimized.

Contributor

etpeterson commented Aug 11, 2016

I was looking at slice 16 but the results look similar in all slices so looking at this is fine. The CSF is not so concerning because that is typically less interesting but the areas around it I would hope would be fairly robust. Just in general I tend to get suspicious if f is more than about 20% because that's getting fairly high compared to reported values (Federau paper) and it just seems like the brain has less flowing volume than that.

I'm looking at your plots and I'm wondering how they both have similar f values. I'm trying to extrapolate by eye from the b>200 points where that should land at b=0 and it seems like the first plot makes sense but the second should be 4.9% or something like that, not 49%.

Would it be possible to plot the D only curve from the estimate? I've attached an example plot which has a f of 4.7% for comparison. Also, making it semilog is helpful because the D becomes a straight line.

Also images of the estimated f values could be helpful too. I'm trying to generate that also but because it's just a single voxel at a time it took some work to wrap that in a loop. I hope I have it right now though.
amygdala_l_mean

@sahmed95

This comment has been minimized.

Contributor

sahmed95 commented Aug 11, 2016

@etpeterson I am trying out the semilog plots. Meanwhile consider these two voxels :
p1
p2

In one of them D = 0, and the other has D* = 0. But both have high f.

You can use this notebook to check some more plots. With noise, f shoots up to 50% :
https://github.com/sahmed95/notebooks/blob/ivim_profiling/ivim_testing.ipynb

@sahmed95

This comment has been minimized.

Contributor

sahmed95 commented Aug 11, 2016

@arokem
In terms of the code and fitting, I really could not find any issue directly. Two things to note from Federaus's paper :

  1. In Federau's paper, values are clipped - "Values with f greater than 0.3 and D* greater than 0.05 were set (arbitrarily) to zero because they are not physiologic and likely to result either from noise or turbulent cerebrospinal fluid flow."
  2. In the procedure that is followed there, D is just fit once from the linear plot of bvals > 200 and by holding it constant, in the next run, f and D* are estimated. We are running the whole fit again, with the initial guess.

Also, @etpeterson in your code : https://github.com/etpeterson/IVIM_fitting/blob/master/IVIM_fit.py#L261 you fit three times :

i. For S0_prime and D for bvals >= 200 
   D is not fit again.

ii. S0, D*prime for bvals < 200 using guess
   S0 is obtained from here, D*prime from linear fit is not used later apart from plotting.
   f = 1 - S0_prime/S0, f is obtained in this step.

iii. Finally fit D* considering S0, f, D obtained in i, ii

Did I miss something ?
Can this be why the plots dont make sense.

@etpeterson

This comment has been minimized.

Contributor

etpeterson commented Aug 11, 2016

Regarding the plots, if D or D* are 0 then it's just a single exponential so the f is no longer true. If D is 0 then it's 100% D* so the f would be 1, and vice versa. I think this may be the problem though, that it's reducing one of the D values to 0 which means that f becomes extraneous (the intercept then becomes S0*f essentially) and that is why the images look bad but the fits look good.

So to me it's looking like the guesses are good but the final fit is still moving very far from the initial guesses. Federau didn't do a full fit at the end, like you mention, so his preserves D in that way. I do a full fit in my approach but I regularize to the initial parameters so it doesn't depart too far from them. You can think about it as ways to enforce a biexponential fit in data that is very slightly biexponential.

It's nothing that you did wrong or even a bug in the code, if I'm right :) I didn't have this experience even before regularization so perhaps leastsq is more aggressive in reducing parameters to 0 than minimize is.

For me the test to see if this theory is right would be to display the D and f map from estimate_x0 and the fit side by side. In my experience the estimates should be fairly clean looking, and that would be the gold standard D so if the fit and the estimate_x0 results differ then your problem is with the fit changing the initial guesses too much.

@etpeterson

This comment has been minimized.

Contributor

etpeterson commented Aug 11, 2016

Actually I just checked my processing and it finished so here are the f estimates and fits (first image left and right). The estimates look smooth and reasonable, but in about half the cases the fitting diverges from that significantly. And in those regions the D maps of the fit went to 0 (second image left and right). It's looking to me like either fixing D (not fitting it) or regularizing D and maybe f is the way to go.

screenshot

screenshot

@sahmed95

This comment has been minimized.

Contributor

sahmed95 commented Aug 11, 2016

I implemented @etpeterson 's method as a quick fix and the results seem to be better.
The fractions are reasonable (not crossing more than 20%)

Please have a look at this notebook : https://github.com/sahmed95/notebooks/blob/ivim_profiling/ivim_testing.ipynb

And the corresponding code is here : https://github.com/sahmed95/dipy/blob/ivim_linear/dipy/reconst/ivim.py#L270

I am following the exact method followed by @etpeterson in his original code. You may pull that branch and try a full brain fitting.

@sahmed95

This comment has been minimized.

Contributor

sahmed95 commented Aug 11, 2016

So is it just a matter of regions in the data where leastsq jumps to D=0 or D*=0 ?

@etpeterson

This comment has been minimized.

Contributor

etpeterson commented Aug 11, 2016

Wow that was fast! You just added a D* fitting right? I'm kind of surprised that was the fix to be honest. Maybe the initial guesses we had weren't that good. The results you're showing do look good though.

I think either D or D* going to 0 was just a sign that it was failing.

@sahmed95

This comment has been minimized.

Contributor

sahmed95 commented Aug 11, 2016

It's looking to me like either fixing D (not fitting it) or regularizing D and maybe f is the way to go.

Ok. That actually makes sense and I think that's what Federau's paper mentions.

Right now, I can think of three ways to do this :

  1. Using Eric's original code :

    i. We do linear fitting for D (bvals > 200) and store S0_prime
    ii. Another linear fit for S0 (bvals < 200).
    iii. Estimate f using 1 - S0_prime/S0.
    iv. Use least squares only to fit for D_star. 
    
  2. Use the current leastsq fitting for all 4 paramters after (1) taking the S0, f, D_star and D obtained as initial guess. This is like a refined fitting even though we have good estimates for the parameters from 1 and might be unnecesarily time consuming but will surely give the best feasible fit.

  3. Keep the routine as it is. After fitting, deal with anomalous cases (D=0 or D* = 0, f > 50%) by doing a linear fit on a single exponential and clipping f values to 0.

@sahmed95

This comment has been minimized.

Contributor

sahmed95 commented Aug 11, 2016

@etpeterson I didn't just add a D* fit. I followed this :

i. We do linear fitting for D (bvals > 200) and store S0_prime.
ii. Another linear fit for S0 (bvals < 200).
iii. Estimate f using 1 - S0_prime/S0.
iv. Use least squares only to fit for D_star. 
@etpeterson

This comment has been minimized.

Contributor

etpeterson commented Aug 11, 2016

OK. I usually do a final fitting for the whole thing at the end, but it's regularized. Mostly because I don't trust the f or S0 estimate completely which affects D* also. I thought this approach ended up being too through and there's no need to fit bvals < 200 and fit only D* separately was just taking extra time and not really helping. But perhaps that kind of attention is necessary?

@sahmed95

This comment has been minimized.

Contributor

sahmed95 commented Aug 12, 2016

Use the current leastsq fitting for all 4 paramters after (1) taking the S0, f, D_star and D
obtained as initial guess. This is like a refined fitting even though we have good
estimates for the parameters from 1 and might be unnecesarily time consuming but will
surely give the best feasible fit.

This doesn't work and least_squares jumps into the D=0 or D*=0 solution nevertheless. I just tried it and the results are same as before.

This leaves us with either option 1 (only fit D once and estimate f from the intercepts) or option 3 (filter out the unfeasible params)

@sahmed95

This comment has been minimized.

Contributor

sahmed95 commented Aug 12, 2016

Update : Using minimize with bounds too doesnot help and the solutions quickly jump to high values of f or D/D*=0. This happens even if I make the tolerance 1e-10.

@arokem
Is this odd ? Since, even though an almost satisfactory solution is fed to optimize or leastsq, it is aggressively jumping to solutions where D/D* = 0 and f is high.

@etpeterson The last images you have posted, which one of them are the results of this PR ? If the left ones are using this PR then I think the fit is working fine. Only in case of unreasonable signal values f is erroneous. This can be simply tackled by clipping such values of f to 0. or selecting a good mask for the data.

Using the linear fitting method mentioned before also leads to fairly low values of f which is making me think that leaving the fit as it is and filtering the anomalous parameters might be a good option.

OK. I usually do a final fitting for the whole thing at the end, but it's regularized. Mostly
because I don't trust the f or S0 estimate completely which affects D*

This might be an issue as just doing a simple linear fitting for D and S0 and using it to get f does not pass the current test cases. One last option I thought was to set up better bounds like f = 0.3 or D/D* = 0.005.

(mismatch 100.0%)

 x: array([[[[ 990.473441,  972.047586,  954.705322,  938.344793,  922.874338,
           894.281925,  868.362186,  844.656519,  822.793097,  802.470279,
           783.443197,  765.512905,  748.517607,  732.325565,  716.829365,...
 y: array([[[[ 1000.      ,   980.862432,   962.744726,   945.566015,
            929.252227,   898.953732,   871.371977,   846.106374,
            822.820585,   801.232115,   781.10359 ,   762.235448,...

Here you can see both one stage and two stage implemented together. This passes all current test cases but again has the issue of jumping to erroneous paramters in two stage : sahmed95@da39858#diff-399e30b2aa12d6577cd208ab46f99082R245

@etpeterson

This comment has been minimized.

Contributor

etpeterson commented Aug 15, 2016

The images are from the older PR, the f and D values. I haven't tested any newer PRs yet. The left images are from the estimated values and the right are from the fitted. So basically the estimates are good but the fits aren't. This is what you're saying too.

I think the current method doesn't have issues with strong signals but it seems to as there is more noise. So maybe using the normal 2 stage fitting with some SNR based regularization, or even simpler just keeping S0 and D the same as the estimate and only fitting f and D* (which would be faster). The first would be more like my approach and the 2nd more like Federau. I'm not sure which would be better or even if it's worth tweaking at this point or that is more of a post-merge PR.

@arokem arokem merged commit e07a989 into nipy:master Sep 27, 2016

4 checks passed

codecov/patch 95.17% of diff hit (target 80.84%)
Details
codecov/project 80.99% (+0.14%) compared to 5edee6b
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
coverage/coveralls Coverage increased (+0.1%) to 83.042%
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment