notebook.community

Edit and run



In [1]:

    
import numpy as np
import nibabel as nib
import dipy.reconst.dti as dti
import dsiadapt as dsi
import dipy.core.gradients as grad
import matplotlib 
%pylab inline
np.set_printoptions(threshold=numpy.nan)









    



Populating the interactive namespace from numpy and matplotlib



In [2]:

    
# Estimate diffusion rates using DTI fitting

fnArr = np.array(['DSI17_exvivo', 'DSI11_invivo_b10k', 'DSI11_invivo_b7k'])
bmaxthrArr = np.array([4500, 2000, 2000]);
adArr = np.zeros(fnArr.shape[0]); # axial diffusion rate

for ii in np.arange(fnArr.shape[0]):
    fn = fnArr[ii];
    data = nib.load('data/' + fn + '_cc.nii.gz').get_data();
    gtab = grad.gradient_table('data/' + fn + '_bvals.txt', 'data/' + fn + '_bvecs.txt');
    mask = gtab.bvals <= bmaxthrArr[ii]; # Use only data with bvals less than 3000
    datamask = data[:, :, :, mask]
    gtabmask = grad.gradient_table_from_bvals_bvecs(gtab.bvals[mask], gtab.bvecs[mask, :])
    disp(gtabmask.bvals.shape) # Number of measurements used for DTI fitting
    disp(gtabmask.bvals.max())
    
    tensormodel = dti.TensorModel(gtabmask)
    tensorfit = tensormodel.fit(datamask)
    adArr[ii] = tensorfit.ad.mean()

adArr = np.concatenate((adArr[0:1], adArr[0:1], adArr)) # Use the ad from ex vivo DSI17 for ex vivo DSI15 and ex vivo DSI11



In [3]:

    
disp(adArr)









    



[ 0.00017951  0.00017951  0.00017951  0.00142466  0.00163717]



In [44]:

    
tensorfit.evecs









    Out[44]:





array([[[[[-0.96167792,  0.20693869, -0.17986648],
          [ 0.26721999,  0.56050859, -0.78385177],
          [ 0.06139255,  0.80187686,  0.5943269 ]],

         [[-0.97242061,  0.13196121, -0.1923133 ],
          [ 0.23205059,  0.6303544 , -0.74081432],
          [-0.02346678,  0.76500952,  0.64359128]]]],



       [[[[-0.98296716,  0.16020355, -0.0900577 ],
          [ 0.15738741,  0.4807815 , -0.86259977],
          [ 0.09489347,  0.86208119,  0.49780644]],

         [[-0.99204213,  0.09891501, -0.0778989 ],
          [ 0.11463644,  0.45375531, -0.88372202],
          [ 0.05206634,  0.88561952,  0.46148364]]]],



       [[[[-0.99974284, -0.02043197,  0.00983773],
          [-0.01676935,  0.37407194, -0.92724806],
          [-0.01526549,  0.92717459,  0.37431837]],

         [[-0.99896407,  0.04293378, -0.01508249],
          [ 0.03005452,  0.37360834, -0.92709953],
          [ 0.03416894,  0.92659241,  0.37451167]]]],



       [[[[-0.98087337,  0.10256257, -0.16543385],
          [-0.1610008 ,  0.05014837,  0.9856794 ],
          [ 0.10939005,  0.99346165, -0.03267654]],

         [[-0.95518342,  0.17954174, -0.23534952],
          [-0.22104474,  0.09619154,  0.97050833],
          [ 0.19688539,  0.97903624, -0.05219381]]]]])



In [4]:

    
# Estimate mean diffusion distance
fnArr = np.array(['DSI11_exvivo', 'DSI15_exvivo', 'DSI17_exvivo', 'DSI11_invivo_b10k', 'DSI11_invivo_b7k'])
DELTAArr = np.array([29.4 * 10**-3, 29.4 * 10**-3, 29.4 * 10**-3, 20.9 * 10**-3, 49.2 * 10**-3]); # in second
deltaArr = np.array([16.7 * 10**-3, 16.7 * 10**-3, 16.7 * 10**-3, 12.9 * 10**-3, 42.3 * 10**-3]); # in second
fovArr = np.zeros(fnArr.shape[0]);
mddArr = np.zeros(fnArr.shape[0]);

for ii in np.arange(fnArr.shape[0]):
    fn = fnArr[ii];
    gtab = grad.gradient_table('data/' + fn + '_bvals.txt', 'data/' + fn + '_bvecs.txt');
    DELTA = DELTAArr[ii]; # second
    delta = deltaArr[ii]; # second
    ad = adArr[ii]; # mm2/s
    
    # Compute FOV
    bmax = gtab.bvals.max(); # s/mm2
    qmax = np.sqrt(bmax / (DELTA - delta / 3)) / (2 * np.pi); # mm-1
    deltaq = qmax / dsi.create_qtable(gtab).max(); # mm-1
    fovArr[ii] = 1 / deltaq; # mm
    
    # Compute MDD
    mddArr[ii] = np.sqrt(6 * ad * (DELTA - delta / 3)); # mm
    
    np.savetxt('data/' + fn + '_stats.txt', (mddArr[ii], fovArr[ii]));



In [5]:

    
disp(mddArr)
disp(fovArr)









    



[ 0.00506661  0.00506661  0.00506661  0.01191204  0.01856845]
[ 0.02797822  0.0391695   0.04476515  0.04047659  0.07034843]



In [51]:

    
# Estimate signal at the q-space edge
fnArr = np.array(['DSI17_exvivo', 'DSI11_invivo_b10k', 'DSI11_invivo_b7k'])
bmaxthrArr = np.array([4500, 2000, 2000]);
adArr = np.zeros(fnArr.shape[0]); # axial diffusion rate

for ii in np.arange(fnArr.shape[0]):
    fn = fnArr[ii];
    data = nib.load('data/' + fn + '_cc.nii.gz').get_data();
    gtab = grad.gradient_table('data/' + fn + '_bvals.txt', 'data/' + fn + '_bvecs.txt');
    ori = np.array([0, 1, 0]);
    mask = np.logical_or(np.dot(gtab.bvecs, ori) == 1, np.dot(gtab.bvecs, ori) == -1);
    mask[0] = True
    oribvals = gtab.bvals[mask];
    oridata = data[:, :, :, mask] * 1.0;
    oridatasum = oridata.sum(axis=0)
    oridatasum = np.squeeze(oridatasum.sum(axis=1))
    oridatasum = oridatasum / oridatasum[0];
    disp(oridatasum)
    disp(oribvals)









    



[ 1.          0.9682532   0.88987076  0.78094232  0.7204265   0.64146775
  0.63224626  0.67720091  0.56481433  0.50314587  0.50718021  0.42072907
  0.42533979  0.46971807  0.41611832  0.3290908   0.45127514]
[     0.    450.    450.   1900.   1900.   4200.   4200.   7500.   7500.
  11700.  11700.  16900.  16900.  22950.  22950.  30000.  30000.]
[ 1.          0.82066344  0.92743607  0.61575674  0.65134762  0.46060815
  0.4975812   0.38908086  0.40255701  0.2950933   0.34208708]
[     0.    400.    400.   1600.   1600.   3600.   3600.   6400.   6400.
  10000.  10000.]
[ 1.          0.86217952  0.72875309  0.71066624  0.67684025  0.53699291
  0.54359198  0.33684531  0.44082347  0.41415536  0.39524052]
[    0.   280.   280.  1120.  1120.  2520.  2520.  4480.  4480.  7000.
  7000.]



In [48]:

    
(0.30603716 + 0.29278132 + 0.39364102 + 0.28989962) / 4









    Out[48]:





0.32058978000000005



In [49]:

    
(0.06185211 + 0.04872149) / 2









    Out[49]:





0.0552868



In [50]:

    
(0.10041282 + 0.12508167) / 2









    Out[50]:





0.112747245



In [52]:

    
(0.46971807 + 0.41611832 + 0.3290908 + 0.45127514) / 4









    Out[52]:





0.4165505825



In [53]:

    
(0.2950933 + 0.34208708) / 2









    Out[53]:





0.31859019



In [54]:

    
(0.41415536 + 0.39524052) / 2









    Out[54]:





0.40469794