WIP: fODF interpolation (help needed)

[AntoineTheb]

Quick description

Context

This PR attempts to implement the log-euclidean framework of ref. [2] for fODF interpolation. For context, it has been demonstrated that diffusion tensors lay on the Riemannian manifold of 3x3 SPD matrices S3++, instead of the euclidean space. Applying standard euclidean metrics (such as interpolation, for example) may lead to the generation of invalid, non SPD matrices. Instead, the log-Euclidean framework for diffusion tensors was proposed (ref. [1]) to map the tensors from S3++ to euclidean space where euclidean metrics may be used, and then back to S3++, guaranteeing the SPD nature of the new tensors. Forgoing this transformation may lead to undesired effects such as non SPD tensors, but also the "swelling" of interpolated tensors (fig. 1).

Some theory

Naturally, ODFs representing diffusion also live on a nonlinear differentiable manifold where the application of euclidean metrics may lead to invalid ODFs. To this end, ref. [2] propose the log-Euclidean framework and even suggest an interpolation method. Here is a brief summary of the log-euclidean framework:

The space of ODFs representing diffusion can be posed as probability density functions

$$ (1) \mathcal{P} = { p : S^2 \rightarrow \mathcal{R} | \forall s \in S^2, p(s) \geq 0; \int_{s \in S^2} p(s)ds = 1 } $$

In other words, ODFs are functions that live on a sphere and, since they represent the probability of diffusion for a given direction, should sum to one. Following ref. [3], $p(x)$ can be reparametrized as $\psi(x) = \sqrt{p(x)}$ s.t.

$$ (2) \Psi = { \psi : S^2 \rightarrow \mathcal{R} | \forall s \in S^2, \psi(s) \geq 0; \int_{s \in S^2} \psi^2(s)ds = 1 } $$

, the positive orthant of a unit Hilbert sphere where the L^2 metric is defined (ref. [3]). Because this is a well studied manifold, the geodesic, exp and log maps are well defined. As well all know, \Psi can be further represented using spherical harmonics, giving rise to the coordinate space

$$ (3) PS_k = { c , || c ||_2 = \sum_i^k c_i^2 = 1, \sum_i^k c_iB_i(s) \geq 0 \forall s \in S^2 } $$

The log and exp maps allow for projection from the riemannian manifold to the euclidean space.

[2,4] defines the log and exp maps as

$$ (4) log_u (c) = v_c = \frac{v_c - u cos \psi}{||v_c - u cos \psi||_2} \psi; \psi = cos^{-1} \sum_i u_i c_i$$

$$ (5) exp_u(v_c) = c = u cos(\psi) + \frac{v_c}{||v_c||_2} sin(\psi); \psi = ||v_c||_2$$

with $u = (1, 0, 0, ... 0) \in R^k$ the isotropic ODF. (yes, $\psi$ is confusingly defined differently in exp and log).

This PR

The authors of [2] did not provide a reference implementation nor do they provide any experiment against a "reference" euclidean baseline. This PR is an attempt at doing both.

For a replication of some of the experiments of [2], see suppl. [1]. As we can see, the GA goes up linearly when interpolating from an isotropic ODF to a single-fiber ODF using the log-euclidean framework. This makes sense, as the ODF is less and less isotropic. When comparing to simply using euclidean metrics, we can notice extra anisotropy during the process. Interestingly, this would mean that not using the log-euclidean framework would have the opposite effect than on diffusion tensor, ie a "slimming" effect.

Problems

Unfortunately, it does not seem to work using real images and several problems arise when dealing with actual MRI data. See the test "fodf_test.sh" in the GDrive files. Interpolating from a 128x128x128 fODF volume to 64x64x64, then back to $128^3$ causes the first coefficient to "swell up". Looking at the Geometric Anisotropy after the 128-64-128 and especially the GA error (GA after reconstruction - GA before) makes it evident that using the log-euclidean framework is worse than not using it.

Some parts of the framework are ambiguous when using real images and I am not sure how to actually deal with them. Notably, the conditions in (1) and (3) are not enforced in real life, ie the SH coefficients do not always (and actually rarely) add up to 1. This is especially true outside the WM/GM. In the PR, you can see I normalized the coefficients, kept the norm and "de-normalized" after the interpolation. Does it make sense to do that ?

Overall, this little side project is starting to take up way too much time and I'm not quite sure it is worth it to invest some more. Maybe someone (@mdesco ?) with a better math background than I can make sense of it all. Right now, my brain turned to mush and I'm not quite sure what to do.

Supplementary
[1]: Riemannian_metrics.pdf
[2]: GDrive

Figures:

[1]: "Swelling" effect when interpolating tensors, from ref. 1

References:
[1]: Arsigny, Vincent, et al. "Log‐Euclidean metrics for fast and simple calculus on diffusion tensors." Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine 56.2 (2006): 411-421.
[2]: Cheng, J., Ghosh, A., Jiang, T., & Deriche, R. (2009). A Riemannian framework for orientation distribution function computing. In International Conference on Medical Image Computing and Computer-Assisted Intervention (pp. 911-918). Berlin, Heidelberg: Springer Berlin Heidelberg.
[3]: Srivastava, A., Jermyn, I., & Joshi, S. (2007, June). Riemannian analysis of probability density functions with applications in vision. In 2007 IEEE Conference on Computer Vision and Pattern Recognition (pp. 1-8). IEEE.
[4]: Anctil-Robitaille, B., Théberge, A., Jodoin, P. M., Descoteaux, M., Desrosiers, C., & Lombaert, H. (2022). Manifold-aware synthesis of high-resolution diffusion from structural imaging. Frontiers in Neuroimaging, 1, 930496.

Type of change

Check the relevant options.

• Bug fix (non-breaking change which fixes an issue)
• New feature (non-breaking change which adds functionality)
• Breaking change (fix or feature that would cause existing functionality to not work as expected)
• This change requires a documentation update

Provide data, screenshots, command line to test (if relevant)

...

Checklist

• My code follows the style guidelines of this project (run autopep8)
• I added relevant citations to scripts, modules and functions docstrings and descriptions
• I have performed a self-review of my code
• I have commented my code, particularly in hard-to-understand areas
• I have made corresponding changes to the documentation
• My changes generate no new warnings
• I moved all functions from the script file (except the argparser and main) to scilpy modules
• I have added tests that prove my fix is effective or that my feature works
• New and existing unit tests pass locally with my changes

[pep8speaks]

Hello @AntoineTheb, Thank you for updating !

There are currently no PEP 8 issues detected in this Pull Request. Cheers! 🍻

Comment last updated at 2024-07-30 18:44:34 UTC

[codecov]

Codecov Report

Attention: Patch coverage is 22.03390% with 46 lines in your changes missing coverage. Please review.

Project coverage is 68.48%. Comparing base (460dac8) to head (faa4105).

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1017      +/-   ##
==========================================
- Coverage   68.61%   68.48%   -0.13%     
==========================================
  Files         426      426              
  Lines       21835    21890      +55     
  Branches     3267     3272       +5     
==========================================
+ Hits        14982    14992      +10     
- Misses       5578     5619      +41     
- Partials     1275     1279       +4     

Components	Coverage Δ
Scripts	69.90% <36.36%> (-0.03%)	⬇️
Library	66.43% <18.75%> (-0.27%)	⬇️

[AntoineTheb]

Comparing with mrtrix' mrtransform, there doesn't seem to be any difference in the results (before the modifications of this PR), using or not the options to modulate and reorient FODs. Therefore I will be closing this PR, (simply) reslicing fODFs seems to be a non-issue and work as-is.

[CHrlS98]

Comparing with mrtrix' mrtransform, there doesn't seem to be any difference in the results (before the modifications of this PR), using or not the options to modulate and reorient FODs. Therefore I will be closing this PR, (simply) reslicing fODFs seems to be a non-issue and work as-is.

Just to make sure I understand, you mean mrtrix already does S3++ interpolation? And that the option to reorient FODs simply does not do anything?

[AntoineTheb]

Comparing with mrtrix' mrtransform, there doesn't seem to be any difference in the results (before the modifications of this PR), using or not the options to modulate and reorient FODs. Therefore I will be closing this PR, (simply) reslicing fODFs seems to be a non-issue and work as-is.

Just to make sure I understand, you mean mrtrix already does S3++ interpolation? And that the option to reorient FODs simply does not do anything?

No, they just do something else entirely, which may not be as mathematically sound but empirically works:

• Raffelt, D.; Tournier, J.-D.; Crozier, S.; Connelly, A. & Salvado, O. Reorientation of fiber orientation distributions using apodized point spread functions. Magnetic Resonance in Medicine, 2012, 67, 844-855
• Raffelt, D.; Tournier, J.-D.; Rose, S.; Ridgway, G.R.; Henderson, R.; Crozier, S.; Salvado, O.; Connelly, A.; Apparent Fibre Density: a novel measure for the analysis of diffusion-weighted magnetic resonance images. NeuroImage, 2012, 15;59(4), 3976-94

Both things amount to making sure the FOD orientations and amplitudes are correct after registration. What I found is that for simple reslicing, both methods don't seem to impact the result much.

So, it may indeed be that mrtrix and scilpy produce the wrong result, but from what I can see, scil_fodf_resample.py and mrtransform (with and without fod modulation and reorientation) produce exactly the same result.