Voxel coordinates of streamline endpoints

Hi,

Sorry to open up again this old thread, with what are probably rather basic questions (to which I could however not find an answer).

I ran seed-based tractography seeding within a subcortical region, and I am specifically interested in the voxel coordinates of the streamlines reaching the cortex.
I was wondering, how are the voxel coordinates of the two endpoints encoded in the .tck file? The read_mrtrix_tracks() gives me access to a bunch of 2x3 cells. How do I know which x,y,z coordinate belongs to the streamline’s endpoint and which to the streamline’s startpoint (if I am even allowed to call it like this. I anyway mean where the streamline was seeded).

The tckconvert -scanner2voxel command requires a reference image; does this need to be in the same space as the tractography’s space (e.g. T1 coregistered to DWI)?

I eventually would like to also translate the obtained voxel coordinates to an image. Is it then more convenient to simply use the tckmap -ends_only? Would there be a way to extract the coordinates of the voxels with the highest streamline density?

Cheers,
Sam

Hi Sam,

Moved your question to its own topic since it expands the scope of the original topic a little.

I was wondering, how are the voxel coordinates of the two endpoints encoded in the .tck file? The read_mrtrix_tracks() gives me access to a bunch of 2x3 cells. How do I know which x,y,z coordinate belongs to the streamline’s endpoint and which to the streamline’s startpoint (if I am even allowed to call it like this. I anyway mean where the streamline was seeded).

So I think I can infer from your description that:

  1. You are looking at the output of tckresample -endpoints, since “2x3 cells” implies that each streamline only has 2 vertices, each with a 3D coordinate;
  2. You have performed unidirectional tracking, as in this scenario, one of the streamline endpoints must be the seedpoint; if you had performed bidirectional tracking, this would not be the case.

In that case, the first row corresponds to the streamline seed point.

The tckconvert -scanner2voxel command requires a reference image; does this need to be in the same space as the tractography’s space (e.g. T1 coregistered to DWI)?

The reference image can be whatever image you want. The point of that command usage is that you want to know, for each streamline endpoint as defined in scanner space, where that corresponds to on a particular voxel grid; it’s that grid you are providing, and what that grid is depends on why you are doing that operation.

I think there may be some confusion about how this works given the question. For instance, if one has both a DWI series, and a T1-weighted image that has been coregistered to the DWI but not resampled onto the same voxel grid as the DWIs (i.e. the higher spatial resolution of the T1-weighted image is retained), then tckconvert -scanner2voxel would give different answers depending on whether you use the DWI or the T1w image as the template. These images can have their anatomical content overlapping in space, but not necessarily be sampled on the same voxel grid. Since the very operation of this command usage is to report the voxel location that corresponds to particular locations in space, the definition of those voxel locations matters.

I eventually would like to also translate the obtained voxel coordinates to an image.

The statement “translate voxel coordinates to an image” is unfortunately also a bit ambiguous. tckconvert -scanner2voxel will give you a voxel location for each streamline endpoint, but it’s unclear what is the “translation” to which you are referring.

Is it then more convenient to simply use the tckmap -ends_only?

If what you ultimately want is an image containing the density of streamlines terminations, then yes, it would likely be more convenient to use tckmap -ends_only (or just tckmap and feeding in the output of tckresample -endpoints; the -ends_only option in tckmap essentially duplicates that operation).

Would there be a way to extract the coordinates of the voxels with the highest streamline density?

There is, but there is more than one way to interpret “the voxels with the highest streamline density”.

If I take the simplest case, which is to find the voxel with the highest streamline density:

tckmap tracks.tck -template template.mif -ends_only endpoint_density.mif
mrthreshold endpoint_density.mif -top 1 max_endpoint_density_voxel.mif
maskdump max_endpoint_density_voxel.mif

You can modify the middle mrthreshold step to do something more advanced if you want; for instance, I could image thresholding at some value, doing a connected-component analysis, selecting the maximal value within each cluster, and then outputting those locations.

Cheers
Rob

1 Like

Hi Rob,

Thanks for your always clear and detailed explanations and interpretations of my rather confusing questions.

Indeed you are right, I was referring to the output of tckresample -endpoints. And yes, in my view I did perform unidirectional tracking, if this conforms to the use of the following call:

tckgen wmfod.mif output.tck \
-seed_image Left_STN_coreg.mif \
-exclude Right_STN_coreg.mif -exclude Right_NAc_coreg.mif -exclude Left_NAc_coreg.mif \
-act 5ttseg.mif -backtrack -crop_at_gmwmi \
-select 10000k -maxlength 250 -info

My aim is to find the streamlines seeded in the left subthalamic nucleus (in this case) that reach the cortex. Do you think then this approach is appropriate?
The tractography call is then repeated 3 other times, using as seed-image each of the exclusion masks listed in the above command (not sure if this is the most time efficient way of proceeding tho?).

tckmap tracks.tck -template template.mif -ends_only endpoint_density.mif
mrthreshold endpoint_density.mif -top 1 max_endpoint_density_voxel.mif
maskdump max_endpoint_density_voxel.mif

This looks like what I need, thanks!
I now understand the meaning of the reference image for tckconvert. Does the same apply to tckmap -template? Or what should the template be in this case?
I do have DWI and T1 images co-registered, but not resampled. Does this mean that, if I eventually want to have the max_endpoint_density_voxel.mif in the original T1 space, I have to apply a diff2struct transform?

Thanks again, I hope this is somewhat clear

Best,
Sam

I just wanted to follow up on this as I changed my approach to tractography, so the tckgen call I reported in my previous reply is no longer valid.

I thought a bit more about this and, because I am essentially identifying endpoints of streamlines based on their density, it seemed more appropriate to employ a method that would give me a reliable indication of this. Hence the decision to use SIFT2.
For this purpose, I am now doing the following:

  1. Whole-brain tractography with dynamic seeding, ACT, 10 million streamlines
  2. SIFT2
  3. Extract pathways of interest, with the following tckedit call:
tckedit ${subID}_wholebrain.tck ${subID}_L_STN.tck \
-tck_weights_in ${subID}_tck_weights.txt -tck_weights_out ${subID}_tck_weights_L_STN.txt \
-include ${subID}_L_STN_coreg.mif -include ${subID}_corticalmask_coreg.mif \
-exclude ${subID}_L_NAc_coreg.mif${subID}_R_NAc_coreg.mif \
-exclude ${subID}_R_STN_coreg.mif

tckedit is repeated for a total of 4 times to extract streamlines passing through left/right STN/NAc (respectively) and the cortex.
I am not sure this is the right place to ask this, but I am gonna throw in just a couple of doubts I have :innocent::

  • Is it always advisable to dilate your subcortical (and cortical ?) masks of 2 voxels to make sure streamlines are selected with tckedit?
  • Being interested in connections between yes, a small nucleus, but to the whole cortex, I decided not to combine whole-brain to targeted tractography before running SIFT2. However, wanting to extract connections from 4 nuclei, do you think it is a sensible choice? Would you then boost the whole-brain tractography to more than 10 million streamlines if I am just extracting the paths from that?

Now thinking back about your answer, I realise that I am not 100% sure I understand what you mean with unidirectional tracking, and how it would be implemented practically speaking in the command. Could you maybe specify this a bit more? Is there a specific option to pass to the tckedit command to do this? Do you generally want to do this?

This question still stands…

Cheers,
Sam

And yes, in my view I did perform unidirectional tracking, if this conforms to the use of the following call:

No, that’s bi-directional. Only way for tracking to be unidirectional is if you explicitly specify -seed_unidirectional, or use -seed_gmwmi and only ever seed from cortical and not sub-cortical GM. So if you run tckresample -endpoints, in most cases the vertices corresponding to the seed locations will no longer be present in the data, and it is not appropriate to refer to one streamline termination as the “start” and the other as the “end”.

My aim is to find the streamlines seeded in the left subthalamic nucleus (in this case) that reach the cortex. Do you think then this approach is appropriate?

Well you’re guaranteed that all of the streamlines will have been seeded in that structure based on the fact that you’ve only provided a single seed image corresponding to that structure. However it won’t be guaranteed that one of the two streamline endpoints will lie within that image, since tracking in both directions bidirectionally may have exited that image. It is at least guaranteed that it won’t be the case that both of those opposing projections from the seed point will exit the sub-cortical GM and enter the white matter, as per ACT prior #6 (the implementation is a little weirder given that the seed point is within sub-cortical GM, and therefore one of the projections has to be permitted to enter the WM, just not both). But following termination of the streamline having tracked, how far tracking is permitted to go in the opposite direction is dictated by the 5TT image, not the seed image.

More likely I suspect that for your experiment you want to be engaging unidirectional tracking.

The tractography call is then repeated 3 other times, using as seed-image each of the exclusion masks listed in the above command (not sure if this is the most time efficient way of proceeding tho?).

While you can provide multiple seed images in a single tckgen call, your current usage guarantees that every streamline intersects exactly one of these four regions, which would not be possible to guarantee if using a single tckgen call. So assuming that requirement is deliberate, no there isn’t a faster way to do it.

I now understand the meaning of the reference image for tckconvert. Does the same apply to tckmap -template ? Or what should the template be in this case?

The command-line options are named identically, but their operations are just as different as are the operations of these two commands. What the template image in the latter case “should” be depends on what you intend to do with those data. In any circumstance in which you want to perform a voxel-wise operation between two images, those images need to be defined on the same voxel grid; so if there is such an operation intended for the output of tckmap, that’s the voxel grid you should be using as the template.

I do have DWI and T1 images co-registered, but not resampled. Does this mean that, if I eventually want to have the max_endpoint_density_voxel.mif in the original T1 space, I have to apply a diff2struct transform?

Firstly, if you take a mask image containing one solitary voxel, and not only apply an affine transformation but also resample that data onto a different voxel grid, you no longer have an image that contains a solitary voxel. Depending on the interpolation method and the relative resolutions, you could have zero voxels, you could have more than one voxel, and you could have floating-point values rather than a binary mask.

Secondly, if the DWI and T1 images are coregistered but not resampled, then as long as you are using MRtrix3 tools, no, you don’t want to be applying any transformation. Such a transformation implies that you are taking information from one spatial location and moving it to a different spatial location; if the anatomical content of your two images already overlap (irrespective of the voxel grids on which they are defined), then shifting that image content in space would break that colocalisation. This is a consequence of MRtrix3 treating everything in “real” / “scanner” space; other softwares may have different demands around such things.

I thought a bit more about this and, because I am essentially identifying endpoints of streamlines based on their density, it seemed more appropriate to employ a method that would give me a reliable indication of this. Hence the decision to use SIFT2.

This is kind of like a switch from “most probable” trajectories based on probabilistic streamlines, to “most dense”, a distinction that often gets overlooked. The latter can be advantageous, but can also be unreliable when dealing with very small nuclei.

Is it always advisable to dilate your subcortical (and cortical ?) masks of 2 voxels to make sure streamlines are selected with tckedit?

It’s cortical masks that are a problem here rather than subcortical. But yes, until one of the better solutions to the problem of streamlines terminating at the GM-WM interface but ROIs being defined within GM voxels gets integrated, something like this is required. The alternative is to generate a parcellation image containing nodes corresponding to your nuclei of interest, one giant node corresponding to the cortex, and potentially another node containing other GM regions not of interest, and use tck2connectome followed by connectome2tck to extract streamlines corresponding to the connections of interest; this will utilise the radial search assignment mechanism that is the current default in tck2connectome, which is essentially akin to dilating ROIs but a fraction more accurate.

Being interested in connections between yes, a small nucleus, but to the whole cortex, I decided not to combine whole-brain to targeted tractography before running SIFT2. However, wanting to extract connections from 4 nuclei, do you think it is a sensible choice? Would you then boost the whole-brain tractography to more than 10 million streamlines if I am just extracting the paths from that?

I don’t think the distinction between 1 nuclei of interest vs. 4 is the deciding factor here. it’s moreso about the size of the smallest nuclei of interest. If a nuclei is so small and your streamline count so small that only a few streamlines intersect it, then no amount of SIFT2-ing will give you a reasonable estimation of the distribution of connectivity to the cortex.

Now thinking back about your answer, I realise that I am not 100% sure I understand what you mean with unidirectional tracking, and how it would be implemented practically speaking in the command. Could you maybe specify this a bit more? Is there a specific option to pass to the tckedit command to do this? Do you generally want to do this?

See tckgen's -seed_unidirectional option. In your original experimental description, its use would be advocated; but if you are instead performing whole-brain tractography, then you don’t want to use it.

Cheers
Rob

Hi Rob,

Thanks for the detailed answer. I followed up on your advice and realised that yes, unidirectional tracking was the most appropriate for my experiment.
I also in the end decided to combine whole-brain tractography with targeted tractography (total of 50M streamlines) → apply SIFT2 on the combined tractogram → apply tckedit to extract my tracts of interest.

I am eventually interested in finding the voxel with the highest streamline density, within different masks. For this purpose, I used tckmap -ends_only endpoint_density.mif, feeding the SIFT2 weights extracted for each of the extracted tracts.

Probably very basic question but, does the endpoint_density.mif image embed information about the weights then?
If I want to perform some operations on this image (e.g. subtract the endpoint_density.mif of two different tracts to identify unique connections, and then apply a threshold to find the voxel with the highest streamline density within certain masks), do I have to use the SIFT2_weights.txt in any way?
If yes, how?

Thanks in advance!
Cheers,
Sam

Hi Sam,

I am eventually interested in finding the voxel with the highest streamline density, within different masks.

May I suggest mrthreshold -top 1 | maskdump.

Probably very basic question but, does the endpoint_density.mif image embed information about the weights then?

No; even the fact that streamlines were used in the generation of that image is more-or-less lost at that point. Plus, even if you had all the streamline weights within that file, you don’t have any information about the streamlines to which they were attributed, so it’s not clear how they would actually be of any use.

If I want to perform some operations on this image (e.g. subtract the endpoint_density.mif of two different tracts to identify unique connections, and then apply a threshold to find the voxel with the highest streamline density within certain masks), do I have to use the SIFT2_weights.txt in any way?

This all consists entirely of operations that are based on voxel-wise intensities. So as previously, it’s wholly agnostic to the fact that streamlines were used at all in the calculation of those input voxel-wise intensities.

Cheers
Rob

Hi Rob,

Thanks for your quick reply.

Yes, this is exactly what I was doing, thanks!

Mmm ok :pensive: I performed SIFT2 to be able to make “quantitative inferences” with a bit more confidence about the “strength” of the connections between a subcortical seed and the cortex. My end goal is to eventually apply neuromodulation techniques on the cortical target that I deem to be more “strongly / probably / densely” (whatever term is the most appropriate in this case) connected to that subcortical seed.
I am just wondering how this top -1 voxel should be interpreted then, if the information conveyed by the weights computed with SIFT2 is more or less lost.

I guess that I am back to what is the voxel with the most probable trajectories to the seed nucleus?

If I wanted to make quantitative assessment of the connections between my subcortical nuclei of interest and the cortex (taking advantage of the use of SIFT2), should I use a tck2connectome kind of approach?

Cheers,
Sam

Mmm ok :pensive: I performed SIFT2 to be able to make “quantitative inferences” with a bit more confidence about the “strength” of the connections between a subcortical seed and the cortex.

I think this may be the result of me misunderstanding your original question.
I had interpreted “embed information” as retaining information about the weights of individual streamlines; which would be pointless, as you would have per-streamline weights but no streamlines information within that image.
But if by this you instead meant something like the endpoint density being “influenced by” the calculated per-streamline weights, then this is the case, as long as you use the -tck_weights_in option. Instead of the intensity of each voxel being the number of streamline terminations in that voxel, it will instead be the sum of the weights of those streamlines terminating in that voxel.

I guess that I am back to what is the voxel with the most probable trajectories to the seed nucleus?

The “most probable” connection from any specific location is based on generating a large number of probabilistic streamlines seeded from that location, and quantifying the fraction of those streamlines that intersect any particular target; i.e. the PICO method. This is a common usage of FSL’s probtrackx I think. I am slightly critical of this interpretation when the source location is a nucleus rather than an infinitessimal point, as there’s a spatial averaging operation hidden in there, but it’s still a legitimate approach.

If I wanted to make quantitative assessment of the connections between my subcortical nuclei of interest and the cortex (taking advantage of the use of SIFT2), should I use a tck2connectome kind of approach?

I think in general that’s what I’d advise; assuming that you are content with parcel targets rather than voxel targets. If you need the higher spatial resolution, a surface-based cortical target would be preferable; in the absence of a software implementation of such, you can use voxel-wise streamline endpoint density, but there will be some bias as the effective surface area available to cause streamline terminations will vary across different voxels.

Rob