In this notebook, we analyze the accuracy of the diffusion tensor model (DTM).
We start by importing the necessary modules:
import nibabel as ni
import osmosis.model.analysis as oza
import osmosis.model.dti as dti
import osmosis.viz.mpl as mpl
We set the path to point to where the data is installed:
import os
import osmosis as oz
import osmosis.io as oio
oio.data_path = os.path.join(oz.__path__[0], 'data')
We create strings that point to the data for this subject:
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 will only analyze data within a white-matter masked region. Here, we load this mask and generate a masking array:
wm_mask = np.zeros(ni.load(data_1k_1[0]).shape[:3])
wm_nifti = ni.load(oio.data_path + '/%s/%s_wm_mask.nii.gz'%(subject, subject)).get_data()
wm_mask[np.where(wm_nifti==1)] = 1
We generate the model class instances, which will be used to analyze accuracy. Model parameters are not saved to file:
TM_1k_1 = dti.TensorModel(*data_1k_1, mask=wm_mask, params_file='temp')
TM_1k_2 = dti.TensorModel(*data_1k_2, mask=wm_mask, params_file='temp')
TM_2k_1 = dti.TensorModel(*data_2k_1, mask=wm_mask, params_file='temp')
TM_2k_2 = dti.TensorModel(*data_2k_2, mask=wm_mask, params_file='temp')
TM_4k_1 = dti.TensorModel(*data_4k_1, mask=wm_mask, params_file='temp')
TM_4k_2 = dti.TensorModel(*data_4k_2, mask=wm_mask, params_file='temp')
Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_1.bvals Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_1.bvecs Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_2.bvals Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_2.bvecs Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.bvals Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.bvecs Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.bvals Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.bvecs Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_1.bvals Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_1.bvecs Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_2.bvals Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_2.bvecs
We calculate relative RMSE, based on average from cross-prediction between the models fit on each data-set in each b value:
rrmse_1k = oza.cross_predict(TM_1k_1, TM_1k_2)
rrmse_2k = oza.cross_predict(TM_2k_1, TM_2k_2)
rrmse_4k = oza.cross_predict(TM_4k_1, TM_4k_2)
Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_1.nii.gz Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_2.nii.gz Predicting signal from TensorModel Fitting TensorModel params using dipy Predicting signal from TensorModel Fitting TensorModel params using dipy Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.nii.gz Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.nii.gz Predicting signal from TensorModel Fitting TensorModel params using dipy Predicting signal from TensorModel Fitting TensorModel params using dipy Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_1.nii.gz Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_2.nii.gz Predicting signal from TensorModel Fitting TensorModel params using dipy Predicting signal from TensorModel Fitting TensorModel params using dipy
We generate a figure that shows the distribution of the rRMSE measure across the voxels in the white-matter:
# Gray-scale color choices, based on http://colorbrewer2.org
fig = mpl.probability_hist(rrmse_1k[np.isfinite(rrmse_1k)], label='b=%s'%str(1000), color=[0.8, 0.8, 0.8])
fig = mpl.probability_hist(rrmse_2k[np.isfinite(rrmse_2k)], fig=fig, label='b=%s'%str(2000), color=[0.59, 0.59, 0.59])
fig = mpl.probability_hist(rrmse_4k[np.isfinite(rrmse_4k)], fig=fig, label='b=%s'%str(4000), color=[0.32, 0.32, 0.32])
ax = fig.axes[0]
ax.set_xlim([0.6, 1.4])
fig.axes[0].plot([1,1], [ax.get_ylim()[0], ax.get_ylim()[1]], '--k')
fig.axes[0].plot([1/np.sqrt(2),1/np.sqrt(2)], [ax.get_ylim()[0], ax.get_ylim()[1]], '--k')
plt.legend()
fig.savefig('figures/Figure2_histogram.svg')
We also print to the notebook the proportion of voxels with rRMSE that is adequate (smaller than 1):
for this in [rrmse_1k,rrmse_2k, rrmse_4k]:
isfin = this[np.isfinite(this)]
print "The proportion of voxels with rRMSE<1.0 is %s"%(100 * len(np.where(isfin<1)[0])/float(len(isfin)))
The proportion of voxels with rRMSE<1.0 is 98.0485232068 The proportion of voxels with rRMSE<1.0 is 98.8458674609 The proportion of voxels with rRMSE<1.0 is 97.4280218416
We examine a particular voxel. This is used to demonstrate how rRMSE is calculated:
# Look at the same corpus callosum voxel as in Figure 1:
vox_idx = (40, 73, 32) # This doesn't make sense in HT's brain.
First, we scatter-plot data-to-data reliability in this voxel:
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],
color=[0.5, 0.5, 0.5])
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])
fig.savefig('figures/Figure2_scatter_self.svg')
Next, we scatter plot the accuracy of model predictions from one data-set, relative to the data from the other data-set in each b value:
fig, axes = plt.subplots(1,3,squeeze=True)
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],
color=[0.5, 0.5, 0.5])
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}$')
axes[0].set_ylabel(r'Predicted $\frac{S}{S_0}$')
fig.set_size_inches([10,10])
fig.savefig('figures/Figure2_scatter_cross.svg')
print rrmse_1k[vox_idx],rrmse_2k[vox_idx], rrmse_4k[vox_idx]
0.768360252996 0.76873348285 0.863855189387