<jupyter_start><jupyter_text>DIPY S0s`02-diff-dipy-S0s.ipynb` ver. 20200422<jupyter_code>import os
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import dipy
from dipy.io.image import save_nifti
from dipy.io.image import load_nifti
from dipy.io import read_bvals_bvecs
from dipy.core.gradients import gradient_table
dipy.__version__
from os.path import expanduser, join
home = expanduser('~')<jupyter_output><empty_output><jupyter_text>Define data and result directories<jupyter_code>dname = join(os.getcwd(), 'data','bids_toro')  # BIDS directory
rname = join(os.getcwd(), 'data','dipy_analysis')  # Results directory<jupyter_output><empty_output><jupyter_text>Preliminary exploration of first subject and first session<jupyter_code>sub = 101
ses = 1

fn = 'sub-%d/ses-%d/dwi/sub-%d_ses-%d_dwi' % (sub, ses, sub, ses)
fdwi = join(dname, '%s.nii.gz' % fn)
fbval = join(dname, '%s.bval' % fn)
fbvec = join(dname, '%s.bvec' % fn)
data, affine, img = load_nifti(fdwi, return_img=True)<jupyter_output><empty_output><jupyter_text>`data` is a 4D array where the first 3 dimensions are the i, j, k voxel coordinates and the last dimension is the number of non-weighted (S0s) and diffusion-weighted volumes.<jupyter_code>print(data.shape)<jupyter_output>(128, 128, 50, 36)<jupyter_text>Check the dimensions of each voxel<jupyter_code>print(img.header.get_zooms()[:3])
bvals, bvecs = read_bvals_bvecs(fbval, fbvec)
# Inspect data

plt.rcParams["figure.figsize"] = (14,7)

b1 = 0
b2 = 10
axial_middle = data.shape[2] // 2
plt.figure('Showing the datasets')
plt.subplot(1, 2, 1).set_axis_off()
plt.imshow(data[:, :, axial_middle, b1].T, cmap='gray', origin='lower')
plt.title('dir=%d: bval=%.0f, bvec=[%.3f,%.3f,%.3f]' % (b1, bvals[b1], bvecs[b1][0],bvecs[b1][1], bvecs[b1][2]), size=14)
plt.subplot(1, 2, 2).set_axis_off()
plt.imshow(data[:, :, axial_middle, b2].T, cmap='gray', origin='lower')
plt.title('dir=%d: bval=%.0f, bvec=[%.3f,%.3f,%.3f]' % (b2, bvals[b2], bvecs[b2][0],bvecs[b2][1], bvecs[b2][2]), size=14)
plt.suptitle('DWI: sub-%d, ses-%d, slice no %d/%d' % (sub, ses, axial_middle, data.shape[2]), size=16)
plt.savefig('./assets/dipy_data.png', bbox_inches='tight')
plt.show()<jupyter_output><empty_output><jupyter_text>In DIPY, we use an object called `GradientTable` which holds all the acquisition specific parameters, e.g. b-values, b-vectors, timings and others.<jupyter_code>gtab = gradient_table(bvals, bvecs)
print(gtab.bvals)
print(gtab.bvecs[:11, :])<jupyter_output>[[ 0.     0.     0.   ]
 [ 0.     0.     0.   ]
 [ 0.     0.     0.   ]
 [ 0.     0.     0.   ]
 [ 0.     0.     0.   ]
 [ 0.     0.     0.   ]
 [ 1.    -0.    -0.   ]
 [ 0.166 -0.986 -0.   ]
 [-0.11  -0.664 -0.74 ]
 [ 0.901  0.419  0.11 ]
 [-0.169  0.601 -0.781]]<jupyter_text>`gtab` can be used to tell what part of the data is the S0 volumes (volumes which correspond to b-values of 0)<jupyter_code>S0s = data[:, :, :, gtab.b0s_mask]
print(S0s.shape)
print(S0s.shape[-1])<jupyter_output>(128, 128, 50, 6)
6<jupyter_text>We save this in a new Nifti file<jupyter_code>fnS0s = join(rname, 'sub_%d_tp%d_dwi_S0.nii.gz' % (sub, ses))
save_nifti(fnS0s, S0s, affine)<jupyter_output><empty_output><jupyter_text>Define function `extract_S0s`<jupyter_code>from dipy.io.image import load_nifti
from dipy.io import read_bvals_bvecs
from dipy.core.gradients import gradient_table
from dipy.io.image import save_nifti

def extract_S0s(inpdir, sub, ses, outdir):
    """
    Extract S0s from dwi file
    Usage:
    Ex: 
    """
    
    fn = 'sub-%d/ses-%d/dwi/sub-%d_ses-%d_dwi' % (sub, ses, sub, ses)
    fdwi = join(inpdir, '%s.nii.gz' % fn)
    #print(fdwi)
    fbval = join(inpdir, '%s.bval' % fn)
    fbvec = join(inpdir, '%s.bvec' % fn)
    
    bvals, bvecs = read_bvals_bvecs(fbval, fbvec)
    gtab = gradient_table(bvals, bvecs)
    data, affine, img = load_nifti(fdwi, return_img=True)
    S0s = data[:, :, :, gtab.b0s_mask]
    nS0s = S0s.shape[-1]
    
    fnS0s = join(outdir, 'sub_%d_tp%d_dwi_S0.nii.gz' % (sub, ses))
    save_nifti(fnS0s, S0s, affine)
    
    return S0s, affine, nS0s
# Testing function for single subject, single examination
from os.path import expanduser, join
home = expanduser('~')
inpdir = join(home, 'prj','BrainGut','GE_Signa3T','bids_toro')
sub = 101
ses = 1
outdir = join(home, 'prj','BrainGut','GE_Signa3T','dipy_analysis')

S0s, affine, nS0s = extract_S0s(inpdir, sub, ses, outdir)
print('affine:\n', affine)
print('nS0s:', nS0s)
S0s_mean = S0s.mean(axis=-1)
nS0s = S0s.shape[-1]
print('S0s.shape:', S0s.shape)
print('S0s_mean.shape:', S0s_mean.shape)
print('nS0s:', nS0s)

fnS0s_mean = join(outdir, 'sub_%d_tp%d_dwi_S0_mean_of_%d.nii.gz' % (sub, ses, nS0s))
save_nifti(fnS0s_mean, S0s_mean, affine)
# !fsleyes ~/prj/BrainGut/GE_Signa3T/dipy_analysis/sub_101_tp1_dwi_S0.nii.gz<jupyter_output><empty_output><jupyter_text>Run for all subjects<jupyter_code>from os.path import expanduser, join
import sys
sys.path.append('.')   # path to utils.py
import utils
import time

inpdir = dname
ses = 1
outdir = rname
%%time
nsubj = 1 # 30
ses=1
for i in utils.progressbar(range(nsubj), "Computing: ", 40):
    time.sleep(0.1)
    sub = i+101
    S0s, affine, nS0s = extract_S0s(inpdir, sub, ses, outdir)
    
    S0s_mean = S0s.mean(axis=-1)
    nS0s = S0s.shape[-1]
    fnS0s_mean = join(outdir, 'sub_%d_tp%d_dwi_S0_mean.nii.gz' % (sub, ses))
    save_nifti(fnS0s_mean, S0s_mean, affine)<jupyter_output>Computing: [########################################] 1/1
CPU times: user 6.1 s, sys: 127 ms, total: 6.23 s
Wall time: 698 ms
{'fraction_non_alphanumeric': 0.16891527253103075, 'fraction_numerical': 0.04515200575643101, 'mean_word_length': 4.571142284569138, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 37, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '69205', 'license': 'MIT', 'n_tokens_mistral': 2588, 'n_tokens_neox': 2442, 'n_words': 496, 'path': 'diffusion/02-diff-dipy-S0s.ipynb', 'repo_name': 'computational-medicine/BMED360-2020', 'size': '5559'}