[ENH] Compute fiber density and fiber spread from ODF using Bingham distributions

[CHrlS98]

As discussed with @villalonreina @RafaelNH @mdesco,
This is an attempt at migrating code from scilpy for fitting Bingham distributions to ODF and computing fiber metrics such as the fiber density and fiber spread [1]. The idea was that this could serve as a basis for GSoC 2023 Project 2: Generalized along tract analysis of fiber orientation dispersion.

The PR is WIP. It still misses unit tests and stuff. I would mostly like some feedback at this point:

1. What would be the preferred data structure for representing Bingham distributions and their derived metrics. Right now I simply use python lists of tuples, but maybe we would like to use a numpy array like we do for tensors... It seems clearer as tuples because each parameter of the distribution is a separate output.
2. Do we want to support arrays of SF instead of a single SF at this point? For instance, given some ODF image, we could output a "Bingham image" like it's done in our scilpy script here. Of course it would be cool, but maybe not needed in a first PR.
3. This does not feel like a reconst module. It's more of a voxel based analysis or something. Should this go in another module? Is it worth creating a new one?

Checklist

• Write unit tests.
• Support array of SF instead of single SF?

References
[1] Riffert TW, Schreiber J, Anwander A, Knösche TR. Beyond fractional anisotropy: Extraction of bundle-specific structural metrics from crossing fiber models. NeuroImage. 2014 Oct 15;100:176-91.

[codecov]

Codecov Report

Attention: Patch coverage is 81.48148% with 40 lines in your changes are missing coverage. Please review.

Project coverage is 83.60%. Comparing base (eec0ad8) to head (cb403e1).

❗ Current head cb403e1 differs from pull request most recent head 25ef2fc. Consider uploading reports for the commit 25ef2fc to get more accurate results

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2826      +/-   ##
==========================================
- Coverage   83.76%   83.60%   -0.16%     
==========================================
  Files         153      153              
  Lines       21340    21468     +128     
  Branches     3445     3455      +10     
==========================================
+ Hits        17876    17949      +73     
- Misses       2608     2646      +38     
- Partials      856      873      +17     

Files	Coverage Δ
dipy/viz/plotting.py	7.33% <7.69%> (+0.04%)	⬆️
dipy/direction/bingham.py	86.20% <86.20%> (ø)

... and 12 files with indirect coverage changes

[villalonreina]

Hi @CHrlS98, thanks for initiating this pull request. Regarding your questions:

1. We can keep it as lists and tuples for now, but at the end we must convert it to a numpy array to build the Nifti file for the full brain image.
2. We will add the feature to loop over the whole brain image, that is to do it over an array of ODFs (á la scilpy).
3. @RafaelNH and I discussed it and we think it should go into the direction module.

[pep8speaks]

Hello @CHrlS98, Thank you for updating !

Cheers ! There are no PEP8 issues in this Pull Request. 🍻

Comment last updated at 2024-05-13 13:49:54 UTC

[skoudoro]

Hello @CHrlS98,

What is the status of this work? Do you know when you will have time to progress on this?

ps: next release should be on beginning of March, it would be great to have this work

[CHrlS98]

Hi @skoudoro ,
It can be done before the next release. I will get back to it in the upcoming weeks.

[villalonreina]

Hi @skoudoro. @CHrlS98 and @RafaelNH and I have been working on this PR. I am trying to do a pull request to Charles' PR but Charles' Dipy fork does not show up in Dipy's Pull request interface. See screenshot attached. Do you know what might be the issue causing this?

[skoudoro]

Hi @villalonreina,

No idea, but when it does not show up, just go directly to his repo and it will show up: https://github.com/CHrlS98/dipy/pulls

[CHrlS98]

@villalonreina @RafaelNH
So I worked on the code and here are a few comments we might want to discuss in a future meeting.

• On the fiber density
  I changed the way the fiber density is estimated. Instead of creating a (theta, phi) grid with constant theta and phi, I use a DIPY repulsion sphere: this way, we are sure all points are at equal distance and cover the same surface area. The sphere can be passed as a parameter. I also added some tests for the new method.

• On the maximum separation angle
  I tried lowering it using the tests like we discussed. It would seem there is a difference between odf_to_bingham and _bingham_fit_peak because for the former I can decrease the angle to 35 while for the later it won't work below 45.

• odf_to_bingham and bingham_from_odf
  Both these methods exist, and from my understanding, bingham_from_odf takes a whole ODF volume as input whereas odf_to_bingham takes a single ODF as input. We might want to change the name because it is quite confusing.

[CHrlS98]

@villalonreina @RafaelNH

Notes from our last conversation. To do:

• T-design for default bingham_fiber_density()
• peaks_direction_nl or subdivide the sphere even more for _bingham_fit_peak and odf_to_bingham to give the same result in tests.
• always use multi-voxel methods for single and multi-voxel. Names bingham_to_sf, sh_to_bingham, sf_to_bingham.

[CHrlS98]

Hi @RafaelNH, is it possible t-design is slightly biased? I create my sphere this way (it doesn't matter if I subdivide it or not, I tried both) :

# default sphere is a subdivision of 45 points spherical t-design (11010 points)
verts = np.loadtxt(get_fnames('t-design'))
sphere = Sphere(xyz=verts).subdivide(4)

And then I added a test where I check if the integral of a Bingham distribution along the y-axis is equal to the integral along the x-axis:

# Rotating the Bingham distribution should not bias the FD estimate
xfits = [(f0_lobe1, k1, k2, axis0, axis1, axis2),
         (f0_lobe1, k1, k2, axis1, axis2, axis0)]
xpars = _convert_bingham_pars(xfits, 2)
xfd = bingham_fiber_density(xpars)

assert_almost_equal(xfd[0], xfd[1])

And sure it's close, but not even that close:

E           AssertionError: 
E           Arrays are not almost equal to 7 decimals
E            ACTUAL: 6.6747531900990253
E            DESIRED: 6.690992893100228

However when I use a unit icosahedron and subdivide it 5 times (10242 vertices), the test passes. I'll keep the unit icosahedron for now because it is better distributed than repulsion spheres, and add an option to control the level of subdivision.

Is the t-design a whole sphere or a hemisphere? Could it be the reason for this bias?

[CHrlS98]

@RafaelNH If I replace peak_dir by -peak_dir I can lower the max angle to lesser than 45 and the test comparing Mus and Mus_ref` will still pass.