ENH: Add gaussian process DWI signal representation notebook

[jhlegarreta]

Add gaussian process DWI signal representation notebook.

[jhlegarreta]

Depends on #188. @esavary If you did changes to the notebook, maybe you can do a PR to this branch with such changes, and I would merge it right away.

With the default kernels, I get this error plot, so the default kernels are not a good choice for the DWI signal:

Not sure if the signal not being normalized plays a role here.

[esavary]

I made a pull request against your branch with some suggestions. Essentially, it involves not fitting one voxel at a time. This ensures that the same hyperparameters are used for the kernel across all voxels. It also improves their estimation when done with more voxels (this is done inside the fit method). These changes result in a slight improvement. However, even with these changes, the optimization is still failing for the parameter k2__noise_level (from the WhiteNoise kernel). It hits the bounds, leading to a significant error, as shown in this plot (I used the default hyperparameters optimization from the fit method of sklearn not the custom one I coded):

Additionally, I tried fitting a larger cube, but it seems to be similar as with the SphericalKernel. See the image below showing the prediction:

An additional comment I have is that we might consider using the angles between the directions instead of bvecs directly, as suggested in Andersson (2015). However, I was unsure if that was only relevant when using the Spherical or Exponential kernel, so I did not include that change.

[esavary]

Depends on #188. @esavary If you did changes to the notebook, maybe you can do a PR to this branch with such changes, and I would merge it right away.

With the default kernels, I get this error plot, so the default kernels are not a good choice for the DWI signal:

Not sure if the signal not being normalized plays a role here.

Actually, I've just noticed something in your plot, I think the ground truth signal should not be zero (but not sure if it causes the problem). Could it be that the data you used for the test is empty except for the first gradient direction? ( I used this data for my plot: https://osf.io/download/6at98/)

[jhlegarreta]

Merged your PR into the branch.

Maybe we can now uncomment the orthoview calls and provide with the same values you used to generate the figures you pasted.

An additional comment I have is that we might consider using the angles between the directions instead of bvecs directly, as suggested in Andersson (2015). However, I was unsure if that was only relevant when using the Spherical or Exponential kernel, so I did not include that change.

In my understanding we want to predict on bvecs. We want to represent a given DWI volume using a GP, e.g. as if we were using another basis (e.g. SH) to represent it and be able to recreate it within some tolerance. The ISMRM paper/poster figure is helpful in this.

Actually, I've just noticed something in your plot, I think the ground truth signal should not be zero (but not sure if it causes the problem). Could it be that the data you used for the test is empty except for the first gradient direction? ( I used this data for my plot: https://osf.io/download/6at98/)

Maybe the volume I chose included some background voxels.

[jhlegarreta]

Maybe we should normalize the DWI signal (divide each voxel signal by the B0 signal in that location) prior to fitting the data. Not sure if this would have a large influence on the results, but maybe worthwhile trying.

[jhlegarreta]

Added a new notebook with a simulated signal using a multi-tensor model.

@arokem In an effort to get something like fig 4 in Andersson 2015 NIMG, I am getting the below figure for the signal I am simulating in docs/notebooks/dwi_simulated_gp.ipynb.

I would expects the diffusion signals to roughly describe the two crossing tensors, but what I get is a blob. Tensors are at 90 degrees, and the fODF I get from them shows so:

What am I missing? Thanks.

[arokem]

Sorry for the slowness. One quick thought is that Andersson et al. used a b=3,000 in their simulations and I see that you are using b=1,000. Another is that I suspect that the Andersson figures were generated with a transformation to diffusion distance, rather than in the original signal domain. Indeed, a crossing fiber signal would look like the linear combination of two "donuts" rather than the linear combination of two "cigars", because the signal looks like a donut (see e.g., figure 1 in my 2014 paper), and when you add these donuts together, especially in low b-values, you get something that more closely resembles a "blob". To get "cigars" you need to do something like this: https://github.com/arokem/osmosis/blob/72d848ddc172ce6a72baba4c0d4cf3675db1a381/osmosis/tensor.py#L233-L248

[jhlegarreta]

Removed the signal plot (#202 (comment)) for now so that it does not distract us as we progress with the kernels, etc. Will add it in a separate PR on top of this branch.

[oesteban]

@jhlegarreta can you commit with the cells' output so that we can use github to visualize the results? It's tough to look at the diff, but we workaround the problem of the HCP data.

[jhlegarreta]

@jhlegarreta can you commit with the cells' output so that we can use github to visualize the results? It's tough to look at the diff, but we workaround the problem of the HCP data.

Brought back the cell in the last push force:
https://github.com/nipreps/eddymotion/blob/40baac4d3225060bc83ab5daf95a026c10f115c9/docs/notebooks/dwi_simulated_gp.ipynb

Note that in PR #217 I had tried with the approach Ariel suggested but did not get the expected shape.

[jhlegarreta]

Shows two disks, so not quite right yet.

Installed ipympl to be able to interact with (rotate) the mpl figures.

[jhlegarreta]

Have not succeeded in getting this right.

[jhlegarreta]

The plot looks perfect because I am predicting on the very same data that I used for fitting.

[jhlegarreta]

A few notes:

• Have not investigated the optimization convergence error shown on the notebok (ABNORMAL_TERMINATION_IN_LNSRCH.): when training and predicting on all bvecs (except b0) the warning was being shown, but not the error.
• I think I have figured out how to properly create the surface plots of the diffusion signal (which could be used for the prediction) and have it in a separate piece of code. When we take the time to discuss about how confident we are about the GP (e.g. the differences wrt the previous implementation, etc.), and other issues (example organisation, data and code to generate the empirical covariance plot) we can bring that part.

[oesteban]

Is there anything I can do here before the TechMon meeting?

[jhlegarreta]

@oesteban I apologise for not having worked on this in the last few weeks; I expect to be able to have another look around the tech meeting date. Sorry.

The single point prediction looks to be working, but I am not sure about the source of the error message mentioned in #202 (comment), so I need to investigate that. I would like to pick other random b-vecs to see if the prediction is OK or whether the convergence message persists.

Additionally, I would like to investigate what has changed wrt the implementation that was not working. I need to prepare
multi-fiber settings, and add the predicted point visualization part. So all this is on my roof; if you feel like, you can run the notebook and see if my impression that it is working seems OK (despite the convergence message).

Probably we can talk about all these experiments/design (esp. when we will add experiments with acquired data) in our eventual meeting.