%pylab inline

import os # Module with commands to interact with the operating system

import numpy as np # this is the Python package for handling arrays

import matplotlib.pyplot as plt # This is for making 2D plots like MATLAB

anat_fname = 'ds107_sub001_highres.nii.gz'
anat_fname

import nibabel as nib

anat_img = nib.load(anat_fname)
anat_img

anat_img.shape

anat_arr = anat_img.get_data() # get the data array
type(anat_arr)

anat_arr.shape

middle_slice = anat_arr[:, :, 96]
middle_slice.shape

plt.imshow(middle_slice)

middle_slice

middle_slice[122:133, 122:133]

anat_arr.max()

anat_arr.min()

anat_arr.mean()

anat_arr_1d = anat_arr.ravel() # as 1D vector
tmp = plt.hist(anat_arr_1d, 100) # 100 histogram bins for nice detail

percentile_95 = np.percentile(anat_arr_1d, 95)
clipped_vector = anat_arr_1d[anat_arr_1d <= percentile_95]
tmp = plt.hist(clipped_vector, 100)

hdr = anat_img.get_header()
print hdr

aff = anat_img.get_affine()
aff

array_coords = [0, 0, 0]
homogenous_array_coords = [0, 0, 0, 1]
mm_coordinates = aff.dot(homogenous_array_coords)
mm_coordinates

homogenous_array_coords = [127, 127, 46, 1] # note the extra 1 again for homogenous coordinates
mm_coordinates = aff.dot(homogenous_array_coords)
mm_coordinates

bold_fname = 'ds107_sub001_run01.nii.gz'
bold_fname

bold_img = nib.load(bold_fname)
bold_img

bold_img.shape

bold_arr = bold_img.get_data()
bold_vol0_arr = bold_arr[:, :, :, 0]
bold_vol0_arr.shape

bold_middle_slice = bold_vol0_arr[:, :, 17]
plt.imshow(bold_middle_slice)

time_course_mean = bold_arr.mean(axis=3) # the last axis (0-based numbering)
time_course_mean.shape


plt.imshow(time_course_mean[:, :, 17])

time_course_std = bold_arr.std(axis=3) # the last axis (0-based numbering)
middle_slice_std = time_course_std[:, :, 17]
plt.imshow(middle_slice_std)

max_index = np.argmax(middle_slice_std) # flattens the array first
max_coord_2d = np.unravel_index(max_index, middle_slice_std.shape) # get 2D coordinate equivalent
max_coord_2d

max_coord_3d = max_coord_2d + (17,) # add the slice number as third coordinate
max_coord_3d

time_course_std[max_coord_3d]

middle_slice_std.max()

max_std_over_time = bold_arr[max_coord_3d]
plt.plot(max_std_over_time)

import nipy.algorithms.diagnostics as nads

tdiff = nads.time_slice_diffs(bold_arr, time_axis=-1, slice_axis=-2)
nads.plot_tsdiffs(tdiff)