In this figure, we show the spatial distribution of the performance of the SFM model in predicting the signal.
We start by importing necessary modues:
import tempfile
from IPython.display import Image, display
import nibabel as ni
import osmosis.viz.maya as viz
import osmosis.utils as ozu
import osmosis.model.sparse_deconvolution as ssd
import osmosis.model.analysis as oza
We set the path to the data-files:
import os
import osmosis as oz
import osmosis.io as oio
oio.data_path = os.path.join(oz.__path__[0], 'data')
And read the file-names to be used in the analysis:
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)
A white-matter mask is used (see Methods):
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 have determined that these are good regularization parameters:
# 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)
The axial and radial diffusivity in each b value are read from stored values (calculated as described in Methods).
Then, a model class instance is initialized for every set of data.
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)
Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_1.bvals Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_1.bvecs Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_2.bvals Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_2.bvecs Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.bvals Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.bvecs Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.bvals Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.bvecs Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_1.bvals Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_1.bvecs Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_2.bvals Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_2.bvecs
To provide anatomical context to the display of the results, a T1-weighted was resampled to the DWI data resolution. We load it here:
vol_anat = oio.get_t1(subject, resample=ni.load(oio.data_path + '/%s/%s_wm_mask.nii.gz'%(subject, subject)))
We calculate the rRMSE between the two instances of models (the two independent data sets) in each b value
rrmse_data = [oza.cross_predict(SD_1k_1, SD_1k_2), oza.cross_predict(SD_2k_1, SD_2k_2), oza.cross_predict(SD_4k_1, SD_4k_2)]
SparseDeconvolutionModel.model_params [****************100%******************] 68693 of 68694 complete
/home/arokem/usr/local/lib/python2.7/site-packages/numexpr/necompiler.py:716: DeprecationWarning: using `oa_ndim == 0` when `op_axes` is NULL is deprecated. Use `oa_ndim == -1` or the MultiNew iterator for NumPy <1.8 compatibility return compiled_ex(*arguments, **kwargs) /home/arokem/usr/local/lib/python2.7/site-packages/numexpr/necompiler.py:716: DeprecationWarning: using `oa_ndim == 0` when `op_axes` is NULL is deprecated. Use `oa_ndim == -1` or the MultiNew iterator for NumPy <1.8 compatibility return compiled_ex(*arguments, **kwargs) /home/arokem/usr/local/lib/python2.7/site-packages/numexpr/necompiler.py:716: DeprecationWarning: using `oa_ndim == 0` when `op_axes` is NULL is deprecated. Use `oa_ndim == -1` or the MultiNew iterator for NumPy <1.8 compatibility return compiled_ex(*arguments, **kwargs)
What follows are calls to generate the visualization. This is done using 3-D visualization in Mayavi.
%gui wx
fn = []
for vol_rmse in rrmse_data:
fn.append('%s.png'%tempfile.NamedTemporaryFile().name)
viz.plot_cut_planes(vol_anat,
overlay=vol_rmse,
vmin=0.5,
vmax=1.5,
overlay_cmap="RdYlGn",
invert_cmap=True,
slice_coronal=40,
slice_saggital=15,
slice_axial=30,
view_azim=-40,
view_elev=60,
file_name=fn[-1])
Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x19cbe590> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x1a18b110> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x1a19f470> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x1a1a8c50> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x1a1b8b90> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x1a1c2e90> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x1a201c50> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x1a20dfb0> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x1ae072f0> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x1ae0d5f0> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x1ae198f0> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x1ae25bf0> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x22277350> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x22286650> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x2228c950> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x22295c50> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x2313af50> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x2314e290>
for this_fn in fn:
i = Image(filename=this_fn, width=1280, height=1024)
display(i)
fn = []
for vol_rmse in rrmse_data:
fn.append('%s.png'%tempfile.NamedTemporaryFile().name)
viz.plot_cut_planes(vol_anat,
overlay=vol_rmse,
vmin=0.5,
vmax=1.5,
overlay_cmap="RdYlGn",
invert_cmap=True,
slice_coronal=40,
slice_saggital=15,
slice_axial=45,
view_azim=40,
view_elev=60,
file_name=fn[-1])
Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x231716b0> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x245869b0> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x24593cb0> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x245a0fb0> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x245b52f0> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x25c195f0> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x25c34d10> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x25c40dd0> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x2691a350> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x26926650> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x2692f950> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x26938c50> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x282633b0> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x2826d6b0> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x282799b0> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x28284cb0> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x28290fb0> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x295652f0>
for this_fn in fn:
i = Image(filename=this_fn, width=1280, height=1024)
display(i)
viz.plot_cut_planes(vol_anat,
overlay=rrmse_data[-1],
vmin=0.5,
vmax=1.5,
overlay_cmap="RdYlGn",
invert_cmap=True,
slice_coronal=None,
slice_saggital=None,
slice_axial=45,
view_azim=0,
view_elev=0)
Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x295867d0> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x2958fad0>
<mayavi.core.scene.Scene at 0x2957b590>
%gui wx
viz.plot_cut_planes(vol_anat,
overlay=rrmse_data[-1],
vmin=0.5,
vmax=1.5,
overlay_cmap="RdYlGn",
invert_cmap=True,
slice_coronal=None,
slice_saggital=None,
slice_axial=30,
view_azim=0,
view_elev=0)
Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x2ab881d0> Cannot set extents for <mayavi.modules.image_plane_widget.ImagePlaneWidget object at 0x2ab95530>
<mayavi.core.scene.Scene at 0x2ab828f0>