dipy versus mrtrix CSD results

[ebrahimebrahim]

I am having trouble getting CSD to give reasonable results in dipy. Any thoughts would be much appreciated!

I am just using ConstrainedSphericalDeconvModel, following the examples in the dipy docs:

# (load dwi, WM mask, gradient table, response function)
csd_model = ConstrainedSphericalDeconvModel(gtab, response_mrtrix, sh_order=8)
csd_fit = csd_model.fit(data, mask=mask_data)
csd_shm_coeff = csd_fit.shm_coeff
save_nifti('output.nii.gz', csd_shm_coeff, affine, img.header) # (where affine and img.header are taken from the original dwi)

When I use mrtrix's mrview tool and zoom in on the corpus callosum I see these funny shaped and misalgined FODs that have been generated by dipy ConstrainedSphericalDeconvModel:

If I use mrtrix for CSD rather than dipy, then the image looks much more plausible:

Both images were made from the same region of the same underlying DWIs and using the same response function.

Any thoughts what could be wrong with my approach?

EDIT:

I tried also the more dipy native method of visualization:

from dipy.data import default_sphere
from dipy.viz import window, actor
from dipy.sims.voxel import single_tensor_odf
csd_odf = csd_fit.odf(default_sphere)
fodf_spheres = actor.odf_slicer(csd_odf, sphere=default_sphere, scale=0.9,
                                norm=False, colormap='plasma')
scene = window.Scene()
scene.add(fodf_spheres)
window.snapshot(scene, fname="scene.png", size=(3000, 3000))

Below we see the same slice, the same FOD image, and roughly the same region. DIPY visual on the left and mrtrix's mrview on the right:

They are quite different.
This suggests that the issue I'm facing is not necessarily a discrepancy between dipy CSD and mrtrix CSD, because there seems to be a problem at the level of visualization. Saving the array csd_fit.shm_coeff generated by dipy seems to generate something that mrview is not interpreting correctly. It's still a mystery why. As @skoudoro suggested it seems to be a likely basis mismatch, but in my replies below I don't see how there could be a mismatch.

[skoudoro]

Hi @ebrahimebrahim,

Multiple possible reasons that comes in my mind:

• the basis: you can look at this document - https://dipy.org/documentation/1.7.0/theory/sh_basis/ - and check what is the basis use in each software. for sure the default one is not the same in DIPY and MRTRIX
• shm_coeff are they saved in the same order as mrtrix ? otherwise, the viewer will show the wrong thing. For example, the tensor format for Dipy's order is [Dxx, Dxy, Dyy, Dxz, Dyz, Dzz]. However, MRTRIX's order is : [Dxx, Dyy, Dzz, Dxy, Dxz, Dyz] which will display the wrong result in a viewer without a conversion.

Let us know when you check those point.

[ebrahimebrahim]

Hmm I just realized something. While I verified that the spherical harmonic basis order is the same between dipy and mrtrix, I did not verify that the actual basis functions are identical. I noticed here there is a default of legacy=True that isn't controlled from higher level callers... I wonder if this is makes the basis functions differ. Let me look more into this and report back.

[ebrahimebrahim]

Sources of disagreement in SH basis functions are explained here:

• Sometimes an extra factor of $\sqrt{2}$ is included to make the basis orthonormal
• Let $Y_l^m$ denote the complex spherical harmonics. Sometimes, when $m<0$, the function $(l,m)^\text{th}$ real basis function is given in terms of $\text{Im}[Y_l^{|m|}]$ while sometimes it is given in terms of $\text{Im}[Y_l^{m}]$.
• (EDIT There is also the choice of whether to take the real part at $m>0$ and imaginary part at $m<0$ or whether to do the opposite. I had missed this but it turned out to be the critical choice here-- see below.)

Now let's look at what choices are made in mrtrix versus dipy:

• In mrtrix we can see here that
  • the $\sqrt{2}$ is included, and
  • we use $\text{Im}[Y_l^{|m|}]$ when $m<0$
• In dipy
  • the multiplication by $\sqrt{2}$ is included in the basis used by dipy as we can see here, and
  • because legacy=True is forced in the use of real_sh_descoteaux_from_index when doing dipy CSD, we can see here that dipy uses $\text{Im}[Y_l^{|m|}]$ when $m<0$

And so the basis functions seem to be the identical.
Could they be using different definitions of the complex spherical harmonic functions $Y_l^m$?
No, because dipy uses scipy.special.sph_harm whose documentation page shows the same exact formula as the mrtrix documentation just mentioned.

So I am stumped 😪
Why does the FOD image from dipy look wrong when viewed with mrview...

[ebrahimebrahim]

After another round of debugging, I finally caught why the two bases are not the same.

The difference is in the m>0 here.
If it were m<0 then it would exactly match the mrtrix basis.

[skoudoro]

I am not sure to understand, because you are pointing the descoteaux basis and not the tournier basis here.

mrtrix basis should be equal to dipy + tournier basis + legacy=False if I am not wrong. in all the other cases, it should be different.

[ebrahimebrahim]

mrtrix basis should be equal to dipy + tournier basis + legacy=False

Exactly right. I will post an answer to that effect soon and mark this question resolved.

[arokem]

Thanks for all the debugging! You know what would be really useful? A script that converts between the two software libraries basis sets. Would you like to contribute that to DIPY? It seems like you know all the ins and outs of this issue at this point, so you'd be really well-prepared to do that.

…

On Sun, Nov 5, 2023 at 9:53 AM Ebrahim Ebrahim ***@***.***> wrote: After another round of debugging, I finally caught why the two bases are *not* the same. The difference is in the m>0 here <https://github.com/dipy/dipy/blob/618b39525ad5276f774b9920da983920d59cf340/dipy/reconst/shm.py#L357> . If it were m<0 then it would exactly match the mrtrix basis <https://mrtrix.readthedocs.io/en/latest/concepts/spherical_harmonics.html#formulation-used-in-mrtrix3> . — Reply to this email directly, view it on GitHub <#2959 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAA46NUOP5WDIIXGV3AWNOTYC7AB7AVCNFSM6AAAAAA6U7OHHCVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM3TIOBQGYZTK> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[ebrahimebrahim]

I was thinking the same 😀

I'll include a snippet here that does the conversion once I figure it out. Then I'll see if there's a good place for it to fit into DIPY.

[ebrahimebrahim]

I posted a conversion snippet as an answer to this question. Where do you think this would best fit into DIPY? A standalone script, or a utility function in dipy.reconst.shm, or somewhere else?

[arokem]

I think the most useful would be to add a function in reconst/shm.py and then use this function in a workflow to provide a cli for users who prefer the command line. I think the function code will also provide a good basis for solid unit testing.

[arokem]

The cli could potentially go into workflows/io.py with the input being a file with output from one software and the output compatible with the other software (might be interesting to enable mif io for full interoperability?)

[ebrahimebrahim]

PR here!

[ebrahimebrahim]

As suggested by @skoudoro it was a basis mismatch.

Difference in spherical harmonic bases

The basis used by dipy is the first (legacy) descoteaux07 basis described here, with code here.

The basis used by mrtrix can be found here. It exists in dipy as dipy.reconst.shm.real_sh_tournier with legacy=False.

However, dipy (v1.7.0 at the time of writing) is locked into the legacy descoteaux07 basis:

• In SphHarmModel, basis function values are sampled from the legacy descoteaux07 basis. This means that when we use dipy to visualize FODs based on spherical harmonic coefficients, it will always assume we were using this basis.
• In ConstrainedSphericalDeconvModel, the legacy descoteaux07 basis is forced to be used.

Converting between dipy and mrtrix representations

Following the convention in the dipy and mrtrix documentation we write $Y^m_l$ for the complex spherical harmonics and $Y_{lm}$ for various versions of real spherical harmonics. Let's say $Y^M_{lm}$ for mrtrix and $Y^D_{lm}$ for dipy. The mrtrix basis is given as

$$Y^M_{lm} = \left\{ \begin{array}{ll} \sqrt{2}\ \textrm{Im}[Y^{-m}_l] & \text{if $m<0$}\\\ Y^0_l & \text{if $m=0$}\\\ \sqrt{2}\ \textrm{Re}[Y^{m}_l] & \text{if $m>0$} \end{array} \right.$$

and the dipy basis is given as

$$Y^D_{lm} = \left\{ \begin{array}{ll} \sqrt{2}\ \textrm{Re}[Y^{-m}_l] & \text{if $m<0$}\\\ Y^0_l & \text{if $m=0$}\\\ \sqrt{2}\ \textrm{Im}[Y^{m}_l] & \text{if $m>0$}. \end{array} \right.$$

The relation could be viewed as a simple re-ordering of basis functions:

$$Y^M_{l,m} = Y^D_{l,-m}.$$

To be clear: the ordering of the $l,m$ indices used by dipy and mrtrix is the same. It's just that the basis functions are defined differently in such a way that the difference can be viewed as an effective change in ordering.

Here is a function that provides the re-ordering as a permutation (a list of indices):

from dipy.reconst.shm import sph_harm_ind_list

def get_dipy_to_mrtrix_permutation(sh_order):
    m,l = sph_harm_ind_list(sh_order)
    basis_indices = list(zip(l,m)) # dipy basis ordering
    basis_indices_permuted = list(zip(l,-m)) # mrtrix basis ordering
    permutation = [basis_indices.index(basis_indices_permuted[i]) for i in range(len(basis_indices))]
    return permutation

The CSD snippet in my original question can then be modified as follows to save the CSD results in a format that mrtrix mrview will interpret correctly:

# (load dwi, WM mask, gradient table, response function)
csd_model = ConstrainedSphericalDeconvModel(gtab, response_mrtrix, sh_order=8)
csd_fit = csd_model.fit(data, mask=mask_data)
csd_shm_coeff = csd_fit.shm_coeff
csd_shm_coeff_mrtrix = csd_shm_coeff[:,:,:,get_dipy_to_mrtrix_permutation(8)]
save_nifti('output.nii.gz', csd_shm_coeff_mrtrix, affine, img.header) # (where affine and img.header are taken from the original dwi)

Since the change of basis is idempotent (it's just a permutation swapping $m\leftrightarrow -m$), the same transformation can be done to go in either direciton: dipy-->mrtrix, or mrtrix-->dipy.