This notebook displays the signal in one voxel across different b values. We also demonstrate the fit of the tensor model to the signal and the process of predicting back the signal from the model. This process is at the core of our cross-validation approach.
We import the modules we will use:
import numpy as np
import os
import tempfile
from IPython.display import Image, display
import nibabel as ni
import osmosis.model.dti as dti
import osmosis.viz.maya as maya
import osmosis.viz.mpl as mpl
import osmosis.utils as ozu
import osmosis.volume as ozv
import osmosis.tensor as ozt
We set the data path to the installation-specific location:
import os
import osmosis as oz
import osmosis.io as oio
oio.data_path = os.path.join('/biac4/wandell/data/osmosis/')
Get data file-names:
subject = 'SUB1'
data_1k_1, data_1k_2 = oio.get_dwi_data(1000, subject)
data_2k_1, data_2k_2 = oio.get_dwi_data(2000, subject)
data_4k_1, data_4k_2 = oio.get_dwi_data(4000, subject)
We choose a voxel by its XYZ coordinates.
To use only this voxel, we define a mask around that voxel, which we will pass to mask the model class instances
# A corpus callosum voxel:
vox_idx = (40, 74, 34)
mask = np.zeros(ni.load(data_1k_1[0]).shape[:3])
mask[vox_idx[0]-1:vox_idx[0]+1,
vox_idx[1]-1:vox_idx[1]+1,
vox_idx[2]-1:vox_idx[2]+1]=1
We initialize TensorModel
class instances
TM_1k_1 = dti.TensorModel(*data_1k_1, mask=mask, params_file='temp')
TM_1k_2 = dti.TensorModel(*data_1k_2, mask=mask, params_file='temp')
TM_2k_1 = dti.TensorModel(*data_2k_1, mask=mask, params_file='temp')
TM_2k_2 = dti.TensorModel(*data_2k_2, mask=mask, params_file='temp')
TM_4k_1 = dti.TensorModel(*data_4k_1, mask=mask, params_file='temp')
TM_4k_2 = dti.TensorModel(*data_4k_2, mask=mask, params_file='temp')
Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0009_01_DWI_2mm150dir_2x_b1000_aligned_trilin.bvals Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0009_01_DWI_2mm150dir_2x_b1000_aligned_trilin.bvecs Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0011_01_DWI_2mm150dir_2x_b1000_aligned_trilin.bvals Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0011_01_DWI_2mm150dir_2x_b1000_aligned_trilin.bvecs Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0005_01_DTI_2mm_150dir_2x_b2000_aligned_trilin.bvals Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0005_01_DTI_2mm_150dir_2x_b2000_aligned_trilin.bvecs Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0007_01_DTI_2mm_150dir_2x_b2000_aligned_trilin.bvals Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0007_01_DTI_2mm_150dir_2x_b2000_aligned_trilin.bvecs Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0005_01_DWI_2mm150dir_2x_b4000_aligned_trilin.bvals Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0005_01_DWI_2mm150dir_2x_b4000_aligned_trilin.bvecs Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0007_01_DWI_2mm150dir_2x_b4000_aligned_trilin.bvals Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0007_01_DWI_2mm150dir_2x_b4000_aligned_trilin.bvecs
The next few notebook cells use 3D visualization techniques from our mayavi module. We plot an interpolation of the measured signal in each of the b values
%gui wx
# Interpolate the signal from the measurement b-vectors using radial basis functions and display in 3D:
fig = maya.plot_signal_interp(TM_1k_1.bvecs[:,TM_1k_1.b_idx], TM_1k_1.relative_signal[vox_idx], cmap='hot', origin=[0,0,0], vmin=0, vmax=1)
fig = maya.plot_signal_interp(TM_2k_1.bvecs[:,TM_2k_1.b_idx], TM_2k_1.relative_signal[vox_idx], cmap='hot', origin=[0,0,2], vmin=0, vmax=1, figure=fig, offset=-2)
# Put it into a temporary file:
fn = '%s.png'%tempfile.NamedTemporaryFile().name
#orientation_axes = mlab.orientation_axes(xlabel='', ylabel='').axes
fig = maya.plot_signal_interp(TM_4k_1.bvecs[:,TM_4k_1.b_idx], TM_4k_1.relative_signal[vox_idx], cmap='hot', origin=[0,0,4], colorbar=True, vmin=0, vmax=1, figure=fig, offset=-4, file_name=fn, roll=90)
/home/arokem/usr/lib/python2.7/site-packages/osmosis/model/base.py:383: RuntimeWarning: divide by zero encountered in divide signal_rel = self.signal/np.reshape(self.S0, (self.S0.shape + (1,))) /home/arokem/usr/lib/python2.7/site-packages/osmosis/model/base.py:383: RuntimeWarning: invalid value encountered in divide signal_rel = self.signal/np.reshape(self.S0, (self.S0.shape + (1,))) /usr/lib/python2.7/dist-packages/scipy/interpolate/rbf.py:122: RuntimeWarning: divide by zero encountered in log result = r**2 * log(r)
Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0009_01_DWI_2mm150dir_2x_b1000_aligned_trilin.nii.gz Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0005_01_DTI_2mm_150dir_2x_b2000_aligned_trilin.nii.gz Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0005_01_DWI_2mm150dir_2x_b4000_aligned_trilin.nii.gz
/usr/lib/python2.7/dist-packages/scipy/interpolate/rbf.py:122: RuntimeWarning: invalid value encountered in multiply result = r**2 * log(r)
i = Image(filename=fn)
display(i)
The following cells are commented out per default, but they can be useful if you want to make animations interrogating the signal from different angles:
# Make figures for a movie:
#for ang in range(360):
# fn = '/home/arokem/Dropbox/osmosis_paper_figures/movie_figures_sig_2k/frame%03d.png'%ang
# if not os.path.exists(fn):
# fig = maya.plot_signal_interp(TM_2k_1.bvecs[:,TM_2k_1.b_idx], TM_2k_1.relative_signal[vox_idx], cmap='hot', colorbar=True, vmin=0, vmax=1, elevation=ang, file_name=fn)
#for ang in range(360):
# fn = '/home/arokem/Dropbox/osmosis_paper_figures/movie_figures_sig_2k_w_pts/frame%03d.png'%ang
# if not os.path.exists(fn):
# fig = maya.plot_signal_interp(TM_2k_1.bvecs[:,TM_2k_1.b_idx], TM_2k_1.relative_signal[vox_idx], cmap='hot', colorbar=True, vmin=0, vmax=1, elevation=ang, file_name=fn, points=True,scale_points=False)
This is the model prediction
fig = maya.plot_signal_interp(TM_1k_1.bvecs[:,TM_1k_1.b_idx], TM_1k_1.fit[vox_idx]/TM_1k_1.S0[vox_idx], origin=[0,0,0], cmap='hot', vmin=0, vmax=1)
fig = maya.plot_signal_interp(TM_2k_1.bvecs[:,TM_2k_1.b_idx], TM_2k_1.fit[vox_idx]/TM_1k_1.S0[vox_idx], origin=[0,0,2], cmap='hot', figure=fig, offset=-2, vmin=0, vmax=1)
# Put it into a temporary file:
fn = '%s.png'%tempfile.NamedTemporaryFile().name
#orientation_axes = mlab.orientation_axes(xlabel='', ylabel='').axes
fig = maya.plot_signal_interp(TM_4k_1.bvecs[:,TM_4k_1.b_idx], TM_4k_1.fit[vox_idx]/TM_4k_1.S0[vox_idx], origin=[0,0,4], cmap='hot', colorbar=True, figure=fig, offset=-4, file_name=fn, roll=90, vmin=0, vmax=1)
Predicting signal from TensorModel Fitting TensorModel params using dipy Predicting signal from TensorModel Fitting TensorModel params using dipy Predicting signal from TensorModel Fitting TensorModel params using dipy
i = Image(filename=fn)
display(i)
The following is a visualization of the ADC surface of the tensor model in the same coordinate frame
ten_1k = ozt.tensor_from_eigs(TM_1k_1.model_params[vox_idx][3:].reshape(3,3),
TM_1k_1.model_params[vox_idx][:3],
TM_1k_1.bvecs[:, TM_1k_1.b_idx], TM_1k_1.bvals[TM_1k_1.b_idx])
fig = maya.plot_tensor_3d(ten_1k, mode='ADC', cmap='hot', origin=[0,0,0], colorbar=True, vmin=0, vmax=2.5)
ten_2k = ozt.tensor_from_eigs(TM_2k_1.model_params[vox_idx][3:].reshape(3,3),
TM_2k_1.model_params[vox_idx][:3],
TM_2k_1.bvecs[:, TM_2k_1.b_idx], TM_2k_1.bvals[TM_2k_1.b_idx])
fig = maya.plot_tensor_3d(ten_2k, mode='ADC', cmap='hot', origin=[0,0,2], vmin=0, vmax=2.5, figure=fig, offset=-2.5)
ten_4k = ozt.tensor_from_eigs(TM_4k_1.model_params[vox_idx][3:].reshape(3,3),
TM_4k_1.model_params[vox_idx][:3],
TM_4k_1.bvecs[:, TM_4k_1.b_idx], TM_4k_1.bvals[TM_4k_1.b_idx])
fn = '%s.png'%tempfile.NamedTemporaryFile().name
fig = maya.plot_tensor_3d(ten_4k, mode='ADC', cmap='hot', origin=[0,0,4], vmin=0, vmax=2.5, figure=fig, offset=-5, file_name=fn, roll=90)
i = Image(filename=fn, width=1280, height=1024)
display(i)
We can also display the ellipsoid mean diffusion distance implied by the model
ten_1k = ozt.tensor_from_eigs(TM_1k_1.model_params[vox_idx][3:].reshape(3,3),
TM_1k_1.model_params[vox_idx][:3],
TM_1k_1.bvecs[:, TM_1k_1.b_idx], TM_1k_1.bvals[TM_1k_1.b_idx])
fig = maya.plot_tensor_3d(ten_1k, mode='ellipse', cmap='Greens', colorbar=True)
ten_2k = ozt.tensor_from_eigs(TM_2k_1.model_params[vox_idx][3:].reshape(3,3),
TM_2k_1.model_params[vox_idx][:3],
TM_2k_1.bvecs[:, TM_2k_1.b_idx], TM_2k_1.bvals[TM_2k_1.b_idx])
fig = maya.plot_tensor_3d(ten_2k, mode='ellipse', cmap='Greens', offset=-2.5, figure=fig)
ten_4k = ozt.tensor_from_eigs(TM_4k_1.model_params[vox_idx][3:].reshape(3,3),
TM_4k_1.model_params[vox_idx][:3],
TM_4k_1.bvecs[:, TM_4k_1.b_idx], TM_4k_1.bvals[TM_4k_1.b_idx])
fn = '%s.png'%tempfile.NamedTemporaryFile().name
fig = maya.plot_tensor_3d(ten_4k, mode='ellipse', cmap='Greens', offset=-5, figure=fig, file_name=fn, roll=90)
/home/arokem/usr/lib/python2.7/site-packages/osmosis/tensor.py:248: RuntimeWarning: invalid value encountered in sqrt dist = np.diag(1 / np.sqrt(sphADC))
i = Image(filename=fn, width=1280, height=1024)
display(i)
In what follows, we compare the signal in two separate measurements
fig, axes = plt.subplots(1,3,squeeze=True)
axes[0].set_ylabel(r'$\frac{S}{S_0}$')
for ax, models in zip(axes,([TM_1k_1,TM_1k_2],[TM_2k_1, TM_2k_2],[TM_4k_1, TM_4k_2])):
ax.scatter(models[0].relative_signal[vox_idx], models[1].relative_signal[vox_idx])
ax.grid()
ax.set_xlim([0,1])
ax.set_ylim([0,1])
ax.set_aspect('equal')
ax.set_xlabel(r'$\frac{S}{S_0}$')
fig.set_size_inches([10,10])
Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0011_01_DWI_2mm150dir_2x_b1000_aligned_trilin.nii.gz Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0007_01_DTI_2mm_150dir_2x_b2000_aligned_trilin.nii.gz Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0007_01_DWI_2mm150dir_2x_b4000_aligned_trilin.nii.gz
In comparison, here is the model prediction relative to the signal in the other measurement. Relative RMSE is defined to be the ratio between the noise in these two displays of the data.
fig, axes = plt.subplots(1,3,squeeze=True)
axes[0].set_ylabel(r'Predicted $\frac{S}{S_0}$')
for ax, models in zip(axes,([TM_1k_1,TM_1k_2],[TM_2k_1, TM_2k_2],[TM_4k_1, TM_4k_2])):
ax.scatter(models[0].relative_signal[vox_idx], models[1].fit[vox_idx]/models[1].S0[vox_idx])
ax.grid()
ax.set_xlim([0,1])
ax.set_ylim([0,1])
ax.set_aspect('equal')
ax.set_xlabel(r'Measured $\frac{S}{S_0}$')
fig.set_size_inches([10,10])
Predicting signal from TensorModel Fitting TensorModel params using dipy Predicting signal from TensorModel Fitting TensorModel params using dipy Predicting signal from TensorModel Fitting TensorModel params using dipy