How to compute AFD?

1

Hi there,

I am new to the MRtrix environment. I am interested in computing AFD maps for my data. What is the easiest way to do that? Is there a specific command for that? I looked into the documentation and I saw the afdconnectivity command. Is this the one I am supposed to use?

Thanks in advance for your help.

Julio Villalon
University of Southern California

2

Hi Julio,
Welcome to the new forum. Our first customer! No prizes sorry.

In our original work we computed the AFD as the FOD amplitude for each direction in a set of directions equally distributed on a sphere. You can do this with sh2amp. However, more recently we compute the integral of the AFD over a FOD ‘lobe’ (where each lobe represents a unique fibre population within a voxel, i.e. a fixel). This can be computed using fod2fixel, and you can visualise the *.msf file output with the vector plot tool in mrview.

Very soon we will release code for FOD-based registration, which will enable you to perform a whole-brain fixel-based analysis as described in the CFE paper.

A word of warning though, to ensure the AFD is quantitative and comparable across groups you need to make sure your images are intensity normalised (across subjects), and bias field corrected. You also need to perform the CSD using a single response function (we normally use the group average). I’ll add a tutorial on how to perform these steps in the coming weeks.

With respect to the afdconnectivity command. As explained here, it was originally created it as a quick tool for obtaining a ROI-based AFD measure, with the idea being that this would be a psuedo measure of the pathway’s ‘connectivity’. As explained, using SIFT is more appropriate if the images can be EPI-distortion corrected. While not as appropriate as SIFT from a theoretical point of view, the afdconnectivity command has been used in three clinical studies so far with plausible findings (see here, with the others under review). We hope to write this up as a technical note at some point. Note that when using afdconnnectivity, your images also need to be intensity normalised, bias corrected, you need to use a single response for the whole group.

Hope this clears things up,
Dave

3

Hey there @Julio_Villalon ,

Welcome indeed! :thumbsup:

If you’re really after a “classic” scalar map, you can check the first volume (the DC term of the SH representation) of the FOD file (i.e., what gets output by dwi2fod). To see it, just open up the FOD file directly in mrview. That contrast shows (up to a globally constant factor) the total AFD per voxel. But, as should be quite clear from @Dave’s comprehensive explanation, an added benefit of AFD is the angular dimension; which provides you with more specific AFD measures per fibre present in a voxel (i.e., what we refer to when we use the term “fixel”).

Note that, even when simply using such a “total” AFD map for any purposes, the need for bias field correction, intensity normalisation and a single response for all subjects for which you need comparable AFD maps, remains valid (and very important).

4

Hi @Dave and @ThijsDhollander,

Thanks a lot for your response. I have tried the commands per your advice and I have some further questions if you don’t mind.

So, yes, as Thijs said, my aim is to compute “classic” scalar maps for my scans and do population statistics. So it seems that I can compute two types of AFD scalar maps. One is a “dixel-based” AFD and the other one is a “fixel-based” AFD. Am I right?

  • The DC term of the SH representation (first volume of the FOD): that would be the “dixel-based” AFD?

  • The output of the fod2fixel is supposed to be the “fixel-based” AFD?

One more thing. After computing the Spherical Harmonics (SH) with dwi2fod, I am using them as an input to the fod2fixel function but it is giving me an error:

[ERROR] cannot create sparse image buffer for accessing non-sparse image

Am I missing a step in-between dwi2fodand fod2fixel?

I greatly appreciate your help.

Thanks!

Julio

5

Looks like we didn’t get all the messages across correctly yet
 :wink:

Let’s see if this overview works better:

  • If you’d name the DC term of the SH representation anything along that line of terminology, it should probably be the voxel-based AFD: it represents the total AFD for each voxel. It is the only one that can be fully represented by a “classical” scalar image, where one voxel comes with a single intensity. As a small side note, if you want the “real” total AFD (i.e., the integral of the FOD), you should divide the DC term by about 0.282
 Of course, this won’t change the contrast of the image (since it’s a constant factor for the entire image). So to obtain it in a single scalar map, you could use the following string of MRtrix commands (where fod.mif is your FOD image, as output by dwi2fod: <br><br> mrconvert fod.mif - -coord 3 0 | mrcalc - 0.282 -divide afd.mif <br><br> The units of afd.mif in this case could in principle be referred to as “times your single fibre response function (integral)”. Formulating things this way, makes clear the FOD from ©SD techniques expresses/models/
 your data in terms of a certain response function. Does this make sense? :slight_smile:

  • A dixel is a “direction in a voxel”. It’s a concept that was initially introduced to extend the 3D spatial domain with 2 extra angular dimensions (e.g., angles theta and phi that parametrise the angular domain). Just as voxels provide a discrete (gridded) sampling of the spatial domain, something that discretises the angular domain is needed for this purpose. A (very) dense grid of orientations was typically used (e.g., a set of 1000 or more directions). The AFD value per dixel in this case would be the FOD amplitude in a given voxel, along a given direction (multiplied by the angular “area” for the dixel). Obviously, this would be a very large file (basically 1000+ times the size of any of your regular (voxel-wise) scalar maps. Also, to compare AFD for corresponding dixels, the registration would need to be super accurate in the angular domain. Even the smallest deviations could cause huge discrepancies in values, since a single FOD lobe is super sparse in the angular domain. So we don’t really use dixels any more in practice. You can still do track-mapping (TDI) to dixel-space though; but I can’t see much practical use for that any more as well


  • A fixel is a “fibre population (with an orientation) in a voxel”. It’s a concept introduced to overcome the problems with dixels mentioned before, and AFD per fixel makes more sense from an anatomical point of view. A fixel is still something that lives in the angular domain for a given voxel, but the catch is that the fixel “domain” directly relates to the anatomy at hand: rather than having, e.g., 1000+ predefined dixels per voxel, you’ll just have a fixel for each fibre population per voxel. So for instance, in a voxel somewhere in the middle of the corpus callosum, each voxel should have a single fixel, whose orientation is aligned with the local fibre bundle passing through the voxel. The “total AFD” for that fixel in this specific case will/should be the same as the total AFD for the corresponding voxel; since there is only a single fixel in said voxel. So the real angular domain of the fixel is not limited to a single discrete orientation (that orientation could also be seen as the average orientation for the fixel if you will), but rather extends to the angular area/domain that the fibre occupies. In case the bundle bends a bit, or has some dispersion (but could still be seen as a single population), the angular distribution is represented in the FOD. The AFD for the fixel is the integral of the associated FOD lobe. In case of a crossing of fibre populations in a voxel, you will have an equal amount of fixels in the associated voxel. In principle, their combined total AFDs (FOD lobe integrals) should sum to the total AFD of the voxel (full FOD integral). The command fod2fixel will try to segment your FODs for each voxel in “lobes” (this is a tricky and surprisingly ambiguous task, mind you), and store the average orientation for each lobe as a fixel. Via the -afd option of fod2fixel, it will also store the associated FOD lobe integral (“total AFD”) for each fixel. Doing statistical analysis on fixel images requires a dedicated framework. We happen to have (the only) one of those lying around, after years and years of @Dave’s hard yakka. The public release of this gem (diamond) is soon upon us!

As to your error at hand: that one’s probably caused by you suggesting the wrong extension for the fixel file. We should fix this, since we’ve had heaps of users already running into this
 but what you should do is the following:

fod2fixel fod.mif -afd my_fixel_image.msf


it’s all about that .msf extension. That’s the (clearly still quite obscure) extension we use for our sparse image types; of which a fixel image is one. A redesign of how (sparse) fixels are represented in such a file is still due, so beware this might still change at some point.

I hope all of the above makes some (more) sense now
 :sunglasses: :sun_with_face: (we’re having a nice 33° day out here; I’m trying to enjoy it in smileys)

1 Like
6

I have been following this conversation closely and the detailed info is very useful. I am eagerly waiting for @Dave’s work on fixel-based statistics.

In the meantime, I was wondering if it is possible to translate the .tsf into voxel-wise maps. This would be done after segmenting the FOD with fod2fixel and supplying a track file to fixel2tsf. Now I would like to have a volume in which each voxel informed me of the average AFD of the fixel related to the track in question.

fixel2voxel does not seem to be a good match for this job, as it does not consider the track and will provide information on all fixels within each voxel. On the other hand, tsfinfo has an -ascii switch that nicely outputs the required scalar per each point in each streamline, as a series of txt files. This is handy, but over-represents voxels with several streamlines passing through.

I am imagining that the scalars computed by fixel2tsf come from interpolated FODs at each point of the streamlines, is this correct? Is the option of having an average track-specific fixel scalar possible?

Thanks.

7

Thanks @ThijsDhollander for you thorough explanation.

Regarding the “directions in a voxel” or dixels, I assume these are the ones that are computed with the dirgen command that needs to be run before running sh2amp. Does the number of directions change the AFD output? How many are good enough? @ThijsDhollander mentioned 1000, in the original paper I saw some examples done with 200. Any suggestions?

And one last question for @Dave. When computing the AFD for a whole population, Dave mentioned that an average response function should be used. What I am planning to do is to compute the response function for all the subjects in the population and then compute the average for each coefficient of the response function across the population. Is this the right approach?

Thanks again for your help!

Julio

8

Hi Guys,
Just at add to what Thijs said about dixels, this is what we previously performed analysis on in our 2012 AFD paper (although we did not call them dixels back then). We started using the word dixels instead of the word ‘direction’ since it eliminated some ambiguities when discussing neighbouring dixels in adjacent voxels.

@Julio_Villalon, yes we used 200 directions in that paper, and they were computed using dirgen. However, as Thijs has pointed out we no longer perform population analyses using dixels for a number of reasons:

  • Residual misalignments in orientation (even after FOD reorientation) can cause additional variance as well as false positives. In the 2012 paper we mitigated this somewhat by smoothing in the angular domain (truncating the SH series to lmax = 4). However in practice this is still an issue.
  • Differences in FOD lobe ‘spead’ between populations may be detected in differences in AFD (when the total FOD lobe integral is in fact the same).
  • Spatial smoothing is not ideal. In the 2012 paper we smoothed the dixel AFD using an anisotropic Gaussian kernel oriented along the dixel direction. However what you really want is smoothing along the fibre tract itself (see here), and minimise smoothing across different tracts.
  • With 200 directions we ended up with very large files. So processing large datasets was very time consuming.
  • The cluster-based statistics used in the 2012 paper were threshold based, and not tract-specific (unlike the CFE method).

So for the above reasons, I would recommend performing a fixel-based analysis instead. @lconcha, the statistical analysis command (fixelcfestats) is currently available in the updated_syntax branch, which we plan to merge with master very soon. You may also want to perform FOD registration to a population-specific template. For that you will need the registration tools that @maxpietsch and I have recently finished coding up. We will be merging it in the next couple of weeks after some more testing. I’ll post back here when we do.

@lconcha, to get a tract-specific AFD voxel map, you might want to try the afdconnectivity tool with the -afd_map option. This essentially uses the input streamlines to identify the fixels belonging to your tract-of-interest. It then computes the fixel-specific AFD (lobe integral) and outputs it as a 3D scalar map.

Hopefully this helps,
Cheers,
Dave

9

Thank you, @Dave, afdconnectivity -afd_map is exactly what I was looking for.

Looking forward for all the updates that you and @ThijsDhollander have commented!

10

Is it possible to put the afd of all different subjects in a template? fixelcorrespondence uses the wm template and the reoriented afd, but the output is one for each subject right?

11

Do you mean in a single file?

Yes, fixel correspondence will output 1 file, however this file will have the same number (and orientation) of fixels as the template (just different fixel-specific quantitative values, e.g AFD). So they are all in the same ‘space’ as the template.

As a side note (and hopefully not to confuse things), the next MRtrix release will use a new fixel image directory format, which is much more flexible and can store any number of quantitative values per fixel. This enables you to put the AFD for multiple subjects in the same ‘fixel directory’, which you can then feed into fixelcfestats, or visualise in mrview (and switch between patients).