affine matrix issue while saving trackvis files

[daniel-ge]

I am attempting to write files in Trackvis-Format (trk). So loading a trk-file and saving it without any further changes should work. But unfortunately it fails in the following example:

Failing example

import nibabel.streamlines.trk as nitrk 
from numpy.testing import assert_equal, assert_almost_equal
import numpy as np
from numpy.linalg import inv
import urllib

trkname = "/tmp/myfile.trk"
trkname_clone = "/tmp/myfile_clone.trk"

# please let me know if the link of the example file is expired
url = "https://gigamove.rz.rwth-aachen.de/d/id/4hd7gVeFc4GweS/dd/100"
urllib.URLopener().retrieve(url, trkname)

trk = nitrk.TrkFile.load(trkname) # just loading...
trk.save(trkname_clone) # ... and saving
trk_clone = nitrk.TrkFile.load(trkname_clone) # let's load the content for comparison

# the contents of the files should be equal:
assert_equal(trk.header,trk_clone.header)
for i in range(len(trk.streamlines)):
    assert_almost_equal(trk.streamlines[i],
                        trk_clone.streamlines[i], decimal=4)
"""
AssertionError: 
Arrays are not almost equal to 4 decimals

(mismatch 100.0%)
 x: array([[-111.3645,  160.6524,   40.247 ],
       [-112.3261,  161.0878,   40.6915],
       [-112.4026,  161.1225,   40.7269],...
 y: array([[-4416.592 ,  4464.1333,  4995.547 ],
       [-4387.707 ,  4434.735 ,  4961.6387],
       [-4385.4077,  4432.3945,  4958.94  ],...
"""

Working example

I figured out that this error is related to the creation of the affine transformation in get_affine_rasmm_to_trackvis. After redefinition of this function the saved file contains the same data as the original one.

def f(hdr):
    T = nitrk.get_affine_trackvis_to_rasmm(hdr)
    T2 = np.eye(4)
    T2[:3,:3] = inv(T[:3,:3])
    T3 = np.eye(4)
    T3[:3,3:4] = T[:3,3:4]*-1
    return np.dot(T2,T3)

nitrk.get_affine_rasmm_to_trackvis = f

trk = nitrk.TrkFile.load(trkname)
trk.tractogram.affine_to_rasmm = np.eye(4) # because it was only ROUGHLY eye(4)
trk.save(trkname_clone)
trk_clone = nitrk.TrkFile.load(trkname_clone)

# the contents of the files should be equal:
assert_equal(trk.header,trk_clone.header)
for i in range(len(trk.streamlines)):
    assert_almost_equal(trk.streamlines[i],
                        trk_clone.streamlines[i], decimal=4)

"""
No Error :)
"""

[effigies]

Hi @daniel-ge, thanks for the report. I don't have any experience with these file formats, so I'll leave evaluation to others. Would you have any interest in submitting a pull request? We almost certainly have round-trip tests, so adding your failing file to the battery would be a useful addition, along with fixes to ensure that your file is correctly loaded (if valid) or rejected (if actually invalid).

@MarcCote Do you have time to check on this and potentially review a PR?

[MarcCote]

Hi, I can have a closer look next week. I'd be happy to review a PR fixing that issue.

[MarcCote]

@daniel-ge the link to the data is now invalid.

[effigies]

Hi @daniel-ge, can you update this issue with a fresh link?

[effigies]

@daniel-ge Just a bump on this. Can the data be made available again? I can grab it so we're less sensitive to delays next time.

[sitek]

I may be having a related issue where nibabel.streamlines.save appears to be overwriting the given affine with np.eye(4). I described the issue on neurostars, with help from @arokem : https://neurostars.org/t/dipy-streamlines-not-aligning-with-image-space/3889/5

The key part is this:

I set the DWI-derived affine as the new affine:

tractogram = Tractogram(sl, affine_to_rasmm=affine)
tractogram.affine_to_rasmm                                                                 
Out[25]: 
array([[  -1.04999995,    0.        ,    0.        ,   90.        ],
       [   0.        ,    1.04999995,    0.        , -126.        ],
       [   0.        ,    0.        ,    1.04999995,  -72.        ],
       [   0.        ,    0.        ,    0.        ,    1.        ]])

and I saved it:

save(tractogram, 'tractogram_affine-fa.trk')

but when I load this new .trk file, the affine is back to an eye:

In [17]: sl_aff_fa, hdr_aff_fa = load_trk('tractogram_affine-fa.trk')                     
In [18]: hdr_aff_fa['voxel_to_rasmm']                                                               
Out[18]: 
array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.]], dtype=float32)

EDIT: my issue is not a nibabel bug – I resolved it by creating a header and feeding it to save() along with the tractogram object (see this dipy issue)

[MarcCote]

Thanks for letting us know. There is some work in progress about making tractogram manipulation more robust in Dipy (dipy/dipy#1812). Hopefully, we'll be able to bring most of the checks in Nibabel.