In [ ]:

    

import nibabel, nipy, nilearn
from nipy.labs.mask import compute_mask
from nilearn.masking import compute_epi_mask
from nilearn.image import mean_img, new_img_like
import pylab as plt
import seaborn as sns
from nilearn.plotting import plot_epi, plot_anat, plot_roi
import numpy as np
%matplotlib inline


    

In [ ]:

    

compute_mask?


    

In [ ]:

    

in_file = "data/sub-ben01_task-unknown_bold.nii.gz"


    

In [ ]:

    

in_4d_nii = nibabel.load(in_file)


    

In [ ]:

    

mean_nii = mean_img(in_4d_nii)
plot_anat(mean_nii)


    

In [ ]:

    

#mask_nii = compute_epi_mask(in_4d_nii)

mask_data = compute_mask(mean_nii.get_data())
mask_nii = new_img_like(mean_nii, mask_data)

plot_roi(mask_nii, mean_nii)


    

In [ ]:

    

a = np.where(mask_nii.get_data() != 0)
bbox = np.max(a[0]) - np.min(a[0]), np.max(a[1]) - np.min(a[1]), np.max(a[2]) - np.min(a[2])
print(bbox)
print(np.argmax(bbox))
longest_axis = np.argmax(bbox)


    

In [ ]:

    

from scipy import ndimage

# Input here is a binarized and intersected mask data from previous section
dil_mask = ndimage.binary_dilation(mask_nii.get_data(), iterations=int(mask_nii.shape[longest_axis]/8))

# Now, we visualize the same using `plot_roi` with data being converted to Nifti
# image. In all new image like, reference image is the same but second argument
# varies with data specific
dil_mask_nii = new_img_like(
    mask_nii,
    dil_mask.astype(np.int))
plot_roi(dil_mask_nii, mean_nii)


    

In [ ]:

    

dil_mask.shape


    

In [ ]:

    

rep = list(mask_nii.shape)
rep[longest_axis] = -1
new_mask_2d = dil_mask.max(axis=longest_axis).reshape(rep)
new_mask_2d.shape


    

In [ ]:

    

rep = [1,1,1]
rep[longest_axis] = mask_nii.shape[longest_axis]
new_mask_3d = np.logical_not(np.tile(new_mask_2d, rep))
#new_mask_3d =


    

In [ ]:

    

new_mask_nii = new_img_like(
    mask_nii,
    new_mask_3d.astype(np.int))
plot_roi(new_mask_nii, mean_nii)


    

In [ ]:

    

from scipy.stats.mstats import zscore

data4d = in_4d_nii.get_data()
for slice_i in range(in_4d_nii.shape[2]):
    slice_data = data4d[:,:,slice_i,:][new_mask_3d[:,:,slice_i]].mean(axis=0)
    slice_data = zscore(slice_data)
    plt.plot(slice_data)


    

In [ ]:

    

data4d[:,:,slice_i,:][new_mask_3d[:,:,slice_i]].shape


    

In [ ]:

    

data4d.shape


    

In [ ]:

    

compute_mask?


    

In [ ]:

    

def plot_spikes(in_file, skip=0):
    in_4d_nii = nibabel.load(in_file)
    mean_nii = mean_img(in_4d_nii)
    plot_anat(mean_nii)
    mask_data = compute_mask(mean_nii.get_data())
    mask_nii = new_img_like(mean_nii, mask_data)

    plot_roi(mask_nii, mean_nii)
    
    a = np.where(mask_nii.get_data() != 0)
    bbox = np.max(a[0]) - np.min(a[0]), np.max(a[1]) - np.min(a[1]), np.max(a[2]) - np.min(a[2])
    print(bbox)
    print(np.argmax(bbox))
    longest_axis = np.argmax(bbox)
    
    from scipy import ndimage

    # Input here is a binarized and intersected mask data from previous section
    dil_mask = ndimage.binary_dilation(mask_nii.get_data(), iterations=int(mask_nii.shape[longest_axis]/8))

    # Now, we visualize the same using `plot_roi` with data being converted to Nifti
    # image. In all new image like, reference image is the same but second argument
    # varies with data specific
    dil_mask_nii = new_img_like(
        mask_nii,
        dil_mask.astype(np.int))
    plot_roi(dil_mask_nii, mean_nii)
    
    rep = list(mask_nii.shape)
    rep[longest_axis] = -1
    new_mask_2d = dil_mask.max(axis=longest_axis).reshape(rep)
    new_mask_2d.shape
    
    rep = [1,1,1]
    rep[longest_axis] = mask_nii.shape[longest_axis]
    new_mask_3d = np.logical_not(np.tile(new_mask_2d, rep))
    
    new_mask_nii = new_img_like(
    mask_nii,
    new_mask_3d.astype(np.int))
    plot_roi(new_mask_nii, mean_nii)
    
    from scipy.stats.mstats import zscore

    data4d = in_4d_nii.get_data()[:,:,:,skip:]
    plt.figure()
    for slice_i in range(in_4d_nii.shape[2]):
        slice_data = data4d[:,:,slice_i,:][new_mask_3d[:,:,slice_i]].mean(axis=0)
        #slice_data = zscore(slice_data)
        plt.plot(slice_data)
    plt.title(in_file)


    

In [ ]:

    

plot_spikes("D:/example_artifacts_dataset/sub-ben01/func/sub-ben01_task-unknown_bold.nii.gz")


    

In [ ]:

    

from glob import glob


    

In [ ]:

    

for file in glob("D:/*/sub-*/func/sub-*_task-*_bold.nii.gz"):
    plot_spikes(file)


    

In [ ]:

    

plot_anat?


    

In [ ]:

    

glob("D:/*/sub-*/func/sub-*_task-*_bold.nii.gz")