Visualization problem with tensors ?

[skoudoro]

From @carpre on January 4, 2018 22:11

Hi,

I did an experiment with 2 gradient tables. The second gradient table is the same than the first but with x and y flipped (-x and -y).

I computed the tensors and tried to plot 3 slices (coronal, axial and sagittal) using vtk as shown in the examples.

I can see that I get different eigen vectors for the 2 gradient tables when I print them.
Original evec [[[-0.71158536 0.21049986 0.67032535]
[ 0.69739829 0.09573807 0.71026041]
[ 0.08533406 0.97289467 -0.21492805]]]
Flipped evec [[[-0.71158536 -0.21049986 0.67032535]
[ 0.69739829 -0.09573807 0.71026041]
[-0.08533406 0.97289467 0.21492805]]]
However, I get the exact same tensors for the 3 planes when plotting the same slices.

I am not sure to understand how I can get the same tensors oriented the same way when my evecs are not in the same direction ? If I would have a flip in the 3 axes, it would obviously be the same direction, but not if only 2 out of 3 are flipped ? I did the same processing with MRtrix and there is a clear flip in sagittal when I plot the tensors with mrview.

I assume it's the way I am plotting my tensors, but I can't figure out what's wrong.
This is how I am doing it:
tenmodel = TensorModel(gtab)
tenfit = tenmodel.fit(data)
RGB = color_fa(tenfit.fa, tenfit.evecs)
sphere = get_sphere('symmetric724')
ren = fvtk.ren()
fvtk.add(ren, fvtk.tensor(tenfit.evals[:,:28], tenfit.evecs[:,:28], RGB[:,:28]/RGB[:,:,28].max(), sphere))
fvtk.record(ren, n_frames=1, out_path='out.png', size=(1000, 1000))

For the other planes, I am making sure to give it as a NxMx1x3 or NXMx1x3x3 like so:
evals = tenfit.evals[:, 40:41, :].squeeze()
evecs = tenfit.evecs[:, 40:41, :].squeeze()
rgb = (RGB[:, 40:41, :]/RGB[:, 40:41, :].max()).squeeze()
ren = fvtk.ren()
fvtk.add(ren, fvtk.tensor(evals[:, :, np.newaxis, :], evecs[:, :, np.newaxis, :, :], rgb[:, :, np.newaxis, :], sphere))

Doing that, the tensors look exactly the same in the 3 planes, even though the gradients are suppose to be flipped in x and y and that the evecs are not the same (some signs are different). Any idea what I am doing wrong ?

Also, is there a way to put the corresponding FA slice under the tensors so it's easier to look at ? I tried adding an actor using vtk.add(ren, fvtk.slicer(fa[:, :, 28])) before the tensors, but the FA is not at the right place under the tensors.

Thank you very much for your help :)

Copied from original issue: dipy/dipy#1396

[samcoveney]

I don't think this issue was ever sorted out.

If you plot tensors using Fury, and then plot the first eigenvectors using arrows also in Fury, you can see that there is an issue with the direction. The tensors do not point in the same direction as the eigenvectors.

This seems like quite a serious bug?

[samcoveney]

Error can be reproduced as follows

import numpy as np
from dipy.io.image import load_nifti, save_nifti
from dipy.io.gradients import read_bvals_bvecs
from dipy.core.gradients import gradient_table
import dipy.reconst.dti as dti
from dipy.data import get_fnames
from dipy.segment.mask import median_otsu
from dipy.reconst.dti import fractional_anisotropy, color_fa
from dipy.data import get_sphere
from dipy.viz import window, actor
from matplotlib import cm


# load data
hardi_fname, hardi_bval_fname, hardi_bvec_fname = get_fnames('stanford_hardi')
data, affine = load_nifti(hardi_fname)
bvals, bvecs = read_bvals_bvecs(hardi_bval_fname, hardi_bvec_fname)
gtab = gradient_table(bvals, bvecs)
print('data.shape (%d, %d, %d, %d)' % data.shape)

# speed up by discarding some data
data = data[:, :, 29:35, :]

# mask
maskdata, mask = median_otsu(data, vol_idx=range(10, 50), median_radius=3,
                             numpass=1, autocrop=True, dilate=2)
print('maskdata.shape (%d, %d, %d, %d)' % maskdata.shape)


# fit tensors
print("fitting")
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)

MD2 = tenfit.md
MD2 = (MD2 - MD2.min()) / (MD2.max() - MD2.min())
turbo = cm.get_cmap('turbo', 256)
RGB = turbo(MD2)[..., 0:3]
print(RGB.shape)

sphere = get_sphere('repulsion724')


# Enables/disables interactive visualization
print("plotting")

scene = window.Scene()

evals = tenfit.evals[:, :, 2:3]
evecs = tenfit.evecs[:, :, 2:3]
RGB = RGB[:, :, 2:3, :]

scene.add(actor.tensor_slicer(evals, evecs, scalar_colors=RGB, sphere=sphere,
                              scale=0.3))


[Y, X] = np.meshgrid(range(evals.shape[1]), range(evals.shape[0]))
centers = np.hstack([X.flatten()[:,None], Y.flatten()[:,None]])
centers = np.hstack([centers, np.zeros([centers.shape[0], 1])])

directions_1 = evecs[:, :, :, 0].reshape(-1, 3)
directions_2 = evecs[:, :, :, 1].reshape(-1, 3)
heights = np.ones(centers.shape[0])

axes = actor.axes(scale = (10, 10, 10))
scene.add(axes)

arrow_actor = actor.arrow(centers, directions_1, RGB.reshape(-1, 3), heights)

scene.add(arrow_actor)

window.show(scene)

[samcoveney]

False alarm @skoudoro , I forgot that evecs was shape [Nx, Ny, 1, 3, 3], therefore I was not plotting the eigenvectors after all (I forgot the third dimension of size 1 was there). Combined with updates to master, I have no problem now.

[skoudoro]

I think it has been fixed on #552, #555, #551.

Can you confirm using FURY master?

[samcoveney]

wow these are quite recent fixes so I probably not using them (this discussion should have been mentioned in these fixes). I installed fury with pip install fury or something similar, do you suggest to install from source then?

[samcoveney]

okay, I have just cloned the repo and installed, and the bug remains. Any clues? (See my example code above)

[skoudoro]

no clue, for now, I need to try it. I will have time to do that this afternoon