Are the bvecs incorrect for PAR / REC files

[matthew-brett]

See discussion at the end of #275.

We are looking at the DTI data in https://bitbucket.org/nipy/nitest-parrec and in particular at the PAR / REC converted bvectors for the file DTI_15DIR_SENSE_4_1.PAR in that repo.

With the code at the current time of writing, parrec2nii --bvs DTI_15DIR_SENSE_4_1.PAR gives bvectors of:

0.0 -1.0 0.0 0.0 0.179 0.064 -0.71 -0.619 -0.242 0.259 0.817 0.844 0.263 0.0 -0.745 -0.973 0.0 
0.0 0.0 -1.0 0.0 0.111 -0.377 -0.052 0.438 -0.784 0.618 -0.17 -0.526 -0.955 -0.969 -0.666 -0.232 0.0 
0.0 0.0 0.0 1.0 -0.978 -0.924 -0.701 -0.651 -0.571 -0.742 -0.551 -0.106 -0.139 0.248 0.024 0.021 0.0

Running dcm2nii on the DICOM folder of this repository gives a converted file x20141118_192645WIPDTI15DIRSENSEs401a1004.nii.gz which appears to be equivalent. This file is flipped in the second (Y) direction compared to DTI_15DIR_SENSE_4_1.PAR. The converted bvectors are:

0 1 0 0 -0.17889153957366 -0.06349743902683 0.71040332317352 0.61909401416778 0.24240906536579 -0.25890484452247 -0.8168773651123 -0.84379339218139 -0.26261380314827 0.00009999634494 0.74529498815536 0.97256475687027 
0 0 -1 0 0.11129473149776 -0.37668478488922 -0.05162931606173 0.43849575519561 -0.78432935476303 0.61801159381866 -0.16969530284404 -0.52609586715698 -0.95485013723373 -0.9688646197319 -0.66629552841186 -0.23169161379337 
0 0 0 -1 0.97755372524261 0.92416268587112 0.70189851522445 0.65149372816085 0.57102131843566 0.74231392145156 0.55128473043441 0.10599917173385 0.13890728354454 -0.24759095907211 -0.02419983781874 -0.02089924365282 

These are the same as the PAR / REC bvecs multiplied by -1, but with the Y sign inverted, as you'd expect from the different voxel orientations of the two images.

When I run FSLs dtifit on the dcm2nii data, and view the first eigenvector lines, these look correct, and aligned with the corpus callosum.

If I do the same on the converted PAR / REC files, the Y direction of the tensor looks flipped.

If I flip the Y axis of the PAR / REC converted image, and flip the sign of the Y bvec value, save this out and then run dtifit, run fslview the Y direction of the tensor looks correct.

Is it possible that FSL insists on canonical P-A image voxel orientation?

@mrjeffs - how did you check the tensor orientation?

[matthew-brett]

Checking with dipy suggested that the directions are correct - so it looks as though FSL does not tolerate A-P flipped images. We could either ship as is, or reorient the images / bvecs to be compatible with FSL.

[larsoner]

+1 for outputting data that is properly interpreted by FSL, even though it's annoying :(

[matthew-brett]

Script to view tensors. Run on fixed_dti.nii, a version of DTI_15DIR_SENSE_4_1.PAR from nitest-parrec with the last spurious volume removed:

import numpy as np
import nibabel as nib

from dipy.core.gradients import gradient_table
from dipy.segment.mask import median_otsu
import dipy.reconst.dti as dti
from dipy.reconst.dti import fractional_anisotropy, color_fa
from dipy.data import get_sphere
from dipy.viz import fvtk

gtab = gradient_table('fixed_dti.bvals', 'fixed_dti.bvecs')
img = nib.load('fixed_dti.nii')

data = img.get_data()
maskdata, mask = median_otsu(data)
tenmodel = dti.TensorModel(gtab)
tenfit = tenmodel.fit(maskdata)
FA = fractional_anisotropy(tenfit.evals)
FA[np.isnan(FA)] = 0
FA = np.clip(FA, 0, 1)
RGB = color_fa(FA, tenfit.evecs)
sphere = get_sphere('symmetric724')
ren = fvtk.ren()
evals = tenfit.evals[42:86, 61:88, 26:27]
evecs = tenfit.evecs[42:86, 61:88, 26:27]
cfa = RGB[42:86, 61:88, 26:27]
cfa /= cfa.max()
fvtk.add(ren, fvtk.tensor(evals, evecs, cfa, sphere))
fvtk.show(ren)

[matthew-brett]

More info - when using the current master branch parrec2nii conversion, and running dtifit followed by fslview, the tensors look flipped in Y - the original observation. Reading the V1 tensor orientation for a voxel in the corpus callosum that should be pointing from right lower to left upper on the slice gives:

In [5]: v1 = nib.load('fsl_base_V1.nii.gz')
In [6]: v1.get_data()[71, 79, 30]
Out[6]: array([-0.65730137, -0.569175  , -0.49395818], dtype=float32)

I don't know what axes these refer to, but if these are using same axes as the bvecs, that is, the voxel directions of the image, then I think this is incorrect, because it points to the left (minus X) lower (minus Y) direction on the image.