ENH: Remove dependency on nibabel.trackvis

[jhlegarreta]

Remove dependency on nibabel.trackvis: the module was deprecated in NiBabel 2.5.0 and removed in version 4.0.0. So this patch set adapts the tractography.trackvis code so that it uses the modern API.

Add a tolerance when comparing streamline data in equal_tracts to account for some precision loss since data is transformed to float32 (vs float64) when applying the operations of the NiBabel API.

Fixes:

tract_querier/tractography/tests/test_tractography.py::test_saveload_trk
  site-packages/nibabel/deprecated.py:35:
 DeprecationWarning:
 The trackvis interface has been deprecated and will be removed in v4.0; please use the 'nibabel.streamlines' interface.
    mod = __import__(self._module_name, fromlist=[''])

and

tract_querier/tractography/trackvis.py:7: in <module>
    from nibabel import trackvis
E   ImportError: cannot import name 'trackvis' from 'nibabel'
 (/opt/hostedtoolcache/Python/3.9.20/x64/lib/python3.9/site-packages/nibabel/__init__.py)

and related errors raised for example in:
https://github.com/demianw/tract_querier/actions/runs/12195905001/job/34022456528#step:6:85 and
https://github.com/demianw/tract_querier/actions/runs/12195905001/job/34022456201#step:6:32

[jhlegarreta]

Commit 8adff56 makes the Python 3.9, 3.10, and 3.11 minimal versions pass:

That commit checks as well the old NiBabel trackvis module-based results and the new NiBabel API results.

To avoid missing the code, I've also put it into a gist:
https://gist.github.com/jhlegarreta/50c1dd689f2482d5837fb74b877ab01d

I will push a commit removing the checks and the dependency on the old NiBabel trackvis module that should make CIs happier.

[jhlegarreta]

Commit c977552 makes CIs happier:

Python 3.12 minimal failure is related to PR #60.

I am not sure whether the image_orientation_patient property,
https://github.com/demianw/tract_querier/pull/64/files#diff-4924744de8b2a10232cf247057a47123480592cf18425a8c237b473d0e89fa60R104

which I had to borrow from NiBabel to match exactly the trackvis-based value since I could not assign it as in
https://gist.github.com/jhlegarreta/50c1dd689f2482d5837fb74b877ab01d#file-tract_querier_tractography_trackvis-py-L91

is currently used in NiBabel's trk implementation: a search in the repository shows two matches:

with no assignments.

I think tract_querier should eventually move to using DIPY's StatefulTractogram which was developed precisely because people at the time had the same issues this patch set brings to surface. Although issues of many sorts still persist with the data and software out there, the implementation is more reliable. Also I was not involved in that effort or the NiBabel transition, so my knowledge on this is limited, so the patch set should be considered with caution.

Following the trackvis vs new NiBabel API checks in the previous gist, I added some checks using DIPY's StatefulTractogram so that it can be available for when revisiting this:
https://gist.github.com/jhlegarreta/36a8ac47b5d179d889a3347ec4d75b1c

That said, DIPY would be an additional, somehow heavy dependency for the purposes of tract_querier, but this where things are at now.

@demianw best would be to try the 93ab73f version against the one in here (c977552) with some known, in vivo data with known results. I do not have any known data/known results, so leaving this up to you @demianw if you can afford the time. Otherwise, this is good to go on my end.

[jhlegarreta]

Converting back to draft; will check one more thing.

[jhlegarreta]

So I converted this back to draft because I wanted to check the solution with a non-identity affine, and as I was suspecting the proposed solution does not work, unfortunately.

The affine I used:

[[ 6.00468478e-01 -4.53037826e-01  5.47165447e-01  2.00000000e+01],
 [ 1.04004191e+00  6.86156843e-01 -3.28425633e-02  5.00000000e+01],
 [-6.93361274e-01  6.36892999e-01  4.24595332e-01 -3.00000000e+01], 
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  1.00000000e+00]]

was one that I created with the following arbitrary values

rotation_angles = (45, 30, 60)  # in degrees
scaling = [2.0, 1.5, 1.0]
translation = [20, 50, -30]

It turns out that debugging tractography_to_trackvis_file I see that the original streamlines and the ones saved with the updated NiBabel version are shifted by some factor:

(I changed the test_tractography.setup_module rng.standard_normal calls to rng.random so that the streamline coordinates are contained in a [0,1) cube and the effect could be visualized more easily).

IMO this has to do with the trackvis data assuming that the origin of the coordinates is at the center of the voxel:
https://github.com/nipy/nibabel/blob/3.0.0/nibabel/streamlines/trk.py#L238-L240

When trying to deal with that and considering what the StatefulTractogram does:
https://github.com/dipy/dipy/blob/db4e8990bbe365c7507cb4d012ba7aac9a1b93b0/dipy/io/stateful_tractogram.py#L749-L770
https://github.com/dipy/dipy/blob/db4e8990bbe365c7507cb4d012ba7aac9a1b93b0/dipy/io/stateful_tractogram.py#L545
https://github.com/dipy/dipy/blob/db4e8990bbe365c7507cb4d012ba7aac9a1b93b0/dipy/io/stateful_tractogram.py#L551

inserting the code below in to the gist (see below for the link):

(...)
affine_voxmm_to_rasmm = nib.streamlines.trk.get_affine_trackvis_to_rasmm(trk_header_new_nbb)

from nibabel.affines import apply_affine
shift = numpy.asarray([0.5, 0.5, 0.0])
tmp_affine = numpy.eye(4)
tmp_affine[0:3, 0:3] = affine_voxmm_to_rasmm[0:3, 0:3]
shift = apply_affine(tmp_affine, shift)
shift *= -1
streamlines = nib.streamlines.ArraySequence(tractography.tracts())
streamlines = nib.streamlines.ArraySequence(map(lambda x: x - shift, streamlines))

# Save the data using the new NiBabel interface
#streamlines = tractography.tracts()
data_per_streamline = None
data_per_point = data
tractogram_nb = nib.streamlines.Tractogram(
    streamlines=streamlines,
    data_per_streamline=data_per_streamline,
    data_per_point=data_per_point,
    affine_to_rasmm=affine_voxmm_to_rasmm,
)

(...)

The two sets of streamlines (original tract_querier code vs the one in the commit that was previously proposed as a solution in this PR c977552) are translated exactly half the size of the bounding box:

I have left the entire modified tract_querier.tracvkis in this gist:
https://gist.github.com/jhlegarreta/a60d96ae1836b31d69e7ae22e021e7e7

I have spent some time trying to fix this without success. So @demianw if you happen to have some insight on what is happening, it would be much appreciated.

[effigies]

@MarcCote I hate to bother you, but if anybody might instantly see the issue here, it would be you.

[MarcCote]

Hey @effigies no problem. I vaguely recall that the old nibabel's trackvis implementation might have not considered that Trackvis (the software, hence the file format) assume (0,0,0) is at the center of the voxel.

If I understand correctly the "issue" here is why the transformation (similar to what is done in the StatefulTractogram) does not bring the two tractograms on top of each other? What software are you using to do the visualization? Some software will take into account the metadata contains in the header of the .trk.

Lighting up the tractogram's beacon @frheault

[jhlegarreta]

Thanks for chiming in @effigies and @MarcCote.

If I understand correctly the "issue" here is why the transformation (similar to what is done in the StatefulTractogram) does not bring the two tractograms on top of each other?

Somehow; I do not fully see where the underlying issue lies. Using the non-identity, arbitrary affine I posted, the StatefulTractogram refuses to read the pre-4.0.0 TRK file as well. So I do not know what is missing in the pre-4.0.0 version that is preventing the file from being read, and what step I may be missing in the pre- to post- transition if the pre- is a perfectly valid TRK file.

What software are you using to do the visualization?

Mi-Brain.

[frheault]

Hello,

Sadly, I am not familiar (and/or I forgot) with how the TRK was written in the older version of Babel. The load_tractogram() function assumes a lot of things.

It is pretty common when I deal with older file that I need to load the streamline with nib.streamlines.load() and then create a StatefulTractogram sft = StatefulTractogram(tractogram.streamlines, 'my_nifti.nii.gz', space.RASMM, origin=origin=Origin.NIFTI) and then try to switch to origin=Origin.TRACKVIS (or vice-versa)

This is common because some software saves the TRK assuming the wrong origin but uses software with the same assumption. This is not a problem with the affine (the affine should be identical to the NIFTI that match with the TRK).

There is an even worse possibility, sometimes a VOX space shift is made in RASMM (adding 0.5 but when the streamlines are in their final RASMM space) or vice versa (A RASMM shift is made while the streamlines are in VOX). This is more complex and requires a lot of tinkering to fix the TRK so they have no negative coordinates in VOX space and align on top of their matching NIFTI.

[jhlegarreta]

Resuming this. Sorry for the long delay.

I wish some migration guide had been written at the time, and all the tinkering involved in these cases could be put into a documented resource that the community could freely and openly access.

Based on the conversation, I will then assume that whatever nibabel post-3.0.0 is what we should expect, and that reproducing the output as written by nibabel pre-3.0.0 is not the right thing to do.

As for the non-identity test, I am contenting myself using an affine that I borrowed from a clinical DWI that I have. After all, all tractography I/O tests in DIPY (at least, those that I considered while working on this patch set) use identity affines.

So I will amend the commit before the end of the week after some more testing and mark the PR as ready.

[MarcCote]

I wish some migration guide had been written at the time, and all the tinkering involved in these cases could be put into a documented resource that the community could freely and openly access.

I wished I would have thought of that.

In the original PR nipy/nibabel#391, I aimed to be as consistent as possible with the trackvis.py implementation with proper ground truth .trk files (e.g. https://github.com/nipy/nibabel/blob/master/nibabel/tests/data/simple.trk). Then, given how little information was provided on the specs for the .trk format (https://trackvis.org/docs/?subsect=fileformat), I reverse engineered how original TrackVis (the software) handles the different .trk properties (e.g. affine, dimension, etc). That's probably when ninabel.streamlines diverged a bit from trackvis.py. Also, as François pointed out, given the ambiguous specs for .trk different softwares implemented the format differently.

[jhlegarreta]

I have force pushed a new commit with what I believe is the right patch set for this. I have also removed the NiBabel upper limit restriction.

As said in #64 (comment), I have contented myself with the identity affine matrix for the synthetic streamlines. When visualizing them before this patch set and after, there is the voxel center/corner translation that was expected asusming M-A's recall in #64 (comment) was correct.

Now, I did try with a non-identity affine matrix for in-vivo data: took participant sub-1012 from the TractoInferno dataset available on OpenNeuro, computed the FreeSurfer wmparc and applied the tract_querier example UF extraction query. For the records, the affine matrix of that data is:

[[ -1.      0.      0.     71.674]
 [  0.      1.      0.    -66.918]
 [  0.      0.      1.    -56.29 ]
 [  0.      0.      0.      1.   ]]

So one rotation and translations.

The results with the code before this patch set and after are shown below:

Pre-/post- (yellow/cyan) patch set overlay:

Pre-/post- patch set overlay and the original (pink) left UF available from OpenNeuro:

There is a perfect overlap between the results before this patch set and after (using NiBabel 3.0.0 on both cases), they are located at the expected anatomical location, and overlap with the original tractography data. There is some streamlines that are not extracted by the new patch set due to the mentioned voxel center/corner difference.

Thus, my conclusion is that the code before this patch set and after provides the same results, and this PR is ready to be reviewed and merged.

[jhlegarreta]

Re #64 (comment) I understand. Thanks @MarcCote.

test_tract_querier_query is now raising a warning about points falling outside the image:
https://github.com/demianw/tract_querier/actions/runs/13992325927/job/39178873335?pr=64#step:6:65

tract_querier/tests/test_tract_querier.py::test_tract_querier_query[inprocess]
65 /home/runner/work/tract_querier/tract_querier/tract_querier/tract_label_indices.py:250: UserWarning: Warning tract points fall outside the image
66 warnings.warn("Warning tract points fall outside the image")

Trying to load the tractogram at issue (IIT3mean_left_hemisphere_small.trk) with DIPY:

from dipy.io.streamline import load_tractogram
sft = load_tractogram(args.in_tractogram, 'same', bbox_valid_check=True)

results in:

ValueError: Bounding box is not valid in voxel space, cannot load a valid file if some coordinates are invalid.
Please set bbox_valid_check to False and then use the function remove_invalid_streamlines to discard invalid streamlines.

Have tried removing (cutting) the invalid streamlines (remove_invalid_streamlines), and this results in all streamlines being removed. Inspecting the data attributes:

----------------------------------------------------------
Looking at /tmp/tract_querier/IIT3mean_desikan_2009.nii.gz
----------------------------------------------------------
Dimensions:      (182, 218, 182)
Voxel order:     LAS
Affine matrix:
[[  -1.   -0.    0.   90.]
 [  -0.    1.   -0. -126.]
 [   0.    0.    1.  -72.]
 [   0.    0.    0.    1.]]
Voxel size:      (1.0, 1.0, 1.0)
----------------------------------------------------------------
Looking at /tmp/tract_querier/IIT3mean_left_hemisphere_small.trk
----------------------------------------------------------------
Number of streamlines:    6783
Affine matrix:
[[1. 0. 0. 0.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]
 [0. 0. 0. 1.]]
Dimensions:               (1, 1, 1)
Origin:                   Origin.NIFTI
Space:                    Space.VOXMM
Voxel order:              RAS
Voxel size:               (1.0, 1.0, 1.0)

So I will not try to tinker with this data. I am happy to upload to OSF the Tractonferno paticipant tractography data and parcellation files and use those for the test.

So I'd change the test data to avoid this issue @demianw.

But the PR is ready on my end.

[demianw]

Hello, I will check this tomorrow. I want to make sure that the affine matrix and other coordinates are well kept

[jhlegarreta]

Were the checks successful @demianw ? Thanks.

[demianw]

It's good thanks ! Sorry for the delay

[demianw]

The Tractonferno idea for the test is actually agreat one @jhlegarreta

[jhlegarreta]

Thanks @demianw. I will have a look at what we can do for the testing data as time permits.