Bug in tcknormalise (when using warpcorrect)

Dear all,

I think I have found a bug in tcknormalise when use with warpcorrect. Below is a description of the bug.

Warp correct:

Replaces voxels in a deformation field that point to 0,0,0 with
nan,nan,nan

It seems that the vertices that now point to nan nan nan are interpreted as empty streamlines (see below).

Command:
tckstats tckwarp_nocorr.tck -nthreads 0 (to keep the same ordering of tract input and output)
Output:

Command:
tckstats tckwarp_corr.tck -nthreads 0
Output:

As you can see, the one that has used the corrected warp has 123 “extra” tracts in the data section of the file.

I probed the file in Matlab:

tckwarp_nocorr.tck data:

tckwarp_corr.tck data:

As you can see the two are equivelent up to cell 1562. Cell 1563 seems to have the tracts that point to [0 0 0] or [nan nan nan]. streamlines that used the corrected warp seem to have added cells (1565 to 1567) that are interpreted as empty.

I have counted the non-empty cells in the corrected tracts file and they do indeed sum to 5000 (the number of streamlines I started) as expected.

I don’t really think that this will cause any problems to calculations but it will produce the annoying warning that I showed above. Do you know about this bug and is there a fix?

Claude

OK, so what’s happening here is that any streamlines vertices outside the deformation field (i.e. voxels that now contain NaNs), get replaced by NaNs rather than zeros. However, a triplet of NaNs is used as the streamline delimiter in the tck file format. So the streamline that strayed outside the deformation field get truncated, and subsequent vertices are interpreted as empty (zero-length) streamlines. I agree this is not ideal behaviour, we should fix that. My only concern is figuring out what to do about streamlines that exit the deformation field, but subsequently re-enter it - I think the only sensible solution there is to split them and keep both valid parts, but I’m not sure… It looks like that wouldn’t be a problem in your case, given that the count is correct for the non-empty streamlines, but we do have to worry about these edge cases…

Would it perhaps be sensible to map the vertices that exit the deformation field to its nearest neighbor inside the field? Would be trickier to code than just replacing zeros with NaNs but could be a solution if you want to ensure that streamlines don’t split… unless you decide that it is better to split the streamlines, in which case my suggestion is a bad idea :slight_smile:

Yes, I suspect there isn’t a ‘right’ answer to this problem - ultimately, there is no information about where these bits of the streamlines should go.

Personally, I reckon anyone affected by this would be better off fixing the problem at source: by making sure the deformation field does cover the entire domain where streamlines can be found. In my experience, this issue arises when the initial no-warp image (generated in template space) is too small, so that once warped to the individual subject’s native space, some edge regions are left undefined (i.e. the deformation would place them outside the region covered in the template). To avoid this, it’s best to use a larger image as the template argument to warpinit, to ensure it covers the full extent of the brain and spine in all of the subjects. A T1 anatomical image is typically good for this if you have one in template space. Otherwise you can use mrpad to generate an larger image than what you’d used previously. Hopefully that’ll avoid the issue in the first place…