Hi all,
We are attempting to create a toolbox for predicting white matter damage as a result of stroke.
For this purpose we were looking for an MRTrix function that can estimate the amplitudes of the spherical harmonics of an SH image along the vectors acquired by the sh2peaks function. Just to be clear, the vectors are acquired from a different SH image (a group mean). Is there a way of doing that?
Kind regards,
Mathijs Raemaekers
Hi Mathijs,
I’m not sure I’ve completely understood what you’re after… Do you need to evaluate the amplitude of the FOD from its SH coefficients along directions provided in a separate image? The sh2peaks command accepts an image of direction via its -peaks option, but that is used to initialise the search so it finds the nearest peak to these initial directions, which I’m not sure is exactly what you’re after. If you need to evaluate these amplitudes along the exact directions provided, and these directiosn change on a per-voxel basis, then there isn’t a command to do that, I’m afraid. It would be relatively trivial to code up though if you were interested. Happy to help, though right now I’m knee deep in teaching-related activities, so it would have to wait a bit…
Cheers,
Donald.
I think what you’re looking for here is for the second argument to the sh2amp command to be generalised, so that rather than being a text file containing a fixed set of directions to be evaluated in all voxels, you could instead provide either a fixel directory or a 3-vectors image, specifying a unique set of directions to sample in each voxel. As @jdtournier said, this would require some coding; but that’s probably where I’d put the functionality.
Hi Donald, Robert,
Thanks for the clear answers and your willingness to help. I meant estimating the amplitudes along the exact directions of the sh2peak output on a per voxel basis. I think I can use the -peaks option from sh2peaks for now, but it would be great if this functionality could be implemented. Btw, do you maybe have a code example where the amplitudes are calculated based on the coefficients (e.g. in the visualization or somewhere else)? I would like to play around with it a bit.
Regards,
Mathijs
From your original description, I don’t think this is what you want. What that option will do is take the directions from the image provided via -peaks, and use those as the starting points when searching for locap peaks within the input image. So what you get will be the peak orientations and amplitudes of the input SH image, not the amplitudes of the SH image along the orientations you specified.
In your case, where you just have one set of SH coefficients and one or a small number of unit vectors on the (half-)sphere in which you want to evaluate the amplitude, this is the function you want.
Rob
Hi Donald, Robert,
Thanks again for you helpfulness. I’m sorry I’m a bit unclear. I will give you some background information that might help to clarify.
We are attempting to validate a toolbox that predicts the damage to white matter fibers based on a lesion volume. The toolbox works by checking in a group of control subjects (in MNI space) which streamlines of a whole-brain tractography pass through voxels marked in the lesion volume.
We are doing the validation of the toolbox in a group of patients of which we also have DWI-data besides the lesion volume. We are warping the streamlines that are predicted to be damaged back to the native space of the patients (still in TCK-format). Using tckmap we can check for changes in e.g. FA and MD (based on the patient’s DTI-data) in the voxels that are predicted by the toolbox to include damaged streamlines. However, ideally we would like to see if there are changes in the FODs along the direction of the streamlines that are predicted to be damaged.
What we are planning to do this by creating an FOD-image of the streamlines that are predicted to be damaged, establish the peak directions in this FOD-image, and check in the data of the patients for changes in FODs along these peak directions. I believe that this would require the functionality that I was looking for. However, I’m happy to hear other suggestions in how to solve this.
Regards,
Mathijs
Hi Mathijs,
I think we may have followed some red herrings due to your original description being interpreted on our part more specifically than would have been ideal.
If we defer for a moment what a “peak” is, how one might be either discovered or represented, and whether it comes from FOD data or something else, and just think about the phrase:
and check in the data of the patients for changes in FODs
, this becomes a much more general question about how one might isolate some subset of an FOD image for quantification, but where that isolation specifically needs to occur in both the spatial and orientation domains.
Expressing the question in this way avoids the ambiguity of what a “peak” is, which is I’m pretty sure the principal source of the confusion. In retrospect our documentation page on “fixels” (which admittedly still needs a good image or two to demonstrate the concepts) would benefit from additionally disambiguating between “fixel” and “peak”…
We explicitly define “peak” as a location at which the FOD amplitude is maximal1. It has a direction on the unit (half-)sphere, and the FOD has an amplitude along that direction. The way that we identify peaks in command sh2peaks is best demonstrated in Figure 1 of this manuscript.
While a fixel also has a direction, and quite commonly some measure of “magnitude” with it, it differs from a “peak” in a number of ways:
For a “peak”, the direction is by definition that in which the amplitude of the FOD is maximal. For a “fixel”, the direction is an estimate of the principal orientation of the fibres within the voxel that are ascribed to that fixel. This could theoretically come from any estimation mechanism. In the context of MRtrix3 command fod2fixel, we by default use the mean orientation of the FOD lobe that was segmented to produce that fixel; this may differ from the peak orientation by quite a bit in cases of complex FOD shapes.
For a “peak”, there is only one scalar value associated with it, being the amplitude of the FOD at that orientation where the amplitude is a local maximum. For “fixels”, we can ascribe any scalar value we wish to each fixel. In a typical processing pipeline, the first such quantitative vvalue that we would derive is the Fibre Density (FD). This is not the peak amplitude of the FOD, but the integral of the FOD lobe in the half-sphere. This is invariant to fibre orientation dispersion, whereas the peak amplitude is not.
The way that we identify fixels in command fod2fixel is described in Appendix 2 of this manuscript.
So, back to the confusion.
When @jdtournier and I interpret very literally a question like this:
What we are planning to do this by creating an FOD-image of the streamlines that are predicted to be damaged, establish the peak directions in this FOD-image, and check in the data of the patients for changes in FODs along these peak directions.
, what we would respond with would be something like:
Use tckmap with the -tod option to produce Track (not Fibre) Orientation Distributions from a set of streamlines;
Run sh2peaks to identify the orientations of the “peaks” (maxima) of these functions in each voxel;
Modify the sh2amp command to optionally accept as input an image containing a set of orientations in each voxel, and in each voxel, evaluate the amplitudes of the subject’s FOD image only in those orientations specified by the image produced by step 2.
So the FOD amplitudes measured as the output of the experiment are not the amplitudes of the subject’s own FOD function peaks; they are the amplitudes of the subject’s own FOD function in those directions defined by the Track Orientation Distributions.
However, what I think you actually want, I would propose expressing something like:
“To evaluate quantitative FOD properties where traversed by a set of streamlines.”
If this is indeed still consistent with your actual intent, then my response would instead be something like:
Map streamlines to fixels (that were produced from segmentation of the subject’s own FOD functions) using the tck2fixel function;
Apply some threshold to these per-fixel streamline counts using mrthreshold to produce a fixel mask;
Evaluate per-fixel quantitative values (e.g. Fibre Density) within those voxels included in the mask.
A re-interpretation of words, but a completely different recommendation; also closer to the conventional approach for such an experiment when dealing with voxel-wise quantitative values (map streamlines to voxels; apply a density threshold; extract voxel-wise values from within the mask).
I hope that gets us all back on the same page… 
Rob
1 Unfortunately what might make this even more confusing is the fact that the way the term “peaks” is used in the context of MRtrix3 command names is different again. We have our convention of frequently naming commands based on the form of the data input / output. There, we currently use “peaks” for commands where the nature of the data is a 4D image, with 3xN volumes where N is the maximum number of discrete orientations in each voxel, and each triplet of volumes encodes an orientation in 3D Euclidean space, optionally with some scalar parameter encoded based on the norm of the vector. We’ve had multiple discussions about this ambiguity, because there’s a discrepancy between what a peak of an orientation distribution function is, and the image data format just described in which such data have conventionally been encoded, which is more like “3-vectors”, but we’ve not justified changing the relevant command names to something better. The fixel directory format can be thought of as a more generalised solution to representing such data, and hence developments around such are of greater priority than how we process—or indeed name—what we hope will over time become a “legacy” format.
Thanks again for your elaborate answer. I was studying a bit to see if I could solve the issue using the approach that you suggested, i.e. through tck2fixel and mrthreshold, as you mention that this would be more conventional.
However, I’m not sure that I understand fully how this should be implemented. I assume I would first do fod2fixel with the patient’s FODs as input, and after that using tck2fixel with as input: (1) the streamlines of the toolbox-output and (2) the fixel-folder created with tck2fixel?
As I understand tck2fixel will establish the directions of the fibre bundles in each voxel of the patient data, and tck2fixel will establish the fibre densities along these directions for the toolbox-ouput. I have the feeling that this approach would evaluate the toolbox-results against the patient data, as the patient data establishes the directions of fibre bundles. Wouldn’t it make more sense to do this the other way around, so establishing the directions (fixels) on the toolbox output, and somehow mapping the results the patients to these fixels? Not sure if I’m making sense here (if this is even possible).
Also, and sorry for not mentioning this before, we have 2 DWI datasets for each patient (several months apart) and we want to measure longitudinal changes al. This would mean establishing fixels twice, which would give 2 sets of directions for each voxel, and thus complicate the comparison between the sessions.
Otherwise I would really appreciate the option with modifying the sh2amp command.
Kind regards,
Mathijs Raemaekers
I have the feeling that this approach would evaluate the toolbox-results against the patient data, as the patient data establishes the directions of fibre bundles. Wouldn’t it make more sense to do this the other way around, so establishing the directions (fixels) on the toolbox output, and somehow mapping the results the patients to these fixels? Not sure if I’m making sense here (if this is even possible).
I was thinking originally that it would be done against patient data, which would therefore involve transforming the streamlines from template to subject space as you had previously mentioned. However particularly given:
we have 2 DWI datasets for each patient (several months apart) and we want to measure longitudinal changes
, I would suggest that it would be preferable to go further along the standard FBA pipeline: transform patient FOD data to template space, perform FOD segmentation, fixel reorientation, and fixel correspondence to the template. That way, once you have projected the subject fixel data to the template fixels independently for each of the two sessions, you then have for each subject the same set of fixels between the two sessions. For a subject you define a fixel mask based on the output of the toolbox, and then you are sampling values from the same fixels for both timepoints.
Rob