In the notebook, we examne the statistics of the accuracy of the Sparse Fascicle Model (SFM) across the voxels of the white-matter.
We start by importing the modules we will use:
import nibabel as ni
import osmosis.model.analysis as oza
import osmosis.model.sparse_deconvolution as ssd
import osmosis.model.dti as dti
import osmosis.viz.mpl as mpl
We set the data-path to point to the installation-specific location of the data-files
import os
import osmosis as oz
import osmosis.io as oio
oio.data_path = os.path.join(oz.__path__[0], 'data')
We generate strings pointing to all of the data-files 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)
Wewill only look at white-matter voxels. This loads the white-matter mask:
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
The following are the parameters defining the performance of the Elastic Net solver , based on the calculations in AppendixA
# This is the best according to rRMSE across bvals:
l1_ratio = 0.8
alpha = 0.0005
solver_params = dict(l1_ratio=l1_ratio, alpha=alpha, fit_intercept=False, positive=True)
Initialize the class instances of the SFM objects, using AD/RD derived from the corpus callosum. Parameters are not saved (hence temp
ad_rd = oio.get_ad_rd(subject, 1000)
SD_1k_1 = ssd.SparseDeconvolutionModel(*data_1k_1, mask=wm_mask, params_file ='temp', axial_diffusivity=ad_rd[0]['AD'], radial_diffusivity=ad_rd[0]['RD'], solver_params=solver_params)
SD_1k_2 = ssd.SparseDeconvolutionModel(*data_1k_2, mask=wm_mask, params_file ='temp', axial_diffusivity=ad_rd[1]['AD'], radial_diffusivity=ad_rd[1]['RD'], solver_params=solver_params)
ad_rd = oio.get_ad_rd(subject, 2000)
SD_2k_1 = ssd.SparseDeconvolutionModel(*data_2k_1, mask=wm_mask, params_file ='temp', axial_diffusivity=ad_rd[0]['AD'], radial_diffusivity=ad_rd[0]['RD'], solver_params=solver_params)
SD_2k_2 = ssd.SparseDeconvolutionModel(*data_2k_2, mask=wm_mask, params_file ='temp', axial_diffusivity=ad_rd[1]['AD'], radial_diffusivity=ad_rd[1]['RD'], solver_params=solver_params)
ad_rd = oio.get_ad_rd(subject, 4000)
SD_4k_1 = ssd.SparseDeconvolutionModel(*data_4k_1, mask=wm_mask, params_file ='temp', axial_diffusivity=ad_rd[0]['AD'], radial_diffusivity=ad_rd[0]['RD'], solver_params=solver_params)
SD_4k_2 = ssd.SparseDeconvolutionModel(*data_4k_2, mask=wm_mask, params_file ='temp', axial_diffusivity=ad_rd[1]['AD'], radial_diffusivity=ad_rd[1]['RD'], solver_params=solver_params)
#TM_4k_1 = dti.TensorModel(*data_4k_1, mask=wm_mask)
#TM_4k_2 = dti.TensorModel(*data_4k_2, mask=wm_mask)
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
CallculaterRMSEfor each b value , from a cross-prediction: fit the model to one data-set in each b value and predict the other data set. Calculations average across both directions (1=>2
rrmse_1k = oza.cross_predict(SD_1k_1, SD_1k_2)
rrmse_2k = oza.cross_predict(SD_2k_1, SD_2k_2)
rrmse_4k = oza.cross_predict(SD_4k_1, SD_4k_2)
#rrmse_tensor_4k = oza.cross_predict(TM_4k_1, TM_4k_2)
SparseDeconvolutionModel.model_params [****************100%******************] 64463 of 64464 complete
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])
# Add one of the tensor curves from Figure 2 and put it in the background as reference:
# fig = mpl.probability_hist(rrmse_tensor_4k[np.isfinite(rrmse_4k)], fig=fig, color='gray', label='Tensor model at b=%s'%str(4000))
#fig.set_size_inches([10, 8])
fig.axes[0].plot([1,1], fig.axes[0].get_ylim(), '--k')
fig.axes[0].plot([1/np.sqrt(2),1/np.sqrt(2)], fig.axes[0].get_ylim(), '--k')
fig.axes[0].set_xlim([0.6,1.4])
plt.legend()
fig.savefig('figures/Figure4_histogram.svg')
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.1198808637 The proportion of voxels with rRMSE<1.0 is 99.9084760486 The proportion of voxels with rRMSE<1.0 is 99.8665922065