In [1]:

    

import numpy as np
import nibabel as nib
import osmosis.model.isotropic as mdm
import osmosis.utils as ozu 
import dsitool as dst
import matplotlib 
font = {'family' : 'normal',
        'weight' : 'bold',
        'size'   : 50}
matplotlib.rc('font', **font)

%pylab inline
np.set_printoptions(threshold=np.nan)


    

    




Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.






    




/home/qytian/usr/local/lib/python2.7/site-packages/osmosis/model/base.py:18: UserWarning: Could not import numexpr. Get it! 
  warnings.warn(e_s)

In [2]:

    

files_fit = np.array(['DSI17_exvivo_cc', 'DSI11_invivo_cc'])
files_mean = np.array(['DSI17_exvivo_wm', 'DSI11_invivo_wm'])
models = [mdm.single_exp_rs, mdm.bi_exp_rs, mdm.single_exp_nf_rs]
labels = ['mono-exp', 'bi-exp', 'mono-exp + noise floor']
colors = ['r', 'g', 'c']


    

In [4]:

    

for ii in np.arange(files_fit.shape[0]):
    filename = files_fit[ii]
    
    data = dst.loaddata(filename);
    data = data.astype(float)
    
    gtab = dst.loadgtab(filename);
    bvals = gtab.bvals
    bvecs = gtab.bvecs
    
    mask = np.ones((data.shape[0], data.shape[1], data.shape[2]))
    mask = mask > 0.5
    
    params = []
    
    for isotropic_model, label in zip(models, labels):
        param_out, fit_out, _ = mdm.isotropic_params(data, bvals, bvecs, mask, isotropic_model, signal = "relative_signal")
        params.append(param_out)
    
    filename = files_mean[ii]
    data_wm = dst.loaddata(filename);
    data_wm = data_wm.astype(float);
    data_wm = np.reshape(data_wm, (data_wm.shape[0]*data_wm.shape[1]*data_wm.shape[2], data_wm.shape[3]))
    data_wm_mean = np.mean(data_wm, axis=0)
    uniqmean, uniqbvals = dst.uniquebvals(data_wm_mean, bvals)
    
    fig = figure();
    fig.set_size_inches([8, 6])
    
    sort_idx = np.argsort(bvals)
    plot(bvals[sort_idx]/1000., data[0,0,0][sort_idx]/data[0,0,0,0], 'o')
    for model, p, l, c in zip(models, params, labels, colors):
        plot(bvals[sort_idx]/1000., model(bvals[sort_idx]/1000., *p[0]), color=c, label=l, linewidth=2)    
        
    plot(uniqbvals/1000., uniqmean/uniqmean[0], '-s', color='y', markersize=10, linewidth=2, label='mean signal over WM')
    
    xlim([0, bvals.max()//1000 * 1.05])
    ylim([0, 1.2])
    
    legend()
    
    fitname = 'figure/figure3/' + filename + '_fit.png';
    savefig(fitname, dpi=600)


    

    




None [*****************88%*************     ]  7 of 8 complete