This notebook separates the annotated JFRC template into masks



In [1]:

    
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy import io
%matplotlib inline 
import pylab
import nibabel as nb
import numpy as np
import scipy.ndimage

Open Data template



In [3]:

    
# from http://stackoverflow.com/questions/3579568/choosing-a-file-in-python-with-simple-dialog
from Tkinter import Tk
from tkFileDialog import askopenfilename

Tk().withdraw() # we don't want a full GUI, so keep the root window from appearing
filename = askopenfilename() # show an "Open" dialog box and return the path to the selected file
print(filename)
img = nb.load(filename)
data = img.get_data()
S=data.shape









    



/home/sophie/Desktop/Registration/862863rby2/862863Rby2MAX_combined.nii

Open Masks



In [4]:

    
# from http://stackoverflow.com/questions/3579568/choosing-a-file-in-python-with-simple-dialog
from Tkinter import Tk
from tkFileDialog import askopenfilename

Tk().withdraw() # we don't want a full GUI, so keep the root window from appearing
filenameM = askopenfilename() # show an "Open" dialog box and return the path to the selected file
print(filenameM)
img1 = nb.load(filenameM)
Masks = img1.get_data()
Sm=Masks.shape
Masks=np.array(Masks)









    



/home/sophie/Desktop/Registration/862863rby2/jrfc862863rby2Transformed.nii

Open Nomenclature



In [5]:

    
filenameM='/home/sophie/Desktop/Registration/RegionList'
with open(filenameM) as f:
    content = f.readlines()
Names=[Line.split('\t') for Line in content]
RegionName=[Names[i][0] for i in range(75)]
Num=[int(Names[i][2]) for i in range(75)]



In [6]:

    
RegionName









    Out[6]:





['AME_R',
 'LO_R',
 'NO',
 'BU_R',
 'PB',
 'LH_R',
 'LAL_R',
 'SAD',
 'CAN_R',
 'AMMC_R',
 'ICL_R',
 'VES_R',
 'IB_R',
 'ATL_R',
 'CRE_R',
 'MB_PED_R',
 'MB_VL_R',
 'MB_ML_R',
 'FLA_R',
 'LOP_R',
 'EB',
 'AL_R',
 'ME_R',
 'FB',
 'SLP_R',
 'SIP_R',
 'SMP_R',
 'AVLP_R',
 'PVLP_R',
 'IVLP_R',
 'PLP_R',
 'AOTU_R',
 'GOR_R',
 'MB_CA_R',
 'SPS_R',
 'IPS_R',
 'SCL_R',
 'EPA_R',
 'GNG',
 'PRW',
 'GA_R',
 'AME_L',
 'LO_L',
 'BU_L',
 'LH_L',
 'LAL_L',
 'CAN_L',
 'AMMC_L',
 'ICL_L',
 'VES_L',
 'IB_L',
 'ATL_L',
 'CRE_L',
 'MB_PED_L',
 'MB_VL_L',
 'MB_ML_L',
 'FLA_L',
 'LOP_L',
 'AL_L',
 'ME_L',
 'SLP_L',
 'SIP_L',
 'SMP_L',
 'AVLP_L',
 'PVLP_L',
 'IVLP_L',
 'PLP_L',
 'AOTU_L',
 'GOR_L',
 'MB_CA_L',
 'SPS_L',
 'IPS_L',
 'SCL_L',
 'EPA_L',
 'GA_L']

Resize



In [7]:

    
Sm









    Out[7]:





(1024, 569, 35, 1)



In [8]:

    
S









    Out[8]:





(90, 50, 35, 1)



In [9]:

    
Masks2=np.zeros((Sm[0],Sm[1],Sm[2],87))
MasksResized=np.zeros((S[0],S[1],S[2],87))
for j in range(87):
    Masks2[:,:,:,j]=np.squeeze((Masks==j))
    for i in range(S[2]):
        MasksResized[:,:,i,j]=scipy.ndimage.zoom(Masks2[:,:,i,j], float(S[0])/Sm[0], order=0)



In [10]:

    
MasksResized.shape









    Out[10]:





(90, 50, 35, 87)

Save the risized maps



In [22]:

    
nim=nb.Nifti1Image(MasksResized[:,:,:,1:87],np.eye(4))
nb.save(nim,'/home/sophie/Desktop/Registration/862/862ResizedMaps.nii')

Now average time series in regions



In [117]:

    
TS=np.zeros((S[3],86))
for i in range(1,86):
    for j in range(S[3]):
        TS[j,i]=np.mean(np.mean(np.mean(np.multiply(MasksResized[:,:,:,i],data[:,:,:,j]))))

might need to do that in matlab parfor



In [184]:

    
pylab.rcParams['figure.figsize'] = (11, 8)
h=1.7
i=0
fig = plt.figure()
axes=plt.gca()
#axes.set_ylim([-0.1,0.4])
axes.set_xlim([-500,10000])
ax = fig.add_subplot(111)

for j in range(70,75):
    i=Num[j]
    plt.plot((TS[1000:11000,i]/np.max(TS[1000:11000,i])-h*j),color=np.squeeze(np.random.rand(3,1))/1.2)
    ax.text(10100, -h*j, RegionName[j], style='italic',fontsize=16,bbox={'alpha':0.5, 'pad':0})

plt.show









    



---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-184-ecc95bce1e1c> in <module>()
     11 for j in range(70,75):
     12     i=Num[j]
---> 13     plt.plot((TS[1000:11000,i]/np.max(TS[1000:11000,i])-h*j),color=np.squeeze(np.random.rand(3,1))/1.2)
     14     ax.text(10100, -h*j, RegionName[j], style='italic',fontsize=16,bbox={'alpha':0.5, 'pad':0})
     15 

IndexError: index 86 is out of bounds for axis 1 with size 86



In [191]:

    
TS2save=np.zeros((10000,74))
for j in range(74):
    i=Num[j]
    TS2save[:,j]=TS[1000:11000,i]/np.max(TS[1000:11000,i])



In [194]:

    
import scipy.io as sio
sio.savemat('/home/sophie/Desktop/100115RegionsTS.mat',{'TS':TS2save,'RegionName':RegionName})



In [ ]: