Calculation of FOD lobe integrals as used in SIFT/SIFT2

1

Hi,

I’m looking to compute the error between the FOD lobe integrals and the corresponding track densities in each voxel after applying SIFT/SIFT2. Would fod2fixel with the -afd option be suitable to compute the lobe integrals? If so how can I then compute the corresponding track densities in these voxel? I understand that the track densities in a voxel are the sum of streamline lengths in the case of SIFT and weighted streamline length in case of SIFT2. Is there a way to obtain the length of each streamline that traverses a given voxel?

Any help is greatly appreciated.

Thanks,
Varsha

2

Hi Varsha,

Easy answer:

I’m looking to compute the error between the FOD lobe integrals and the corresponding track densities in each voxel after applying SIFT/SIFT2.

Both tcksift and tcksift2 have an option -output_debug that will create both voxel-wise and fixel-wise output images containing this information. This may well be exactly what you’re looking for without necessitating any further gymnastics.

Hard answer:

Would fod2fixel with the -afd option be suitable to compute the lobe integrals?

Yes, this is the mechanism that we nowadays use as our quantification of ‘AFD’.

If so how can I then compute the corresponding track densities in these voxel? I understand that the track densities in a voxel are the sum of streamline lengths in the case of SIFT and weighted streamline length in case of SIFT2. Is there a way to obtain the length of each streamline that traverses a given voxel?

So there’s multiple aspects here:

  • You can quantify the sum of streamline lengths rather than the streamline count in each voxel by using the -precise option in tckmap. This activates the ‘precise’ streamline-voxel mapping mechanism described in the SIFT paper appendix, and for each voxel, sums the arc lengths of streamlines within that voxel.

  • tckmap accepts the per-streamline weights calculated by SIFT2 using the -tck_weights_in option. This applies both with and without the -precise option:

    • Using both -precise and -tck_weights_in, the voxel value will be: sum of (length * weight). Since the product of a length and a cross-sectional area is a volume, this is a sum of volumes.
    • Using -tck_weights_in alone, the voxel value will be the sum of streamline weights of those streamlines traversing the voxel.

  • The -output_debug option will give you fixel-wise measures of both fibre density and track density, which is more sensitive & specific than voxel-wise measures. Manually obtaining a fixel-wise measure of track density outside of SIFT & related commands will be possible in the future, but there’s still work for me to do on that code branch.

Cheers
Rob

3

Hi Rob,

Thank you for the detailed explanation. I shall try these out!

Regards,
Varsha

4

Hi Rob,

I see that tcksift and tcksift2 produce a lot of files with the -output_debug option. I have the following questions:

  1. What do the prefixes before and after signify?
  2. Which file corresponds to the voxel/fixel wise FOD lobe integral?
  3. There seems to be a file named after_diff.mif (after_diff_fixel.msf). Is this what measures the error between the FOD lobe integral and the underlying track density? If so, can I know how this measure was computed?
  4. Is there a way to read .msf files in MATLAB, keeping the fixel information intact and obtaining a mapping between the length of a streamline and the corresponding FOD lobe (given a voxel)?

Regards,
Varsha

5

  1. Before & after the application of the algorithm, whether SIFT or SIFT2. So “before” corresponds to the tractogram as it is fed in to the command, whereas “after” corresponds to the reconstruction either containing only the subset of tracks retained by SIFT, or incorporating the streamline weights calculated by SIFT2.

  2. The “target” reconstruction corresponds to the FOD lobe integrals that the algorithms manipulate the streamlines densities towards. (Not my best naming I agree)

  3. Yes, the “diff” or difference image is simply the difference between the fibre density and streamlines density: (u.TD - FD), where u = proportionality coefficient, TD = track density, FD = fibre density.

  4. I recall a user on here having written MatLab code to import .msf files; I can’t find it though. If you can wait a little while, the upcoming MRtrix3 tag update will completely change how fixel information is stored on disk, which will make it easier for users like yourself who want to access and manipulate such information directly.

Your last point (“obtaining a mapping between the length of a streamline and the corresponding FOD lobe (given a voxel)”) may require more clarity to get an accurate answer, as there’s a few different ways to interpret exactly what you’re looking for here.

Cheers
Rob

6

Hi Rob,

Thank you for the clarifications. As for the last point I made, I want to know if given a voxel and a streamline through that voxel, is it possible to find out which fixel (or FOD lobe) that streamline belongs to? This would particularly be useful if one needs to do post-hoc manipulations to the track densities attributed to a particular FOD lobe. It should be possible if there is a way to import the fixel information into a format that is amenable to such manipulations, which is why I wanted to know if .msf files could be imported into MATLAB.

Thanks,
Varsha

7

OK, I thought that might have been the case.

The reason I asked is that there’s two aspects to that question: not only do you need to be able to read fixel data, but you need to be able to determine the set of fixels traversed by any particular streamline (or equivalently the set of streamlines traversing any particular fixel). While the former is a matter of interface (which, given the timing, I wouldn’t even try right now; .msf is just about to become obsolete), the latter is more difficult, since currently commands such as tcksift will only give you the total streamlines density in each fixel, not the set of fixels traversed by each streamline. You could theoretically get at this information by using tckedit to isolate individual streamlines, feeding them into e.g. tcksift and analysing the subsequent fixel images, it’d be very clumsy.

I have a development branch hidden away that, when completed, will give you the information you’re looking for through the tcksample command. Alternatively, sometimes in such experiments it turns out to be much easier and cleaner to implement the actual processing within MRtrix3 itself, rather than trying to export the data & manipulate it externally; you’re welcome to PM me with details if you think this might be the case.

Rob

8

  1. If I do run tcksift on say a single streamline, I assume that the subsequent fixel images contain information (track and fiber densities, say) about the the streamline under consideration alone. Also, the resultant track densities cannot be used for further computation as the weights must be grossly overestimated due to the fact that we’re using a single streamline to reconstruct the FOD lobe integrals. However, the true track densities can be obtained using tckmap with the -precise option. Is that right?

  2. With regards to the voxel2fixel command, I understand that given an image, the command maps the scalar value in a voxel to all the fixels in that voxel. What exactly do you mean my “map”? If I had to do this to the track density fixel image, would it overwrite the fixelwise track densities or would it do some other operation (such as a multiplication with the track densities)?

Many thanks,
Varsha

9

  1. Yes, the main thing this would give you is the streamlines density per fixel, given a single streamline as input. You wouldn’t actually run SIFT /SIFT2 on that data, it’d be a trick purely to get at those streamlines densities. But on second thought, this would be even clumsier than I initially thought, since the FOD segmentation algorithm would have to be re-run on the input image every time you used a different individual streamline as input. So probably not even worth considering.
    You can get at the track densities using tckmap -precise, but this will only give you voxel-wise streamlines densities, rather than fixel-wise.

  2. With voxel2fixel, a new sparse output image is generated where the number of fixels in each voxel, and their directions and sizes, are defined by the input fixel image, but the quantitative value in each fixel is defined by the 3D input image voxel intensities. It won’t manipulate your input images in any way, or perform any kind of hidden mathematical operations that aren’t expected - just like we try to uphold with all MRtrix3 commands.

Cheers
Rob

10

Hi Rob,

Thank you for the clarifications. I think I have a better understanding of how this works now.

I also wanted to know how one could isolate individual streamlines using tckedit. I tried writing them into ascii files and subsequently converting them to .tck but a lot of the header information is lost, in particular the step size, something that seems to be necessary for tcksift2.

Thanks,
Varsha

11

There’s no need to be converting anything to ASCII: tckedit takes as input a tracks file, and outputs a track file. The difference between input and output will depend on the options you provide. So for instance, if you wanted to isolate an individual streamline, you would use -number 1 (one streamline in output) and -skip # (where # is the number of input streamlines to skip before apply any other selection criteria).

12

Thank you for your response.

I have a query regarding SIFT2. I’ve been using the new fixel format lately (I must say these are really nice to work with) and have noticed that there exist certain voxels in my data where there are no fixel elements but the track density in these voxels (as obtained using tckmap with the -precise option) is non-zero. How is it possible to not have any fixels in a voxel and yet have streamlines passing through them? I have noticed though that these track densities are very small relative to those in other voxels. Is there perhaps some sort of criteria being applied to exclude fixels with very small track densities attributed to them? I’m not specifying any such options from my end so I’m not sure what is happening here.

Regards,
Varsha

13

So there’s a couple of reasons why such a discrepancy can arise:

  • The default tracking algorithm (iFOD2) doesn’t operate on segmented fixel data, but raw FOD data, upon which interpolation is performed. If a tracking algorithm operated directly on segmented fixels, and did not perform any form of sub-voxel interpolation (as is typically the case when tracking using discrete mixture models with crossing fibres; see e.g. tckgen -algorithm fact), then yes, one would expect that any voxel with zero fixels would correspondingly contain no track density. But the presence of interpolation, and the potential for differing thresholds between iFOD2 and the segmentation algorithm used to identify fixels, means that a streamline within a voxel with zero fixels may still be able to propagate.

  • Even if a tracking algorithm such as FACT were used, if the integration is performed as first-order, then it is likely that the final streamline step will traverse from a non-empty voxel to an empty one, but that final streamline vertex in the empty voxel will be retained; hence there will be non-zero streamline density in that voxel upon performing the streamlines mapping operation.

There is a mechanism in SIFT2 for excluding fixels with very little streamlines density; however this should only exclude those fixels from contributing to the optimisation rather than ‘removing’ them; so I don’t think this is the effect you’re observing.

14

Hi Rob,

I want to know if there is a way to know what the value of the proportionality constant mu (as mentioned in the papers) is after running SIFT/SIFT2. You mentioned earlier that the image after_diff.mif is an image of the difference (u.TD - FD). I’m guessing the u here is exactly what I’m looking for. I worked my way backward, using the TD and FD from after_tdi.mif and after_target.mif to compute mu. But this value always turns out to be 1. Why is it so?

Additionally, the file named after_scatterplot.csv has four fields in it namely, Fiber density, Track density scaled, Track density unscaled and Weight. Is Track density scaled nothing but mu.TD where TD is the same as Track density unscaled? If so, I should be able to get at mu just by a simple operation. If this is indeed the case, the mu that I get from after_scatterplot.csv should match with what I get from after_diff.mif. But this doesn’t happen. Am I missing something here?

Also, what does the field Weight correspond to?

Thank you for your help.

Regards,
Varsha

15

I want to know if there is a way to know what the value of the proportionality constant mu (as mentioned in the papers) is after running SIFT/SIFT2

If you’ve run tcksift, this will be reported as a key-value entry “sift_mu” in the track file’s header, visible through tckinfo. For tcksift2, you don’t get this since no track file is generated (in the coming update there’s an additional command-line option to write this value to a text file, which can come in handy). The value does however get reported if you run either command with the -info flag (both before and after filtering in the case of tcksift, and just once for tcksift2 since it remains fixed).

So what I’ve done in the past is run tcksift with the -info and -nofilter flags: That’ll redo the FOD segmentation and track mapping, print the value of the proportionality coefficient, but won’t actually remove any streamlines.

I worked my way backward, using the TD and FD from after_tdi.mif and after_target.mif to compute mu. But this value always turns out to be 1. Why is it so?

Because after_tdi.mif is pre-scaled by u. :stuck_out_tongue:

Is “Track density scaled” nothing but mu.TD where TD is the same as “Track density unscaled”? If so, I should be able to get at mu just by a simple operation.

Correct - if you’ve generated that file already.

If this is indeed the case, the mu that I get from after_scatterplot.csv should match with what I get from after_diff.mif. But this doesn’t happen. Am I missing something here?

Suspect this is a continuation of the point above?

Also, what does the field “Weight” correspond to?

This was expressed as PM in the manuscript. I called it a “processing mask” early on and stuck with it, but in retrospect it’s a bad name because the influence of these data is non-binary. Different fixels contribute more or less toward the model on a continuous scale depending on the (estimated from 5TT image) partial volume contamination.