import nibabel as ni
import osmosis as oz
import osmosis.model.analysis as oza
import osmosis.model.dti as dti
import osmosis.viz.mpl as mpl
import osmosis.io as oio
oio.data_path = '/biac4/wandell/biac2/wandell6/data/arokem/osmosis'
subject = 'FP'
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)
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
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: /biac4/wandell/biac2/wandell6/data/arokem/osmosis/FP/0009_01_DWI_2mm150dir_2x_b1000_aligned_trilin.bvals Loading from file: /biac4/wandell/biac2/wandell6/data/arokem/osmosis/FP/0009_01_DWI_2mm150dir_2x_b1000_aligned_trilin.bvecs Loading from file: /biac4/wandell/biac2/wandell6/data/arokem/osmosis/FP/0011_01_DWI_2mm150dir_2x_b1000_aligned_trilin.bvals Loading from file: /biac4/wandell/biac2/wandell6/data/arokem/osmosis/FP/0011_01_DWI_2mm150dir_2x_b1000_aligned_trilin.bvecs Loading from file: /biac4/wandell/biac2/wandell6/data/arokem/osmosis/FP/0005_01_DTI_2mm_150dir_2x_b2000_aligned_trilin.bvals Loading from file: /biac4/wandell/biac2/wandell6/data/arokem/osmosis/FP/0005_01_DTI_2mm_150dir_2x_b2000_aligned_trilin.bvecs Loading from file: /biac4/wandell/biac2/wandell6/data/arokem/osmosis/FP/0007_01_DTI_2mm_150dir_2x_b2000_aligned_trilin.bvals Loading from file: /biac4/wandell/biac2/wandell6/data/arokem/osmosis/FP/0007_01_DTI_2mm_150dir_2x_b2000_aligned_trilin.bvecs Loading from file: /biac4/wandell/biac2/wandell6/data/arokem/osmosis/FP/0005_01_DWI_2mm150dir_2x_b4000_aligned_trilin.bvals Loading from file: /biac4/wandell/biac2/wandell6/data/arokem/osmosis/FP/0005_01_DWI_2mm150dir_2x_b4000_aligned_trilin.bvecs Loading from file: /biac4/wandell/biac2/wandell6/data/arokem/osmosis/FP/0007_01_DWI_2mm150dir_2x_b4000_aligned_trilin.bvals Loading from file: /biac4/wandell/biac2/wandell6/data/arokem/osmosis/FP/0007_01_DWI_2mm150dir_2x_b4000_aligned_trilin.bvecs
fig = mpl.probability_hist(TM_1k_1.mode[np.isfinite(TM_1k_1.mode)])
idx = np.isfinite(TM_1k_1.mode)
fig = mpl.scatter_density(TM_1k_1.mode[idx], TM_1k_1.fractional_anisotropy[idx])
/home/arokem/usr/lib/python2.7/site-packages/osmosis/viz/mpl.py:353: RuntimeWarning: divide by zero encountered in log10 imax = ax.matshow(np.log10(np.flipud(data_arr.T)), cmap=cmap, **kwargs)
fig = mpl.mosaic(TM_1k_1.mode.T, cmap=matplotlib.cm.hot, vmin=-1)
fig.set_size_inches([20,10])
pdd_rel_1k = oza.pdd_reliability(TM_1k_1, TM_1k_2)
pdd_rel_2k = oza.pdd_reliability(TM_2k_1, TM_2k_2)
pdd_rel_4k = oza.pdd_reliability(TM_4k_1, TM_4k_2)
Loading from file: /biac4/wandell/biac2/wandell6/data/arokem/osmosis/FP/0009_01_DWI_2mm150dir_2x_b1000_aligned_trilin.nii.gz Fitting TensorModel params using dipy Loading from file: /biac4/wandell/biac2/wandell6/data/arokem/osmosis/FP/0011_01_DWI_2mm150dir_2x_b1000_aligned_trilin.nii.gz Fitting TensorModel params using dipy Loading from file: /biac4/wandell/biac2/wandell6/data/arokem/osmosis/FP/0005_01_DTI_2mm_150dir_2x_b2000_aligned_trilin.nii.gz Fitting TensorModel params using dipy Loading from file: /biac4/wandell/biac2/wandell6/data/arokem/osmosis/FP/0007_01_DTI_2mm_150dir_2x_b2000_aligned_trilin.nii.gz Fitting TensorModel params using dipy Loading from file: /biac4/wandell/biac2/wandell6/data/arokem/osmosis/FP/0005_01_DWI_2mm150dir_2x_b4000_aligned_trilin.nii.gz Fitting TensorModel params using dipy Loading from file: /biac4/wandell/biac2/wandell6/data/arokem/osmosis/FP/0007_01_DWI_2mm150dir_2x_b4000_aligned_trilin.nii.gz Fitting TensorModel params using dipy
For b=1000, the median PDD reliability is 3.62 For b=2000, the median PDD reliability is 3.80 For b=4000, the median PDD reliability is 5.01
# For comparison with Jones 2003, Figure 3:
fig = mpl.scatter_density(np.hstack([0,1, TM_1k_1.linearity[np.isfinite(pdd_rel_1k)]]), np.hstack([0,1, pdd_rel_1k[np.isfinite(pdd_rel_1k)]]))
fig = mpl.scatter_density(np.hstack([0,1, TM_2k_1.linearity[np.isfinite(pdd_rel_2k)]]), np.hstack([0,1, pdd_rel_2k[np.isfinite(pdd_rel_2k)]]))
fig = mpl.scatter_density(np.hstack([0,1, TM_4k_1.linearity[np.isfinite(pdd_rel_4k)]]), np.hstack([0,1, pdd_rel_4k[np.isfinite(pdd_rel_4k)]]))
/home/arokem/usr/lib/python2.7/site-packages/osmosis/viz/mpl.py:350: RuntimeWarning: divide by zero encountered in log10 imax = ax.matshow(np.log10(np.flipud(data_arr.T)), cmap=cmap)
fit_rel_1k = oza.fit_reliability(TM_1k_1, TM_1k_2)
fit_rel_2k = oza.fit_reliability(TM_2k_1, TM_2k_2)
fit_rel_4k = oza.fit_reliability(TM_4k_1, TM_4k_2)
Predicting signal from TensorModel Predicting signal from TensorModel Predicting signal from TensorModel Predicting signal from TensorModel Predicting signal from TensorModel Predicting signal from TensorModel
fig = mpl.probability_hist(fit_rel_1k[np.isfinite(fit_rel_1k)], label='b=%s'%str(1000))
fig = mpl.probability_hist(fit_rel_2k[np.isfinite(fit_rel_2k)], fig=fig, label='b=%s'%str(2000))
fig = mpl.probability_hist(fit_rel_4k[np.isfinite(fit_rel_4k)], fig=fig, label='b=%s'%str(4000))
ax = fig.axes[0]