

In [1]:

    
%pylab inline
from __future__ import print_function, division









    



Populating the interactive namespace from numpy and matplotlib

A notebook to test the scripts to analyze simpace data



In [2]:

    
import nibabel as nib
import numpy as np
import os
import os.path as osp
import matplotlib.pyplot as plt
import nipy
import json
import scipy.linalg as lin
import scipy.stats as sst
from warnings import warn
import datetime, time
import glob as gb
from six import string_types



In [3]:

    
HOME = osp.expanduser('~')
atlas_rois = osp.join(HOME,'code', 'regions', 'regions','tmp')
DDIR = osp.join(HOME,'data','simpace','data','sub001','res_analysis_03')
img_filename = osp.join(DDIR,'waepi_noton.nii' )



In [4]:

    
aal_files = gb.glob(osp.join(atlas_rois, "aal_*"))



In [5]:

    
print(osp.join(atlas_rois, "aal_*"))
print(aal_files[:3])









    



/home/jb/code/regions/regions/tmp/aal_*
['/home/jb/code/regions/regions/tmp/aal_Frontal_Inf_Tri_L___.nii', '/home/jb/code/regions/regions/tmp/aal_Lingual_R___________.nii', '/home/jb/code/regions/regions/tmp/aal_Cerebelum_Crus1_R___.nii']



In [6]:

    
aal = {}
aal['Frontal_Inf_Tri_L'] = nib.load(aal_files[0])
aal_aff = aal['Frontal_Inf_Tri_L'].get_affine()
print(aal_aff)









    



[[  -2.    0.    0.   90.]
 [   0.    2.    0. -126.]
 [   0.    0.    2.  -72.]
 [   0.    0.    0.    1.]]



In [7]:

    
img = nib.load(img_filename)
img_aff = img.get_affine()
print(img_aff)









    



[[  -2.    0.    0.   78.]
 [   0.    2.    0. -112.]
 [   0.    0.    2.  -50.]
 [   0.    0.    0.    1.]]



In [8]:

    
"""
fn = '/home/jb/data/adhd200/NeuroIMAGE/2756846/sfnwmrda2756846_session_1_rest_1.nii.gz'
fn = '/home/jb/data/ds105/sub001/BOLD/task001_run001/bold.nii.gz'
img = nib.load(fn)
print(img.get_affine())
"""









    Out[8]:





"\nfn = '/home/jb/data/adhd200/NeuroIMAGE/2756846/sfnwmrda2756846_session_1_rest_1.nii.gz'\nfn = '/home/jb/data/ds105/sub001/BOLD/task001_run001/bold.nii.gz'\nimg = nib.load(fn)\nprint(img.get_affine())\n"



In [9]:

    
roi = aal['Frontal_Inf_Tri_L']
ijk_roi = (roi.get_data() == 1)
roi.get_data().sum()









    Out[9]:





memmap(2529)



In [10]:

    
print(ijk_roi.shape)
(roi.get_data()[ijk_roi]).shape
(roi.get_data()[ijk_roi])[:10]









    



(91, 109, 91)






    Out[10]:





memmap([1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=int16)



In [11]:

    
roi2img = lin.inv(img_aff).dot(aal_aff)
#lin.inv(aal_aff).dot(img_aff)



In [12]:

    
roi.get_data().shape









    Out[12]:





(91, 109, 91)



In [13]:

    
img.get_data().shape









    Out[13]:





(79, 95, 68, 196)



In [14]:

    
ijk = np.where(ijk_roi)



In [15]:

    
np.asarray(ijk)









    Out[15]:





array([[60, 60, 60, ..., 75, 75, 75],
       [78, 78, 78, ..., 75, 76, 76],
       [35, 36, 37, ..., 47, 42, 43]])



In [16]:

    
coords = np.vstack((ijk[0], ijk[1], ijk[2], np.ones(len(ijk[0]))))
#print(coords.shape)



In [17]:

    
img_ijk = roi2img.dot(coords)
img_ijk.max(axis=1)









    Out[17]:





array([ 69.,  80.,  40.,   1.])



In [18]:

    
np.asarray(img.shape[:3]+(2,)).reshape(4,1) - img_ijk.max(axis=1).reshape(4,1)









    Out[18]:





array([[ 10.],
       [ 15.],
       [ 28.],
       [  1.]])



In [19]:

    
assert np.all(np.logical_and(img_ijk > 0,
              (np.asarray(img.shape[:3]+(2,)).reshape(4,1) - img_ijk.max(axis=1).reshape(4,1) > 0)))



In [19]:



In [33]:

    
import importlib
import utils._utils as ucr
ucr = reload(ucr)



In [22]:

    
HOME = osp.expanduser('~')
atlas_rois = osp.join(HOME,'code', 'regions', 'regions','tmp')
DDIR = osp.join(HOME,'data','simpace','data','sub001','res_analysis_03')
img_filename = osp.join(DDIR,'waepi_noton.nii' )
aal_files = gb.glob(osp.join(atlas_rois, "aal_*"))
aal_file = aal_files[0]



In [23]:

    
sig, issues = ucr.extract_signals(atlas_rois, "aal_*", img_filename, mask=None)









    



[[  -2.    0.    0.   78.]
 [   0.    2.    0. -112.]
 [   0.    0.    2.  -50.]
 [   0.    0.    0.    1.]]
[[  -2.    0.    0.   90.]
 [   0.    2.    0. -126.]
 [   0.    0.    2.  -72.]
 [   0.    0.    0.    1.]]
translation:  [[ -6.]
 [ -7.]
 [-11.]]
could not sample all roi : /home/jb/code/regions/regions/tmp/aal_Cerebelum_8_L_______.nii
[40 19 -5] [61 39 10]
could not sample all roi : /home/jb/code/regions/regions/tmp/aal_Cerebelum_Crus2_R___.nii
[10 10 -1] [38 39 12]
could not sample all roi : /home/jb/code/regions/regions/tmp/aal_Cerebelum_9_L_______.nii
[39 25 -4] [50 37  9]
could not sample all roi : /home/jb/code/regions/regions/tmp/aal_Cerebelum_Crus2_L___.nii
[38 10 -1] [67 36 14]
could not sample all roi : /home/jb/code/regions/regions/tmp/aal_Cerebelum_7b_L______.nii
[40 17 -3] [64 39  8]
could not sample all roi : /home/jb/code/regions/regions/tmp/aal_Cerebelum_9_R_______.nii
[28 25 -5] [38 37  9]
could not sample all roi : /home/jb/code/regions/regions/tmp/aal_Cerebelum_7b_R______.nii
[13 15 -4] [37 39  5]
could not sample all roi : /home/jb/code/regions/regions/tmp/aal_Cerebelum_8_R_______.nii
[16 18 -6] [36 39 10]



In [24]:

    
_keys = sig.keys()
for k in _keys[:10]:
    if sig[k] != None:
        print(sig[k][1])



In [25]:

    
k = _keys[3]
print(k, sig[k][1])
#print(sig[k][0])
if sig[k][0] != None:
    plt.plot(sig[k][0])
    print(issues[k])









    



/home/jb/code/regions/regions/tmp/aal_Cerebelum_10_R______.nii 159
False



In [26]:

    
arr, keys = ucr._dict_signals_to_arr(sig)



In [27]:

    
arr.shape









    Out[27]:





(196, 108)



In [28]:

    
arr_cl, keys_cl = ucr._check_roi_signals(arr, keys)



In [29]:

    
print(set(keys) - set(keys_cl))
print(len(set(keys) - set(keys_cl)))
print(arr_cl.shape)









    



set(['/home/jb/code/regions/regions/tmp/aal_Cerebelum_10_L______.nii', '/home/jb/code/regions/regions/tmp/aal_Vermis_8____________.nii', '/home/jb/code/regions/regions/tmp/aal_Cerebelum_Crus1_L___.nii', '/home/jb/code/regions/regions/tmp/aal_Vermis_7____________.nii', '/home/jb/code/regions/regions/tmp/aal_Cerebelum_Crus1_R___.nii', '/home/jb/code/regions/regions/tmp/aal_Cerebelum_10_R______.nii', '/home/jb/code/regions/regions/tmp/aal_Vermis_9____________.nii', '/home/jb/code/regions/regions/tmp/aal_Vermis_10___________.nii'])
8
(196, 100)



In [56]:

    
xbf, (h,l) = ucr._create_bandpass_bf(196, (0.1, 0.01), 2.0)



In [57]:

    
fig = plt.imshow(xbf, aspect=1, interpolation='nearest')



In [43]:

    
xbf.shape









    Out[43]:





(196, 7)



In [53]:

    
196./dt









    Out[53]:





98.0



In [46]:

    
nt = 196
dt = 2.
tim = np.linspace(0, (nt-1)*dt, nt)
cosBF = ucr._cosine_drift((2*dt), tim)



In [47]:

    
cosBF.shape









    Out[47]:





(196, 196)



In [ ]: