it's finally here

[etpeterson]

@arokem I finally got this set up and uploaded! I think it works right, though it's still messy. IVIM_fit is the mail file and fitting_diffusion is the supporting file where the important stuff happen. Happy to explain or help, just let me know.

[arokem]

Nice. Thanks for getting that up. What do you think about having a Google Summer of Code student work on integrating this into Dipy? We're in the process of making a list of potential projects over here: https://github.com/nipy/dipy/wiki/Google-Summer-of-Code-2016

[etpeterson]

Sure, it would be great to get it integrated with dipy and it probably would be a good summer project. I could see a few stages of development and even beyond what I've done so far like adding the jacobian to the fitting to speed it up, etc. I can spend some time supervising too but it would be good if you or someone else more familiar with dipy and its development was also involved. How about I write up a description and send it along, kind of like the other examples and you let me know what you think?

[arokem]

Fantastic. If we could put you in as a co-mentor, together with me, that would be great.

On Wed, Feb 3, 2016 at 9:41 AM, Eric Peterson notifications@github.com wrote:

Sure, it would be great to get it integrated with dipy and it probably would be a good summer project. I could see a few stages of development and even beyond what I've done so far like adding the jacobian to the fitting to speed it up, etc. I can spend some time supervising too but it would be good if you or someone else more familiar with dipy and its development was also involved. How about I write up a description and send it along, kind of like the other examples and you let me know what you think?

— Reply to this email directly or view it on GitHub https://github.com/etpeterson/IVIM_fitting/issues/1#issuecomment-179367452 .

[quantshah]

Hi @etpeterson. I have been getting familiar with dipy coding style and diffusion imaging as I started contributing to open source recently. GSoc is a great opportunity to get involved with the open source community and I am interested in contributing to dipy. I worked on #932 and started looking at implementations of Rohde and CHARMED as @arokem suggested. I would like to work on this meanwhile as I gain better understanding of the papers suggested by arokem. I used your code to make IVIM_fit a function and wrote a basic ipython example which outputs an image file - https://github.com/sahmed95/IVIM_fitting/tree/dipy-ivim. Let me know if you can mentor me to integrate this with dipy and if I should open a PR

[etpeterson]

Hi @sahmed95. It would be a huge help to get this going, and I'm happy to mentor and work on this as well. I think a PR would be a great place to start but I'm also fairly new to dipy so @arokem is the real expert with the integration so I'd be happy to hear what he thinks too. IMHO the place to start is a PR and model it after the free water DTI PR even though that is not complete yet (#825).

Regarding the code itself, I think a lot of the options in the code should be removed both because it's excessively complicated and it doesn't fit into dipy as well. I think a two-stage fitting where the high b-values are fit linearly first and then the full model is fit is a good way to start. This two-stage fitting is similar to what #825 has so that's why I mention it.

[quantshah]

Thanks. I have made a separate branch and worked out a simple fit() to run your script as a function. I wanted @arokem and you to know before I started changing the code. @arokem suggested that "it might be worth putting in a few plots. For example, what are the parameters in different parts of the image, how well does the model fit the data, etc."

I will try out a few plots, refactor the code a little more, and open a PR here. Btw the link to #825 is taking me to my own PR. Did u mean #835? RafaelNH's PR?

[etpeterson]

Yes sorry, you're right #835 is the one.

I definitely agree with @arokem on the plotting. It's kind of a delicate fit sometimes so plots are invaluable, especially when the results aren't what you expected.

[quantshah]

Hi, couldn't work much through the weekend but here are some plots - Ipython Notebook using the plotting code @etpeterson wrote in ivim_fit. Each is a plot of the items in np.ndindex(shape3d). Will work more on cleaning up the code this week and then open a PR. @arokem.

[arokem]

The notebook is great -- exactly the kind of reality check we need at this point. I'll defer to @etpeterson to say whether these plots look as they should -- I don't have much experience with this particular model. One thing that occurs to me is that you might want to create some functions that simulate the IVIM model results without noise. This is invaluable for testing the model fitting procedure, and is also a good way to teach/learn about the model in small self-contained scripts. Here, it will also be useful to make sure that you get the same results as you start refactoring the model code to match the dipy API (e.g., inheriting from ReconstModel and ReconstFit).

On Sun, Mar 6, 2016 at 3:42 PM, Shahnawaz Ahmed notifications@github.com wrote:

Hi, couldn't work much through the weekend but here are some plots - Ipython Notebook https://github.com/sahmed95/IVIM_fitting/blob/dipy-ivim/ivim_example.ipynb using the plotting code @etpeterson https://github.com/etpeterson wrote in ivim_fit. Each is a plot of the items in np.ndindex(shape3d). Will work more on cleaning up the code this week and then open a PR. @arokem https://github.com/arokem.

— Reply to this email directly or view it on GitHub https://github.com/etpeterson/IVIM_fitting/issues/1#issuecomment-193016965 .

[etpeterson]

Hey @sahmed95, that is great. The plots are really the best way to see if it's fitting well or not. It looks like most of the fits worked but the b-values are for the most part higher than what we really want for IVIM. Here are the values I used for my experiment: 0, 10, 20, 30, 40, 60, 80, 100, 120, 140, 160, 180, 200, 300, 400, 500, 600, 700, 800, 900, 1000. This is because IVIM is focused on slow flow, and diffusion coefficient of free water (more static or slowly diffusing that slow flow) is ~3e-3 so if you plot that out it drops to near zero by b>1000 which means that flowing blood will decay even faster, hence the focus on lower diffusion values.

@arokem had a good suggestion to simulate data and then try to fit it. If you use the b-values above, and reasonable other parameters would be: 0.01<f<0.2, 1e-2<D*<3e-3, and 1e-4<D<1e-5.

[quantshah]

Thanks @arokem @etpeterson for your inputs. Should I try plotting more with the parameters suggested by Eric?

"create some functions that simulate the IVIM model results without noise"

How?

What can be the next step here? Should it be similar to dti or dki fitting? ReconstModel works with gradient table (gtab). Should the IVIM also inherit this and work with gtab?

[etpeterson]

Basically you want to create diffusion signal curves using the equation S(b)=S0(fe^-bD*+(1-f)e^-bD). If you use the values I suggested above then you'll get a signal value S for each b-value. And this you can fit as if it's any other signal coming from an image series. I think you should basically be able to use the fitting functions you created to fit these curves. And hopefully the fitted values should be similar to what was generated. Does that make sense?

[quantshah]

Thanks. It makes sense now after I could get a quick read of the paper you sent. I will generate the data and try plotting.

[quantshah]

@etpeterson @arokem I generated signals corresponding to the b-values and taking random parameters and got a plot of signal vs bvalues. Link

I am not sure about what S0 values to take so I took random values between 0-1. Next step is fitting.

To fit, I have to get an img from the generated signals similar to an .nii file img as your code takes that as input. Or I can take the part where you fit in begin processing and modify it such that img takes the signal values generated. Not sure how to go about it. Just wanted to know if I am on the right track.

[etpeterson]

You're on the right track.

The S0 value can be anything really, it's just a scaling but I always normalized to the maximum value so S0 was always close to 1, probably 0.95<S0<1.05.

And the other change is that for each parameter set you need to generate the intensities at each b-value. In other words in your generate_signal function you'd input single numbers for the parameters like f, and the vector of b-values so the function would return a vector of values the same size as the b-values.

[quantshah]

@etpeterson Thanks. I have changed the generate_signal function to return a series of values for a b-value array. I am now making a minimalistic fit function which only has the part where processing and fitting is performed. It will take bval, img and mask as input and return the fitted parameters (S0, f , D, D_star). Example :

model = ivim.fit(img, bvals, mask)
S0, f, D, D_star = model.get_fit_params()
model.save('fit.nii')
model.plot()

Does it make sense to do that ? I am thinking from the user's perspective as to what will be the easiest and most flexible.

@arokem For integration to dipy, I think we need to discuss what parameters the function will take and what it gives as output. Is there a similar module in dipy which can be followed?

[etpeterson]

Looks good. I'd say for testing purposes just taking in an array of values and an array of b-values and outputting and array of fits is the cleanest way. Especially if you want to plot the fits and things like that. That would be the internal function that does the fitting and the next function outside of that would be the one that reads and writes images.

[arokem]

@sahmed95 : best to look to @rafaelnh's recent PR here: https://github.com/nipy/dipy/pull/835

That gives you a sense of the API.

In general, we need to follow the structure that is outlines in reconst.base:

https://github.com/nipy/dipy/blob/master/dipy/reconst/base.py

We inherit a ReconstModel (In this case maybe IVIMModel?) and a ReconstFit (IVIMFit?). The Model object carries information about the measurement and the settings of the model, while the Fit object carries the parameters and derived measures (take also a look at the dipy.reconst.dti module to get a sense of how that works).

[quantshah]

@etpeterson had already suggested me this and I am going through the code. This is great as I can follow the same structure for the IVIM model. I will have time over the weekend to read and come up with a basic plan for IVIM fitting based on the tensor fit module for a Gsoc application. Thanks a lot.

[quantshah]

Hi, @arokem @etpeterson. I found some time and worked on writing classes following the diffusion tensor module. model class - IvimModel. I have to change the IvimModel.fit() to take data instead of putting it into IvimModel as pointed out by Garyfallidis. I also have to write tests and work on implementing checks etc

This is just a basic attempt but I feel that I am getting a better understanding of the code as well as dipy's coding style. I have some doubts regarding writing tests and how to use the generated signals for finding errors in fit which I will try to solve by reading the ivim paper once again. In particular given b values and chosen parameters, I will generate signals using ivim_prediction. After that how do I use this array of signals? Do I convert it into a .nii file like this and pass it to an IvimModel again and fit?How to analyse the errors? RMSE ?

The ipython example and usage can be found here - Notebook I made a new branch called dipy similar to reconst : tree

[arokem]

Yeah - this is really good progress! For testing purposes, and for readability, it would help to break down some of the longer functions into sub-functions, so that we can better understand the different parts. I guess your questions about nii files have been answered on the dipy gitter channel?

It's starting to get close to the point where you might want to start working on this as a branch off of dipy, rather than this repo. Do you agree @etpeterson ?

In terms of analysis of the errors: we want to make sure that when we simulate data with this model, with a very high SNR (say with zero noise) and fit the simulated data with the model, that we get back the parameters that we put into the simulation. Or at least very close. We would also want to make sure that the model predictions from this "true model" are very similar. One way to check that is to calculate the RMSE of the model prediction, relative to the simulated data.

[etpeterson]

Agreed, it is looking really good, and I think you're right @arokem, it's time to integrate it closer with dipy.

And @arokem is exactly right about the tests, if it doesn't fit noise free data very well then there's likely something wrong with the code. RMSE is a good metric, but because of the wide range of parameter values between S0, f, D*, and D it probably would be good to normalize to the actual value (the one you chose to generate the signal with) so they're all scaled similarly.

[quantshah]

ivim branch in dipy: https://github.com/sahmed95/dipy/tree/ivim/dipy/ivim. I will work on that branch from now on with a dev folder inside the dipy module for development related things (ipython notebook, quick fits, plots etc)

@etpeterson I have prepared the first draft of a Gsoc proposal for integrating ivim to dipy. Please feel free to edit it and add anything that you feel is important. This is just a first draft: Proposal

[RafaelNH]

GSoC_proposal.pdf I added here the proposal that I wrote last year - I thought that this could be useful as a proposal example.

I think you should explain better the relevance of the work (see my comment to the proposal).

Regarding to your branch, I will suggest you to submit a pull request to dipy instead. You can name your pull request as "WIP - (IVIM fitting)" to inform that the pull request is still work in progress. It will make your work much easier to review.

[quantshah]

@RafaelNH Thank you so much for taking a look. Your comments are really helpful and I am sure with the guidance from all three of you, I will be able to write a solid proposal for GSoc and work on this successfully over the summer. I am very positive about the IVIM project and also feel that this will equip me to work on the eddy correction project later which I had started some time back. I will act on your suggestions in the draft asap. Your last year's proposal is a big help. I will open a PR to dipy today and start working on that branch.

[quantshah]

From Le Bihan's original paper on IVIM - one of the prime technical limitation seems to be eddy currents. Coincidently, I had started with Eddy correction using Rohde as the first project in Dipy. If this limitation still applies to IVIM, it would be interesting to see if an eddy correction followed by fitting an ivim model results in anything significant. I haven't had much time to work on that, but as @arokem told me Bob Dougherty has allowed the use of his matlab code for a python implementation which should make it easier to implement in dipy. I don't know if its possible but if I find that the IVIM work leaves me with time to simultaneously look at Rohde, do you suggest that I mention it in my GSoC application. (Not committing to it obviously)

[etpeterson]

In general eddy currents effect higher b-values more than lower ones, so roughly b>500 I'd say could need correction. So it probably doesn't corrupt IVIM too much these days. That said I always do eddy current correction and motion correction just in case.

[quantshah]

Ok, got it. In any case, I think dipy does not have an eddy correction module. I mentioned it because I was reading the papers @arokem sent me on eddy correction and Le Bihan mentioned it under the section Technical limitation. But that was a long time back I guess. :smile: We have better tech now.

[arokem]

That's all correct. I don't think that mentioning eddy currents as an issue for this method (as well as others) and the potential correction as a "stretch goal", plus pointing to work you've already done on that would hurt your application.

On Tue, Mar 15, 2016 at 4:38 PM, Shahnawaz Ahmed notifications@github.com wrote:

Ok, got it. In any case, I think dipy does not have an eddy correction module. I mentioned it because I was reading the papers @arokem https://github.com/arokem sent me on eddy correction and Le Bihan mentioned it under the section Technical limitation. But that was a long time back I guess. [image: :smile:] We have better tech now.

— You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub https://github.com/etpeterson/IVIM_fitting/issues/1#issuecomment-197069545

[quantshah]

@arokem @etpeterson @RafaelNH I have updated the proposal draft and have filled in most of the details including a rough timeline. If you would take a quick look - Proposal it would be helpful. I have to add references and change some things here and there.

Also, there were some issues with rebasing so I deleted my old dipy repo and created a new one where I will be updating the ivim code under a new branch.

[arokem]

Next time you have issues with rebasing, give me a buzz. In the long run, we'll want you to be able to do this on your own. We can have a hangout together to sort you out.

[quantshah]

Rebasing seemed tricky. I followed this tutorial the first time I had to rebase. This was due to an old pep8 fix which I was working on but someone had already merged it. Issue-930. Thanks. If I run into trouble again, I will contact you over hangouts.

[quantshah]

@arokem I want to generate bvalues and bvectors for testing the ivim fit model. Right now, following the dti code my IvimModel takes a gtab instance. I have the bvalues as suggested by @etpeterson and I would just like to know if I am right about the bvectors. Please correct me if I am wrong.

The bvecs are a set of unit vectors with the shape (N,3) where N is the number of bvalues. Now for artificially generated data, how do I generate these bvecs? Do I randomly make a set of unit vectors?

[arokem]

On Sun, Mar 20, 2016 at 9:45 AM, Shahnawaz Ahmed notifications@github.com wrote:

@arokem https://github.com/arokem I want to generate bvalues and bvectors for testing the ivim fit model. Right now, following the dti code my IvimModel takes a gtab instance. I have the bvalues as suggested by @etpeterson https://github.com/etpeterson and I would just like to know if I am right about the bvectors. Please correct me if I am wrong.

The bvecs are a set of unit vectors with the shape (N,3) where N is the number of bvalues.

Yes - that's correct.

Now for artificially generated data, how do I generate these bvecs? Do I randomly make a set of unit vectors?

You want the gradient directions to be as distinct from each other, while also taking into account the antipodal symmetry of the diffusion signal. One way to do that is to simulate the dispersion of charged particles on a sphere. Fortunately, you won't have to write the code to do that, because we already have that as part of Dipy. Check out this tutorial for details: http://nipy.org/dipy/examples_built/gradients_spheres.html#example-gradients-spheres

— You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub https://github.com/etpeterson/IVIM_fitting/issues/1#issuecomment-198961240

[quantshah]

Thanks. I could make a simple plot for generated data which you can see here

Does it look like a good fit? I found some time over the weekend and went through the dti implementation. As suggested, I made some changes such that Ivim follows dti. I still feel that the code needs work before I open a PR in dipy but with this simple implementation we can create an IvimModel, do a fit with some data and finally plot the results. @etpeterson @arokem

iv = IvimModel(gtab, fit_method="two_stage_fit")
fittedmodel = iv.fit(data, gtab) #fit function needs further change
fittedmodel.plot(1) #Number of plots

[etpeterson]

Hey @sahmed95, the final fit looks fairly good but f and D* seem fairly far off from what you put in. Not sure what that could be but the regularization may just be too high. Also as far as I can tell you aren't using the tensor data and it looks like the fitting function couldn't fit tensors correctly anyway. Is that correct?

On a related note (I'm not sure you want to go here yet), to simulate tensor data for IVIM you probably want to start with two sets of eigenvalues and eigenvectors - one for D* and the other for D - and convert those back to directional D* and D values.

And another thought I had when glancing at the code is that for the initial fit (fitfun_S0_D) you did what I did which is just use the same fitting routine. That's fine and should work but it may be easier and would certainly be faster if you could fit using dipy because that's just a normal linear diffusion fit. The only caveat is that it needs to return a slope and intercept (not assume an intercept of 1) so I'm not sure if dipy can do that or not - I just haven't tried/looked.

[quantshah]

Thanks for taking a look. You are right about using your code, I haven't changed much for now and the core still remains your initial script. I have merely tried to make it closer to dipy. I could only work during the weekends on this due to my upcoming midterm exams. I will start from scratch anyway once my semester ends and will be rewriting the fitting routine.

I am currently not using tensor data and I will look into it.

@etpeterson Even I was thinking of the point you made about regularization and catching poorly fitting regions. It certainly seems important (Atleast that's what this example shows) and I will work on it. I am also exploring dipy's interpolate. For faster fitting routines, I will be trying out the ones in Scipy. Orthogonal Distance regression seems interesting. @arokem It uses a Levenberg-Marquardt-type algorithm as you were discussing on the gitter channel. Once a basic fitting routine is complete, we can look at other routines.

@arokem @etpeterson I have made all the changes in the proposal and if you find some time, it would be great if you could take a look. Proposal Since I will be caught up with exams, I have tried to nearly complete the draft and put it up on the GSoc page. Please let me know if something needs to be added or changed. Thanks.

[etpeterson]

Of course, I totally understand you have school and other things happening. I was just looking at your notebook and ended up with thoughts and comments so I thought I'd write them down.

I'll check out the proposal today or tomorrow.

[quantshah]

:+1: :smile: Thanks a lot. Also, you had mentioned about getting the data from Federau's paper. I have put it as a stretch goal and if we can get the data then it will be a good way to demonstrate the use of the package.

[arokem]

Good to leave something for Summer :-)

On Mon, Mar 21, 2016 at 4:19 PM, Eric Peterson notifications@github.com wrote:

Of course, I totally understand you have school and other things happening. I was just looking at your notebook and ended up with thoughts and comments so I thought I'd write them down.

I'll check out the proposal today or tomorrow.

— You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub https://github.com/etpeterson/IVIM_fitting/issues/1#issuecomment-199531655

[quantshah]

@RafaelNH Thanks for all your help and suggestions. Looking forward to working on this. You may find the final proposal here

[quantshah]

@etpeterson (This might not be completely related to the code, just curious about the model)

1. With bvalues > 200, the plots are nearly a straight line. Approximating the exponential terms in S to 1st order we get a simple linear function S0 (1 - b ( D + f( D* - D ) ). But while plotting for the D, D* values as suggested by you, the signals obtained by this linearization was falling faster that the actual signal. I was thinking of fitting bval <200 and > 200 seperately and thought I could aproximate the equation. Any thoughts?
2. D* = Lv/6 + D where v is blood velocity, L is mean capilary segment. If we have some average values for L and v for a part of the body (say brain/liver) can we feed it in and try to start off with good initial guesses for fitting.
3. Will adding in a function to return the velocity v be of use? Can this model can be used for visualizing blood flow in cells or explore capillary structure

[quantshah]

@arokem Do you suggest a new [WIP]-PR in dipy for this ? I am working on my own repo but there will be lot of changes, dirty code etc. and I want to have a clean documented patch before I open a PR so that it's easier for other's to review. Is it necessary to open a PR in dipy for my proposal to get accepted or a link to my own repository is fine? [https://github.com/sahmed95/dipy/tree/ivim/dipy/ivim]

[arokem]

Link to your own work here is fine for now, especially given you already have links to PRs you made on Dipy.

Open the PR on Dipy when you feel comfortable doing that.

On Thu, Mar 24, 2016 at 1:48 PM, Shahnawaz Ahmed notifications@github.com wrote:

@arokem https://github.com/arokem Do you suggest a new [WIP]-PR in dipy for this ? I am working on my own repo but there will be lot of changes, dirty code etc. and I want to have a clean documented patch before I open a PR so that it's easier for other's to review. Is it necessary to open a PR in dipy for my proposal to get accepted or a link to my own repository is fine? [https://github.com/sahmed95/dipy/tree/ivim/dipy/ivim]

— You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub https://github.com/etpeterson/IVIM_fitting/issues/1#issuecomment-201014460

[quantshah]

Ok. Thanks.

[etpeterson]

@sahmed95 that's great thinking. Individual responses are below.

1. Yes, b>200 are nearly linear because the influence of D* is becoming fairly low at that point which removes that term and you get the D term only. That's why I was doing a linear fit to those values as a first step and this works fairly well to determine D at least (we're just omitting f at that point too). So I'd consider the b>200 side to be fairly easy. The tricky part is the b<200 because that's nonlinear and has 3 variables. Maybe a 2nd or 3rd order approximation would work there? If it at least approximates it well for b<200 that would be good enough!
2. I'm not sure any tissue would be homogeneous enough to approximate like that. I think the easiest approach would be along the lines of bullet point 1. But it would be nice to use neighboring voxels to help improve the fitting. That's much more complicated though.
3. Velocity would be great, but actually the typical perfusion number reported is ml/100g/min (I know it's strange units) but if we have a flow rate, flow fraction, voxel volume and a few assumptions about tissue density we could probably approximate it. Regarding velocity, can you get to this? http://web.stanford.edu/class/rad229/Notes/5b-Diffusion.pdf On page 11 there's the Einstein Relation where s^2=2Dt. We're fitting for D (actually D* is our flow variable) and t is the diffusion time of our acquisition so we can calculate s=sqrt(2Dt) which I believe tells us the average distance moved for our diffusion time. From there we should be able to calculate a velocity and we can see how that compares to known capillary flow.

Great thoughts and I hope that made sense.

[quantshah]

@etpeterson Makes sense but I will read the whole thing again next weekend after my mid semester tests. The class notes were really helpful.

I just wanted to know if adding functionality for velocity would benefit users. Once we get the parameters, it shouldn't be a problem. If it is helpful and this model can help to visualize blood flow, I guess I should look into that.

Thanks.

[etpeterson]

Take your time and let me know if you have questions, and I'll hopefully be able to answer them or at least find something that does.

Regarding the velocity, yes if we can get that then it would be good but the perfusion in ml/100mg/min would be better (maybe that's even the ideal result). But I myself am not totally convinced that we can convert from diffusion to perfusion so easily. I'll have to think about that some more.

[quantshah]

Hi @etpeterson. I opened a PR and will continue working on this over here