MRtrix3 / mrtrix3

MRtrix3 provides a set of tools to perform various advanced diffusion MRI analyses, including constrained spherical deconvolution (CSD), probabilistic tractography, track-density imaging, and apparent fibre density
http://www.mrtrix.org
Mozilla Public License 2.0
294 stars 180 forks source link

difference between scalar and fixel AFD values? #1611

Open maxpietsch opened 5 years ago

maxpietsch commented 5 years ago

The fod2fixel -afd command outputs the "AFD per fixel (integral of FOD lobe)". So when I feed an FOD image to fod2fixel with the -fmls_no_thresholds option and sum all fixel in a voxel, I'd expect to get the total AFD in that voxel.

Ignoring negative lobes, I'd expect this to be the same as the l=0 term of the FOD image. Dave agreed it to contain the "‘total AFD’ measure".

However, the fixel sum is sqrt(4 pi) larger than the l=0 term. So if I was to fixelate the l=0 term of the FOD image with voxel2fixel and extract the mean fixel value per voxel, I'd get a result that is off by that factor.

Why is the fixel-AFD scaled by sqrt(4*pi)? If that is intended, I think we should document that somewhere.

thijsdhollander commented 5 years ago

The fod2fixel -afd command outputs the "AFD per fixel (integral of FOD lobe)". So when I feed an FOD image to fod2fixel with the -fmls_no_thresholds option and sum all fixel in a voxel, I'd expect to get the total AFD in that voxel.

Yep, and that's indeed what you get (ignoring negative lobes and small numerical inaccuracies).

I'd expect this to be the same as the l=0 term of the FOD image.

Nope, of course not. The l=0 term is an SH coefficient.

Dave agreed it to contain the "‘total AFD’ measure".

Keep on reading... 2 posts below that one: http://community.mrtrix.org/t/fod2metric/184/14 . Dave's initial statement (the one you quote) was not intended to be so specific, I think (hence it put the thing in quotes). The thing though is that the constant (0.282...) is indeed constant, because it merely relates to the definition of the SH basis, so for an analysis it may not matter that much. It might start to matter when you're talking absolute magnitudes or absolute effect sizes, but even then you got to be relatively careful with assigning a meaningful interpretation to that absolute number. We normalise it (e.g. mtnormalise) after all. Most accurate way would be what I sometimes say/state as: the unit of AFD is "times your response function", or even just "your response function". So 1 unit of (total or other) AFD in a voxel means the presence of the same amount of signal as contained in exactly 1 response function. If that signal is also distributed in the exact same way angularly (and still has the same amount), your FOD is a delta function (again, with total integral = 1). In an SH basis, for example, that would the be an SH delta function with integral 1. Such a delta function (with integral 1) has a DC term of 0.282... or more exactly 1 / sqrt(4*pi) indeed. This is where that factor pops up.

So that explains why:

However, the fixel sum is sqrt(4 pi) larger than the l=0 term.

And thus...

So if I was to fixelate the l=0 term of the FOD image with voxel2fixel and extract the mean fixel value per voxel, I'd get a result that is off by that factor.

...that would indeed happen, but for what I think you intend to do there, you should've first re-scaled the l=0 term before feeding that into voxel2fixel.

So finally:

Why is the fixel-AFD scaled by sqrt(4*pi)? If that is intended, I think we should document that somewhere.

So it's not the fixel-AFD that is scaled at all. It's the SH coefficient that it scaled (relative to the integral), because it's an SH coefficient. That, I'd argue, is properly documented wherever we define the SH basis.

Lestropie commented 5 years ago

Agreed the factor is best interpreted as scaling the SH basis, not the fixel FDs: dwi2fod | fod2fixel on a voxel containing the response function should give an aPSF with integral 1.0, which makes sense. However:

I think we should document that somewhere.

That, I'd argue, is properly documented wherever we define the SH basis.

Actually I don't think this aspect is adequately documented. There are multiple conventions for normalising the SH basis; our documentation in commands that use SH data only describes the order of volumes, not the normalisation.

(E.g. We know that compared to Dipy we omit the (-1)m factor, but I don't actually know whether or not their normalisation is the same.)

jdtournier commented 5 years ago

OK, so it's all self-consistent as long as we agree that the convention is that fixels integrate to one, and SH coefficients for the corresponding FOD integrate to 1/√4π. I think this does deserve to be clearly documented somewhere (assuming it isn't already), since it has clearly caused some confusion here...

thijsdhollander commented 5 years ago

It's always self-consistent, by definition of SH basis; unless we don't agree with our current SH basis... :stuck_out_tongue_closed_eyes: Independent of SH basis, deconvolving a function with itself always results in a delta function with unit amplitude and integral.

@Lestropie brings up a good point though, that the actual SH basis itself (that MRtrix uses) is not provided in any documentation; I had overlooked that fact indeed.

An idea might be to make use of the https://mrtrix.readthedocs.io/en/latest/concepts/orthonormal_sh_basis.html page to formally define the SH basis we use? The current content of that page might just as well be removed (?), as that issue (compatibility with MRtrix 0.2) is becoming less relevant nowadays (and might even confuse most of our current user base...?). To take that even one step further, with release in mind, this might be the time to start to consider to start to remove constructs referring to the MRtrix 0.2 SH basis? At this point, backwards compatibility with MRtrix 0.2 is probably broken in other ways in any case; it might just be easier to openly break it entirely that state that as such. In this context, there's the shbasis command that I'd argue is starting to become somewhat redundant (and maybe even confusing). Are there other places in the code still where we allow or try to detect MRtrix 0.2 SH basis conventions? (I don't think so, but I'm not sure)

thijsdhollander commented 5 years ago

A relevant statement from that page I linked to above:

We hope that once we have feature completeness in MRtrix3, the old version will no longer be necessary, and therefore this will not be a problem.

Looks like we've reached that stage, I reckon.