DIPY streamlines not aligning with image space

I’m using DIPY’s workflows to estimate DTI and CSD and create streamlines from these inputs. Each node is creating the appropriate inputs, but the streamlines don’t appear in the diffusion image space, despite taking DTI and CSD outputs as inputs: (full code below)

wf.connect(recon_dmri, 'out_pam', tracker, 'pam_files')
wf.connect(recon_dti, 'out_fa', tracker, 'stopping_files')
wf.connect(masker, 'mask_filename', tracker, 'seeding_files')

After loading the streamlines in TrackVis and importing the MD image created in the DTI workflow node, the MD image fails to appear: (this also happens loading both tractogram.trk and md.nii.gz in dipy_horizon)

Viewing the MD and other DTI derivatives on their own is fine, leading me to question the space mismatch during the tractogram generation process.

And actually when I first load the tractogram in TrackVis, it’s hyper-zoomed in, doubling down my space thoughts.

Any ideas how to fix the tractogram space issue? @Ariel_Rokem

Here’s the full code:

1 Like

Hmm. I am not 100% sure. What is the value of the affine stored in the trk file header?

from dipy.io.streamline import load_trk
sl, hdr = load_trk('/path/to/fname.trk')
array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.]], dtype=float32)

I wonder whether this is the problem, and it’s expecting to have the affine inherited from the DWI data nifti image. What happens if you change it and re-save?

Hi @Ariel_Rokem, picking this up again – still haven’t figured out where I’m going wrong.

I originally thought that the issue was with this line in the dipy workflow, which saves the affine as np.eye(4).

I loaded in the stopping image (fa.nii.gz) – here’s its affine:

array([[  -1.04999995,    0.        ,    0.        ,   90.        ],
       [   0.        ,    1.04999995,    0.        , -126.        ],
       [   0.        ,    0.        ,    1.04999995,  -72.        ],
       [   0.        ,    0.        ,    0.        ,    1.        ]])

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

tractogram = Tractogram(sl, affine_to_rasmm=affine)
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 = sl, hdr = load_trk('tractogram_affine-fa.trk')                     
In [18]: hdr_aff_fa['voxel_to_rasmm']                                                               
array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.]], dtype=float32)

Am I doing something wrong here, or is nibabel.streamlines.save doing something unexpected?

EDIT: looking into it further, this may be related to an existing nibabel issue – I’m posting there as well to get their feedback.

@Ariel_Rokem, I think I resolved the issue (in my case at least). LocalFiberTrackingPAMFlow() creates a Tractogram object and then saves this with nibabel.streamlines.save(tractogram, filename).

The save() function takes either a Tractogram object or a TractogramFile object – if it’s given a Tractogram without also giving it a header, it will construct a header (including VOXEL_TO_RASMM) with default values.

The solution, then, is either to first turn the Tractogram object into a TractogramFile object with a correct header and then save the TractogramFile object, or to pass a header to save() along with the Tractogram object.

The header is created like this:

from nibabel.orientations import aff2axcodes
import nibabel as nib

img = nib.load(stopping_path)

affine = img.affine
vox_size = img.header.get_zooms()[:3]
fov_shape = img.get_fdata.shape[:3]

 hdr = {}
hdr[Field.VOXEL_TO_RASMM] = affine.copy()
hdr[Field.VOXEL_SIZES] = vox_size
hdr[Field.DIMENSIONS] = fov_shape
hdr[Field.VOXEL_ORDER] = "".join(aff2axcodes(affine))

and then either option 1:

tractogram = Tractogram(streamlines_generator, affine_to_rasmm=np.eye(4))
trk_file = nib.streamlines.TrkFile(tractogram, header=hdr)
save(trk_file, 'tractogram_correct_affine.trk')

or option 2:

tractogram = Tractogram(streamlines_generator, affine_to_rasmm=np.eye(4))
save(trk_file, 'tractogram_correct_affine.trk', header=hdr)

**note that the Tractogram object is created with a blank affine (np.eye(4)) but the header takes the image affine

Now the streamlines are in the same location and have the same size as the DWI:

I created a github issue here

1 Like