

In [3]:

    
from skimage import io
import nibabel as nib
import os
import numpy as np
import re
import scipy.io as sio
import scipy.optimize
from Tkinter import Tk
from tkFileDialog import askopenfilename
from tkFileDialog import askdirectory
import matplotlib.pyplot as plt
import subprocess
%matplotlib inline

Import data from cluster



In [4]:

    
Tk().withdraw() # we don't want a full GUI, so keep the root window from appearing
foldername = askdirectory() # show an "Open" dialog box and return the path to the selected file
print(foldername)









    



/media/test5/FreeBehaviorPanNeuronalGCaMP6/Removed/100609series/100609series/100614ss2



In [5]:

    
path=foldername



In [6]:

    
Dataname=path.split('/')[-1]



In [7]:

    
Dataname









    Out[7]:





'100614ss2'

Open the images



In [9]:

    
tt = io.imread(path+'/'+Dataname+'-00010.tif') 
S=tt.shape
S









    Out[9]:





(40, 208, 210)



In [10]:

    
file_name_list = sorted(os.listdir(path))

Use the one below to avoid indexing problems



In [11]:

    
%%time
data=np.zeros([S[0],S[1],S[2],len(file_name_list)])
for i,image_file in enumerate(sorted(file_name_list)):
    #print(i)
    tt = io.imread(path+'/'+image_file)
    data[:,:,:,i]=tt[:][:][:]









    



CPU times: user 7min 11s, sys: 52.9 s, total: 8min 4s
Wall time: 15min 3s

%%time

file_name_list = os.listdir(path)

data=np.zeros([S[0],S[1],S[2],len(file_name_list)]) i=0 for j in range(30):

for j in range(len(file_name_list)):

for fn in os.listdir(path):

tt = io.imread(path+'/'+Dataname+'-'+str(i+1).zfill(5)+'.tif') 
if os.path.exists(path+'/'+Dataname+'-'+str(j).zfill(5)+'.tif'):
    #print path+'/'+Dataname+'-'+str(j).zfill(5)+'.tif'
    tt = io.imread(path+'/'+Dataname+'-'+str(j).zfill(5)+'.tif')
    data[:,:,:,i]=tt[:][:][:]
    i=i+1



In [12]:

    
data.shape









    Out[12]:





(40, 208, 210, 4750)



In [ ]:

    
plt.show()
plt.imshow(data[5,:,:,10],interpolation='none')
plt.show()



In [ ]:

    
%%time
#file_name_list = os.listdir(path)
data=np.zeros([S[0],S[1],S[2],len(file_name_list)])
i=0
#for j in range(30):
for j in range(len(file_name_list)):       
#for fn in os.listdir(path):
    #tt = io.imread(path+'/'+Dataname+'-'+str(i+1).zfill(5)+'.tif') 
    
    #if os.path.exists(path+'/'+Dataname+'-'+str(j).zfill(5)+'.tif'):
    if os.path.exists(path+'/'+Dataname+'-'+str(j)+'.tif'):
        #print path+'/'+Dataname+'-'+str(j).zfill(5)+'.tif'
        tt = io.imread(path+'/'+Dataname+'-'+str(j)+'.tif')
        #tt = io.imread(path+'/'+Dataname+'-'+str(j).zfill(5)+'.tif')
        data[:,:,:,i]=tt[:][:][:]
        i=i+1



In [ ]:

    
data.shape

Find end of onset of light and begining of offset (to align to behavior)

Calculate average time series



In [13]:

    
M=np.mean(np.mean(np.mean(data,0),0),0)
Mav=M.mean()
plt.title('Mean_intensity_profile')
plt.plot(M,'+')









    Out[13]:





[<matplotlib.lines.Line2D at 0x7f98691e3b50>]

Get approximate on and off times



In [ ]:

    
liston=[i for i in range(len(M)) if M[i]>Mav*0.8]
liston[0]

Model for fitting onset and offset



In [ ]:

    
def model(x,a,b,c,d):
    if x<a:
        return b
    elif x<c:
        return b+(x-a)*d
    else:
        return (c-a)*d+b

Model onset and find precise onset time



In [ ]:

    
Ms=M[range(liston[0]-8,liston[0]+8)]



In [ ]:

    
def Sq(X):
    return sum([(model(i,X[0],X[1],X[2],X[3])-Ms[i])**2 for i in range(len(Ms))])



In [ ]:

    
res = scipy.optimize.minimize(Sq,x0=[7,0.3,9,0.7])
ON=liston[0]-8+res.x[2]
print(ON)



In [14]:

    
#ONint=np.int(np.ceil(ON))
ONint=1
print(ONint)
plt.plot(np.squeeze(M[range(liston[0]-8,liston[0]+8)]),'+')
plt.plot(np.arange(0,len(Ms),0.1),[model(i,res.x[0],res.x[1],res.x[2],res.x[3]) for i in np.arange(0,len(Ms),0.1)])
plt.title('Fit_light_on')
plt.show()









    



1






    



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-14-0a0741cce6f5> in <module>()
      2 ONint=1
      3 print(ONint)
----> 4 plt.plot(np.squeeze(M[range(liston[0]-8,liston[0]+8)]),'+')
      5 plt.plot(np.arange(0,len(Ms),0.1),[model(i,res.x[0],res.x[1],res.x[2],res.x[3]) for i in np.arange(0,len(Ms),0.1)])
      6 plt.title('Fit_light_on')

NameError: name 'liston' is not defined

Model offset and find precise offset time



In [ ]:

    
Ms=M[range(liston[len(liston)-1]-6,liston[len(liston)-1]+6)]



In [ ]:

    
def Sq(X):
    return sum([(model(i,X[0],X[1],X[2],X[3])-Ms[i])**2 for i in range(len(Ms))])



In [ ]:

    
res = scipy.optimize.minimize(Sq,x0=[6,3,8,-1])



In [ ]:

    
OFF=liston[len(liston)-1]-6+res.x[0]
#OFF=liston[len(liston)-1]
print(OFF)
OFFint=np.int(np.floor(OFF))
print(OFFint)



In [ ]:

    
plt.plot(np.squeeze(Ms),'+')
plt.plot(np.arange(0,len(Ms),0.1),[model(i,res.x[0],res.x[1],res.x[2],res.x[3]) for i in np.arange(0,len(Ms),0.1)])
plt.title('Fit_light_off')
plt.show()



In [ ]:

    
data.shape

Open images times



In [ ]:

    
for f in os.listdir('/home/sophie/Downloads/csvs'):
    if Dataname[:6] in f:
        if f.endswith('csv'):
            TimeFile='/home/sophie/Downloads/csvs/Data'+Dataname[:6]+'_.csv'
            Listfile = open(TimeFile, 'r')
            ListTime = [line.split('\n')[0] for line in Listfile.readlines()]
            Timespl=[float(ListTime[i].split(',')[2]) for i in range(1,len(ListTime))]
        elif "Original" or "Info" in f:
            TimeFile='/home/sophie/Downloads/csvs/'+f
            with open(TimeFile, 'r') as metafile:
                lines = metafile.readlines()
            time_from_start_list = []
            for line in lines:
                if "Time_From_Start" in line:
                    split_line = line.replace('\t', ' ').replace(' = ', ' ').strip().split(' ')
                    time_from_start_list.append((int(split_line[1]),float(split_line[3])))
            Timespl = list(zip(*sorted(time_from_start_list))[1])



In [ ]:

    
S

Get times corresponding to images during light on (excitation light completely on : t=0)



In [ ]:

    
print(ONint)
print(OFFint)
#print(ON)



In [ ]:

    
(11974-498)*0.05



In [ ]:

    
len(Timespl)



In [ ]:

    
TimeOn=[Timespl[i] for i in range(ONint,(OFFint+1))]
Tinit=(ON-(ONint-1))*(Timespl[ONint]-Timespl[ONint-1])+Timespl[ONint-1]
Toff=(OFFint+1-OFF)*(Timespl[OFFint+1]-Timespl[OFFint])+Timespl[OFFint]
Toff-Tinit



In [ ]:

    
Timespl[ONint]-Timespl[ONint-1]



In [ ]:

    
TimeOn[0]-Tinit



In [ ]:

    
Dataname[:6]



In [ ]:

    
TimeOnFinal=np.array(TimeOn)-Tinit



In [ ]:

    
sio.savemat('/home/sophie/Desktop/'+Dataname+'TimeFluoOn.mat', {'TimeFluoOn':TimeOnFinal})



In [ ]:

    
TotalTimeOn=Toff-Tinit



In [ ]:

    
sio.savemat('/home/sophie/Desktop/'+Dataname+'TotalTimeOn.mat', {'TotalTimeOn':TotalTimeOn})

Keep only the frames for which the excitation is on and save



In [ ]:

    
OFFint=4750



In [ ]:

    
D4=np.transpose(data[:,:,:,range(ONint,(OFFint))],(2,1,0,3))
#D4=np.transpose(data[:,:,:,range(ONint,40000)],(2,1,0,3))
nim=nib.Nifti1Image(D4,np.eye(4))
nib.save(nim,'/media/test15/'+Dataname+'on.nii.gz')



In [ ]:

    
del data
del D4

Movement Correction



In [ ]:

    
bashCommand = "3dvolreg -prefix "+"/media/sophie/Elements/"+Dataname+"onreg.nii "+"/mediasophie/Elements/"+Dataname+"on.nii.gz"
import subprocess
process = subprocess.Popen(bashCommand.split(), stdout=subprocess.PIPE)
output, error = process.communicate()



In [ ]:

    
import matlab.engine



In [ ]:

    
eng=matlab.engine.start_matlab()



In [ ]:

    
bashCommand = "matlab -r 'RoughAllSteps.m'"
import subprocess
process = subprocess.Popen(bashCommand.split(), stdout=subprocess.PIPE)
output, error = process.communicate()



In [ ]:

    
%matlab
RoughAllSteps



In [ ]: