In [1]:
"""
Function to compute Avalanche images
Expects full file path to resting state 4D nifti image
Currently, output images will be written to wd
"""
from __future__ import absolute_import, division, print_function
import numpy as np
import nibabel as nib
import os as os
from scipy.ndimage.measurements import label
from scipy.stats import zscore
from nilearn.masking import compute_epi_mask
from nilearn.masking import apply_mask
from nilearn.masking import unmask
from scipy.ndimage.morphology import generate_binary_structure
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline
Load EPI data
In [3]:
func_file = '/Users/dlurie/Dropbox/Projects/avalanche/avalanche/data/ds000133/sub-01/post/sub-01_ses-post_task-rest_run-01_bold_space-MNI152NLin2009cAsym_clean_std_fwhm6.nii.gz'
Load mask
In [4]:
mask_file = '/Users/dlurie/Dropbox/Projects/avalanche/avalanche/data/ds000133/sub-01/post/sub-01_ses-post_task-rest_run-01_bold_space-MNI152NLin2009cAsym_brainmask.nii.gz'
In [5]:
def avalanche(func_file, mask_file, thresh, connectivity, out_dir):
"""
Load and preprocess data.
"""
# Load the functional image.
func_img = nib.load(func_file)
func_data = func_img.get_data()
if mask_file is not None:
# Load the mask image.
mask_img = nib.load(mask_file)
else:
# Create a new brain mask from the EPI data.
mask_img = compute_epi_mask(func_img)
# Mask the functional data
func_data_masked = apply_mask(func_img, mask_img)
# Z-score the masked functional data.
func_data_masked_z = zscore(func_data_masked, axis=0)
"""
Create the initial point-process image and 4D cluster image.
"""
# Apply the threshold and get the point-process data.
pp_data_masked = func_data_masked_z >= thresh
# Unmask the point-process data, creating an image, and save it.
pp_img = unmask(pp_data_masked, mask_img)
pp_img.to_filename(os.path.join(out_dir,'{0}_thresh-{1}_conn-{2}_ppimage.nii.gz'.format(os.path.basename(func_file)[:-7], str(thresh), str(connectivity))))
# Label 3D and 4D clusters.
struct = generate_binary_structure(4,connectivity)
pp_data = pp_img.get_data()
cluster_data, num_clusters = label(pp_data, struct)
# Create a cluster image and save it.
cluster_img = nib.Nifti1Image(cluster_data, pp_img.affine)
cluster_img.to_filename(os.path.join(out_dir,'{0}_thresh-{1}_conn-{2}_clusterimage.nii.gz'.format(os.path.basename(func_file)[:-7], str(thresh), str(connectivity))))
"""
Remove clusters that only exist in a single volume.
"""
# Get the dimensions of the functional image.
(n_x, n_y, n_z, n_t) = cluster_img.shape
# Mask the cluster image to get a 2D array.
cluster_data_masked = apply_mask(cluster_img, mask_img)
# Initialize a matrix to count unique clusters in each volume.
uMat = np.zeros(shape=(num_clusters, n_t))
# loop through each time point and identify cluster members:
for i in range(n_t):
#determine uniques at this time point:
inds = np.unique(cluster_data_masked[i,:])
#Note: 0 will always be a unique value, and we don't
#actually care about zeros, so ignore them:
uMat[(inds[1:].astype('int')-1),i] = 1
#now the uMat matrix is populated with indicators for each of the "avalanches"
#but we need to determine which are truly avalanches (i.e., occur in more than
#one time point), and which are not
avarray= np.zeros((uMat.shape[0], n_t+1))
#fill this with binary data:
avarray[:,1:]= uMat
#and take the difference:
avarray[:,1:] = np.diff(avarray, axis=1)
#and then take it again to identify 1 to -1 (diff==-2) points - these are not
#avalanches!!!
tmparray = np.diff(avarray, axis=1)
#now we can efficiently remove all fake avalanches:
fakes = np.asarray(np.where(np.sum(tmparray==-2,axis=1)))+1
cluster_data[np.isin(cluster_data,fakes)]=0
#now single time point clusters should be removed
#and we want to regenerate our labels:
cluster_data, num_clusters = label(cluster_data, struct)
# Create an image from the new cluster data (now avalanche data) and then save it.
cluster_img = nib.Nifti1Image(cluster_data, header=func_img.header, affine=func_img.affine)
cluster_img.to_filename(os.path.join(out_dir,'{0}_thresh-{1}_conn-{2}_avalancheimage.nii.gz'.format(os.path.basename(func_file)[:-7], str(thresh), str(connectivity))))
""""
img_new = nb.Nifti1Image(labeled_array, header=img.header, affine=img.affine)
# Reconstruct the 4D volume
ava_file = os.path.join(os.getcwd(), 'avalancheLabels.nii.gz')
img_new.to_filename(ava_file)
#make a new 4D array that is 1 time point greater than original data:
transition_array = np.zeros((n_x, n_y, n_z, n_t+1))
#the point process should come from the de-faked avalance data:
signal = labeled_array
signal[signal>0]=1
#fill this transition array with binary data and leave first time point (t[0])
#empty:
transition_array[:,:,:,1:] = signal
#take the difference between points across the 4th (time) dimension:
onset_array = np.diff(transition_array, axis=3)
onset_array[onset_array==-1] = 0
#the onset array is back in the same time order and number of time points
#as the original arrays...
img_new = nb.Nifti1Image(onset_array, header=img.header, affine=img.affine)
# Reconstruct the 4D volume
pp_file = os.path.join(os.getcwd(), 'binary_pp.nii.gz')
img_new.to_filename(pp_file)
Nvoxels = np.unique(labeled_array, return_counts='true')
Nvoxels = np.asarray(Nvoxels).T
#get rid of the 0 count (these aren't avalanches)
Nvoxels = Nvoxels[1:]
#could return num_features, Nvoxels
"""
In [6]:
for thresh in [1, 1.5, 2, 2.5]:
for struct in [1, 2, 3, 4]:
avalanche(func_file, mask_file, thresh, struct, '/Users/dlurie/Dropbox/Projects/avalanche/avalanche/data/ds000133/results/')