First spherical harmonic coefficient Y(0,0) meaning

1

Dear Team,

I have a question regarding the interpretation of the spherical harmonic coefficients as produced with dwi2fod in MRTrix. From my understanding, the shape of the FODs can be represented as a linear combination of the spherical harmonics. In other words, in the simplest case, the basis function Y(l,m) that most resembles the shape of the FOD will have a large coefficient. For example, I have seen this in ACPC-oriented volumes in the white matter corresponding to the corticospinal tract (Corona Radiata, Posterior Limb of the Internal Capsule, and Cerebral Peduncles) for the coefficients of the Y(2,0) basis function which just looks like a vertical dumbell. The FODs are indeed very much vertically oriented here, so a large coefficient in this area makes sense to me.

That being said, I am confused as to the coefficients of the Y(0,0) volume, which is simply a sphere. My intuition tells me that an area with a large coefficient for the sphere would represent an area with high isotropic diffusion (i.e. grey matter and ventricles). However, upon examination of the Y(0,0) volume from MRTrix (i.e. first volume of the spherical harmonic coefficients 4D volume), the contrasts seem to be completely opposite from my expectations. The white matter has a much higher coefficient than the grey matter and the ventricles which are close to being black.

Is this due to some normalization of the coefficients that I am not aware of?
In another thread, (Fod2metric), the first volume was referred to as the total apparent fiber density (AFD) map. But if the basis function for this coefficient is a sphere, why would this represent the fiber density?

My last comment may seem off-topic, but when I calculate the spherical harmonic coefficients with Dipy using the peaks_from_model function, the first volume is simply a single value brain of 0.2821’s which is literally the equation for Y(0,0) = 1/sqrt(4*pi) (http://nipy.sourceforge.net/dipy/reference/dipy.reconst.html). However, shouldn’t there be a different coefficient for every voxel weighting the contribution of the Y(0,0) basis function? I will contact the developers of Dipy separately to try to understand this behavior.

My ideal goal is to process data with Dipy to generate the SH coefficient volumes and use mrregister for registration; however, due to the different behaviors of MRTrix and Dipy, I think it’s important to understand how the SH coefficient volumes are generated before naively performing registration with these volumes.

Thank you for your time and explanation,
Eric

2 Likes
2

Hello Eric,

The FOD file format (sign convention) is described in the help of the command dwi2fod.

Just a quick note on SH as my experience with them is rather limited: The separation into radial and angular part is a notion of spherical coordinates. I am not aware of an intuitive description of the SH coefficient space but you can think of the Y(0,0) component as the spherically symmetric part (monopole moment) and of the other coefficients as the axially symmetric dipole, quadrupole, … moments that together describe an angular function. (Analogue to Fourier transform harmonics in Euclidean space).

The zeroth component, taken by itself, can not encode any directional information. However, if you take higher order harmonics into account, you add the possibility to encode angular information and then the Y(0,0) component describes the average value when integrated over the sphere.

Is this due to some normalization of the coefficients that I am not aware of?

As Dave said, the Y(0,0) coefficients are “normalised” to the diffusion weighted signal intensity. Did you run MSMT_csd? Intuitively speaking it makes sense that, if the spherical harmonics express the angular density of white matter (CSD with WM response) then one would expect to see little WM amplitude in the CSF.

Hope this helps (and was correct…), let us know if not!

You can use shview to get a feeling for spherical harmonics and this and this paper give good introductions to spherical harmonics.

3

Hi Eric,

I think that although your simplistic interpretation of ‘most resembles’ is useful, but it’s the Y(0,0) term where the breakdown of that interpretation is most severe, hence the confusion. I’ll see what more I can contribute given my own limited experience with SH manipulation / logic.

… in the simplest case, the basis function Y(l,m) that most resembles the shape of the FOD will have a large coefficient. For example, …

More-or-less correct. However, I’d encourage a couple of additional observations here:

  • The Y(0,0) term will most likely also have a ‘large’ value here; not necessarily any larger than any other area of WM, but nonetheless non-zero, and probably of a comparable magnitude to the Y(2,0) term.

  • If you look at the actual Y(2,0) term (either from a printed chart or using some shview trickery), it has both positive and negative values: a large positive peak along the Z-axis, but also a large negative ring around the XY plane.

If an FOD expressed in the SH basis contained only the Y(2,0) term, then it would include a substantial number of directions in which the fibre density was negative. Combination of this with a positive Y(0,0) term is necessary in order to produce an FOD profile that both has a peak along the Z-axis, and is non-negative.

Now although one might argue that this FOD does not ‘resemble’ a sphere, the fact that the function is non-negative and also has large positive values means that it does in fact ‘resemble’ a sphere, depending on what brand of glasses you’re wearing. All SH basis functions above l=0 have a spherical mean of zero, therefore addition of a sphere with a positive amplitude is necessary by design in order to express in SH coefficients a function that has a non-zero spherical mean.

That being said, I am confused as to the coefficients of the Y(0,0) volume, which is simply a sphere. My intuition tells me that an area with a large coefficient for the sphere would represent an area with high isotropic diffusion (i.e. grey matter and ventricles).

There’s a couple of potential mistakes here:

  • Thinking in terms of correlations rather than the addition of basis functions. Let’s for example consider one WM image voxel with a very large fibre density, where all fibres are coherently oriented in a single direction, and another image voxel that contains just a tiny speck of GM surrounded by air. The ‘ideal’ FODs, with no noise and using the normalisation interpretation in MRtrix3 (in particular, no normalisation to b=0 signal and no normalisation of ODF to unit integral), the WM voxel will have a gigantic ‘blob’ in the fibre direction, and the GM voxel will be a tiny sphere (‘a little bit of fibre density in all directions’). The Y(0,0) term will be much larger in the first voxel than the second, since the mean value of the FOD function is large. If however you were to calculate a correlation on S2 between each FOD and the Y(0,0) function (a sphere), the correlation coefficient would likely be larger for the second voxel, since it ‘looks like a sphere’; but this is influenced by the normalisation that is intrinsic to a correlation calculation.
  • Thinking in terms of diffusivity rather than fibre density. In a region with ‘high isotropic diffusion’, the diffusion-weighted signal will be isotropic but small. So although the resulting FOD will (ideally) be spherical, it will still be small, since the size of the FOD is directly proportional to the diffusion-weighted signal.

As @maxpietsch suggested, making some fake SH functions in a text editor and viewing them in shview can be useful.

My last comment may seem off-topic, but when I calculate the spherical harmonic coefficients with Dipy using the peaks_from_model function, the first volume is simply a single value brain of 0.2821’s which is literally the equation for Y(0,0) = 1/sqrt(4*pi) (http://nipy.sourceforge.net/dipy/reference/dipy.reconst.html). However, shouldn’t there be a different coefficient for every voxel weighting the contribution of the Y(0,0) basis function?

There’s a chance that this relates back to something I half-skipped before: normalisation. Some diffusion MRI researchers choose to normalise the ODF in each voxel independently to unit integral, such that the ODF represents a probability distribution function (PDF). If such a step were performed, this would produce the result you’re observing. We avoid this step because it introduces a number of effects that we consider undesirable (e.g. the ODF in a voxel with 10% WM and 90% CSF is just as large as that in a 100% WM voxel), and doesn’t really provide any substantial benefits to outweigh those.

Rob

1 Like
4

To further explain/understand why exactly the value of the coefficient that goes with Y(0,0) in this case happens to be the same as the amplitude of Y(0,0), there’s 2 things to consider:

  1. Y(0,0) is, as Rob explained as well, the only one of all the basis functions that has a non-zero integral (but it’s integral is not 1, just non-zero!). All the other SH functions integrate to zero over the spherical (angular) domain. In a signal processing context, Y(0,0) would be referred to as the DC term.

  2. As Rob mentioned too, Dipy chooses to normalise to unit integral.

To compute the integral of an FOD (or any other function), represented by a sum of weighted SH basis functions, it is sufficient to consider only the integral of Y(0,0) times the coefficient that goes with it. Since the basis is linear, we could consider the integral of each function times the coefficient that goes with it independently, and sum those integrals; but since all the other ones integrate to zero, their integrals will be zero too. So we’re indeed just after that one of Y(0,0) times the coefficient that goes with it.

The amplitude of Y(0,0) is indeed 1/sqrt(4pi), and so is the coefficient that goes with it after this normalisation. The amplitude of Y(0,0) times this coefficient is thus (1/sqrt(4pi)) * (1/sqrt(4pi)) = 1/(4pi). A function with that amplitude as a constant over the spherical (angular) domain, can easily be integrated by just multiplying this amplitude with the “size” of the domain. That size is actually an integral too; but in this case it’s simply the area of the surface of the sphere. Now this would be annoying if we worked with degrees, so radians are used here. Furthermore, it’s important in this context to make a choice on what (size of) sphere we use; so the convention dictates we use a unit sphere (i.e. with radius = 1). The surface area of a sphere in general is 4pir^2, so that one of a unit sphere is 4pi. The total spherical integral of our aforementioned constant function with amplitude 1/(4pi) is hence 1/(4pi) * 4pi = 1. So, indeed, normalising any function represented in the SH basis to unit integral, will lead to the coefficient going with Y(0,0) being 1/sqrt(4*pi), the same value as the amplitude of Y(0,0).

I hope this all makes sense. I still remember once coming into this domain without prior knowledge of signal processing, and basically asking myself the exact same questions you do now; so I can appreciate where those questions and thoughts are coming from. :wink: Coincidentally, I had to teach bachelor engineering students about traditional signal processing as part of my TA activities back then; so I actually figured both of these things out in parallel, only later on understanding that the theory was essentially the same. That was a very interesting experience, to be honest. :smiley: (especially taking into account that these bachelor students typically think of you as being an expert in whatever you’re teaching; they probably witnessed the lectures before I started looking into the course materials myself… :rolling_eyes: ah well, “good times”).

5

Thank you all for your extremely detailed explanations and the time you took writing them. Indeed, I can see how my intuition was distorted when it came to understanding the contribution of each SH basis function. I understand now that in order to counter the negative portion of the basis functions, we require a significant contribution of the Y(0,0) term even though the underlying water diffusion may be anisotropic. Thank you also for helping me understand the behaviour of Dipy …that was an added bonus!

I hope this thread will help current and future users of MRtrix.

Thanks again,
Eric

2 Likes