diff --git a/day2_exercises/fMRI_homework_20151107.py b/day2_exercises/fMRI_homework_20151107.py new file mode 100644 index 0000000..613f464 --- /dev/null +++ b/day2_exercises/fMRI_homework_20151107.py @@ -0,0 +1,88 @@ +# -*- coding: utf-8 -*- +""" +Created on Sat Nov 7 09:02:19 2015 + +@author: Elissaios +""" + +#%% Load nii image (after it was unzipped) +import nibabel as nib + +folder = r"c:\Anaconda3\Lib\site-packages\nibabel\tests\data"; +#filePath = path + 'ds107_sub001_highres.nii'; +filePath = folder + r"\example_nifti2.nii"; + +img = nib.load(filePath); +data = img.get_data(); + +#%% Anatomical image +import numpy as np +#import math; +import matplotlib.pyplot as plt + +folder = r"c:\Users\Elissaios\Box Sync\Projects\fMRI\analysis1"; +filePath = folder +r"\anatomical.txt"; +anat = np.loadtxt(filePath); + +zdim=32; +anat_slice = int(len(anat)/zdim) +for ix in range(120,200+1): + anat_mod = anat_slice % ix; + if anat_mod == 0: + ij = anat_slice/ix; + #if ij.is_integer(): + #if np.floor(ij)-ij == 0: + #if int(ij)-ij == 0: + if ij % 1 == 0: + print(ix, int(ij)) + anat3d = anat.reshape(ix,int(ij),zdim) + plt.imshow(anat3d[:,:,16],'gray') + #answer 170 by 156 + + +#%% Working with four dimensional images +import numpy as np +import nibabel as nib +import matplotlib.pyplot as plt + +folder = r"C:\Users\Elissaios\Box Sync\Projects\fMRI\analysis1"; +filePath = folder + r"\ds114_sub009_t2r1.nii"; +img = nib.load(filePath) +data = img.get_data() +volume = data[...,0] +#plt.imshow(volume[...,15]) +var_vol = np.var(data,3) +plt.figure(1) +plt.imshow(var_vol[...,14],'gray') +low_vol = volume[volume<10] +plt.figure(2) +plt.hist(low_vol) +uniq = np.unique(low_vol) + +#%% Into four dimensions +import numpy as np +import nibabel as nib +import matplotlib.pyplot as plt +folder = r"C:\Users\Elissaios\Box Sync\Projects\fMRI\analysis1"; +filePath = folder + r"\ds114_sub009_t2r1.nii"; +img = nib.load(filePath) +data = img.get_data() + + + + + + + + + + + + + + + + + + + diff --git a/day2_exercises/into4D.py b/day2_exercises/into4D.py new file mode 100644 index 0000000..4054ba3 --- /dev/null +++ b/day2_exercises/into4D.py @@ -0,0 +1,113 @@ +#%% Into four dimensions +# -*- coding: utf-8 -*- +""" +Created on Sat Nov 7 13:46:24 2015 +Exercise for Fall 2015 +@author: Elissaios +""" +# Compatibility with Python 3 +from __future__ import print_function, division + +# Standard imports of libraries +import numpy as np +import nibabel as nib +import matplotlib.pyplot as plt + +# Declarations +TR = 2.5 +std_vol = [] +ave_task = [] +ave_rest = [] +outlier_task = [] +outlier_rest = [] + +# Input data +folder = r"C:\Users\Elissaios\Box Sync\Projects\fMRI\analysis1"; +filePath = folder + r"\ds114_sub009_t2r1.nii"; +condPath = folder + r"\ds114_sub009_t2r1_cond.txt"; +img = nib.load(filePath) +data = img.get_data() +#print(img.shape) +task = np.loadtxt(condPath) +task[:,0:2] = task[:,0:2]/TR # change sec to TR for first two columns (2 is for 1+1) +dur = list(img.shape)[3] +nTR = dur +time_course = np.zeros(nTR) +nBlock = task.shape[0] + +for itr in range(0,nBlock): + durBlock = task[itr,1] + blockS = int(task[itr,0]) # block start + blockE = int(blockS+durBlock) # block end + time_course[blockS:blockE] = 1 + +plt.figure(0) +plt.ylim(-0.5,1.5) +plt.plot(time_course) + +# Create two matrices, one for each condition; it includes first 10 s of rest +is_task_tr = [time_course==1] +is_rest_tr = [time_course==0] +taskInd = np.squeeze(np.array([is_task_tr])) +restInd = np.squeeze(np.array([is_rest_tr])) +data_task = data[...,taskInd] +data_rest = data[...,restInd] + +# Identify outlier volumes in each; use 2.5 x interquartile range as cutoff +# You could do it with variance too rather than average signal (take np.std instead of np.average) +# however, for a proper analysis of outliers you probably need both (e.g. either signal +# is high in entire slice but with little variance, or signal is high is certain voxels only +# with low average value but high std) +nBlock_task = int(data_task.shape[3]) +for t in range(0,nBlock_task): + ave_task = ave_task + [np.average(data_task[...,t])] +ave_task_vol = np.average(ave_task) +iqRange_task = np.percentile(ave_task,75) - np.percentile(ave_task,25) +for t in range(0,nBlock_task): + if ave_task[t] > (ave_task_vol + (iqRange_task*2.5)): + outlier_task = outlier_task + [t] + +nBlock_rest = int(data_rest.shape[3]) +for t in range(0,nBlock_rest): + ave_rest = ave_rest + [np.average(data_rest[...,t])] +ave_rest_vol = np.average(ave_rest) +iqRange_rest = np.percentile(ave_rest,75) - np.percentile(ave_rest,25) +for t in range(0,nBlock_task): + if ave_rest[t] > (ave_rest_vol + (iqRange_rest*2.5)): + outlier_rest = outlier_rest + [t] + +cleanInd_task = np.arange(0,nBlock_task) +cleanInd_task = np.setdiff1d(np.array(cleanInd_task),np.array(outlier_task)) + +cleanInd_rest = np.arange(0,nBlock_rest) +cleanInd_rest = np.setdiff1d(np.array(cleanInd_rest),np.array(outlier_rest)) + +clean_task = data_task[...,cleanInd_task] +clean_rest = data_rest[...,cleanInd_rest] + +mean_task = np.average(clean_task,axis=3) # Average voxel volumes per task +mean_rest = np.average(clean_rest,axis=3) # Average voxel volumes per task +dTask_Rest = mean_task - mean_rest +plt.figure(1) +plt.imshow(dTask_Rest[...,10],'gray') + +mean_task = np.average(data_task,axis=3) # Average voxel volumes per task +mean_rest = np.average(data_rest,axis=3) # Average voxel volumes per task +dTask_Rest = mean_task - mean_rest +plt.figure(2) +plt.imshow(dTask_Rest[...,10],'gray') + +# Find std of volumes +dur = list(img.shape)[3] + +for t in range(0,dur): + vol = data[...,t] + #print(vol1.shape) + voxels = list(vol.shape) # find how many voxels in each dimension + slic = int(voxels[2]/2)-1 # define the middle slice + #plt.figure(t) + #plt.imshow(vol[...,slic],'gray') + std_vol = std_vol + [np.std(vol)] + +plt.figure(3) +plt.hist(std_vol) # outline in the 800s