%pylab inline
import nibabel as nib

import nibabel

import numpy as np
import nibabel as nib
import nibabel.nicom as nic
import scipy as sc
import os
from os.path import join as pjoin

os.path.realpath(nib.__file__)

pwd

CURDIR = os.path.abspath(os.path.curdir)
ddata = os.path.join(CURDIR,"bolddata")
print os.listdir(ddata)
print os.listdir(CURDIR)

%alias
# install mricron ... 
# change to directory with epi_01 : os.chdir( pjoin(CURDIR,bolddata) )
#!dcm2nii epi_01
# move back : os.chdir( CURDIR )

ls


# !mv epi_01/20090421_155120EPI9613TRSAGs005a001.nii.gz dimg.nii.gz

fimg = pjoin(ddata, "bold.nii.gz")
print fimg
img = nib.load(fimg)

print img.shape
arr = img.get_data()
print arr.shape

nib.nifti1.Nifti1Image?

small_arr = arr[:,:,:,:42]
print small_arr.shape

aff = img.get_affine()
print aff
#
fname_small_img = pjoin(ddata, "smallbold.nii.gz")
small_img = nib.nifti1.Nifti1Image(small_arr, aff)
nib.save(small_img, fname_small_img)

plt.gray()
print fname_small_img
simg = nib.load(fname_small_img)
imshow(simg.get_data()[:,:,17,5])

plot(simg.get_data()[24,12,17,:])

arr = simg.get_data()
n_scans = arr.shape[3]
#
temporal_mean = arr.mean(axis=3)
print temporal_mean.shape
print temporal_mean.mean()
print arr.mean()
#
imshow(simg.get_data()[:,:,17,5] - temporal_mean[:,:,17])

print "temporal_mean[:,:,1].shape : ", temporal_mean[:,:,1].shape

t = 12
arr12 = temporal_mean - arr[:,:,:,t]
print arr12.shape
arr12sq = arr12**2
print arr12sq.shape
print np.sqrt(arr12sq.mean())

rms = np.zeros(n_scans) # initialize an array to hold our rms values
for t in range(n_scans):
    squared_difference = (temporal_mean - arr[:,:,:,t])**2
    rms[t] = np.sqrt(np.mean(squared_difference))
plt.plot(rms)

ashape = np.asarray(arr.shape)

print ashape, 64*64*35, ashape[:3].prod()
print ashape[:3].prod(),ashape[3]


# another way ?
ashape = np.asarray(arr.shape)
arr_reshaped = arr.reshape(ashape[:3].prod(),ashape[3])
print arr_reshaped.shape
#plot(arr_reshaped.std(axis=0))

arr_demeaned = arr - arr.mean(axis=3)


tmp = np.arange(24)

tmp

tmp.shape = 2, 3*4
print tmp
print tmp.ravel()

# reshape 
print tmp.reshape(3*4, 2)

tmp.shape = 2, 3*4
print tmp

# transpose 
print tmp.T
print tmp.T.ravel()

print tmp

# has our array changed ? 
tmp.shape = 2, 3, 4
print tmp
print tmp.ravel()

print tmp[0].shape, "\n",  tmp[0],  "\n",  tmp[-1] 
print tmp[:,:,0].shape, tmp[:,:,-1].shape, "\n", tmp[:,:,0] 


tmp.shape = 2, 3, 4
print tmp

tmp = tmp.reshape(4, 3, 2)
# could also be done with tmp.shape = 4, 2, 3 
print tmp
print tmp.ravel()


tmp = tmp.reshape(2, 3, 4)
# could also be done with tmp.shape = 4, 2, 3 
print tmp
print tmp.ravel()

# what does transpose do ? 
tmp.shape = 2, 3, 4
print tmp

print tmp.T.shape
print tmp.T
print tmp.T.ravel()

print tmp.shape, tmp.ravel()

print np.rollaxis(tmp, 2, 0).shape, np.rollaxis(tmp, 2, 0).ravel()
print np.rollaxis(np.rollaxis(tmp, 2, 0), 2, 1).shape, \
      np.rollaxis(np.rollaxis(tmp, 2, 0), 2, 1).ravel()

#if more than 3 dim
tmp4d_shape = (2,3,4,5)
tmp4d = np.arange(np.prod(tmp4d_shape)).reshape(tmp4d_shape)
print tmp4d[0,0:2,:,:]
print np.prod(tmp4d_shape)
print tmp4d.T.shape
print tmp4d.T[0,0:2,:,:]
print tmp4d.T.ravel()[:21]

mean_img = arr.mean(axis=3)
mean_img_shape = np.asarray(mean_img.shape)
print mean_img.shape, mean_img_shape

tximg = arr.reshape(mean_img_shape.prod(),arr.shape[3]).T
print tximg.shape

# last dimension is the same: broadcasting
tximg_centred = tximg - mean_img.ravel()

print tximg_centred.shape
print mean_img.ravel().shape

#fig = figure()
ax = subplots(1,1)

ax[1].plot(tximg_centred.std(axis=1))
ax[1].set_xlabel('TR')


tximg = tximg_centred + mean_img.ravel()
arr2 = tximg.T.reshape(arr.shape)
print (arr - arr2).max(), (arr - arr2).min()
print np.allclose(arr,arr2)
print np.all((arr-arr2) == 0)
print np.any(arr-arr2)


arr_std = arr.std(axis = 3)
(fig,axes) = subplots(3,3)

for i in range(3):
    for j in range(3):
        axes[i,j].imshow(arr_std[:,:,(3*i + j)*4])

diff_arr = arr[:,:,:,:-1] - arr[:,:,:,1:] 

print arr.shape, diff_arr.shape

imagedim = np.prod(arr.shape[:3])
print imagedim

squarediff = (diff_arr.reshape(imagedim, arr.shape[3]-1))**2
print squarediff.shape
m2diff = squarediff.mean(axis=0)

figure(num=None, figsize=(12,2))
plot(m2diff)

print np.argmax(m2diff)

nuisance = np.zeros(n_scans)
nuisance[np.argmax(m2diff)] = 1
print nuisance

# get the variance per slice - no difference between adjacent slices
ash = arr.shape
print ash

# slice std : mean taken within slice
slice_std = arr.reshape((ash[0]*ash[1], ash[2], ash[3])).std(axis=0)
slice_mean = arr.reshape((ash[0]*ash[1], ash[2], ash[3])).mean(axis=0)
print slice_std.shape
#(fig, axes) = subplots(1,2)
#axes[0].plot(slice_std.T)
#axes[1].plot(slice_mean.T)


# display slice std across TR: mean taken across slice
m = arr.mean(axis=3)
m_shape = np.asarray(m.shape)
tximg = arr.reshape(m_shape[0:3].prod(),-1)
tximg_T = np.rollaxis(tximg, 0, 2)
print np.rollaxis(tximg, 0, 2).shape

arr_demean = tximg_T - tximg_T.mean(axis=0)
arr_demean_reshaped = arr_demean.reshape((ash[-1],)+ash[:3]) 
print (ash[-1],)+ash[:3]
print 'arr_demean_reshaped.shape',  arr_demean_reshaped.shape

print arr_demean_reshaped[:,32,32,17].mean()

arr_demean_reshaped_std = arr_demean_reshaped.std(axis=0)

plot(arr_demean_reshaped_std.reshape(64*64, ash[2]).mean(axis=0))

# do it for each time ?

# do it for the difference between adjacent scans ?
# print diff_arr.shape


import nipy.algorithms.diagnostics as nads

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

import scipy.linalg as lin

print ash
tximg = arr.reshape(np.prod(ash[:3]),ash[3]).T
print tximg.shape
aTa = tximg.dot(tximg.T)
u,s,vh = lin.svd(aTa)


plot(s)

plot(u[:,0])

tximg_centred = tximg - tximg.mean(axis=0)
print tximg_centred.shape
aTa = tximg_centred.dot(tximg_centred.T)
v,s,vh = lin.svd(aTa)


plot(s)

u1 = tximg_centred.T.dot(v[:,0])

print u1.shape
print tximg_centred.shape
u1 = u1.reshape(ash[:3])
print u1.shape
imshow(u1[:,:,17])

# load the image again - to be sure
fimg = pjoin(ddata, "smallbold.nii.gz")
print fimg
img = nib.load(fimg)
arr = img.get_data()

# increase the mean at volume 18
arr[:,:,:,18] = arr[:,:,:,18] + 30

def std_diff_ts(a4d, timeaxis = 'last'):
    """
    Takes a 4d numpy array, take the diff between adjacent times, 
    compute and returns the std of the diff volume time series

    Inputs
    ------
    arr4d: numpy array
        The array on which we compute the std of the diff time series
    timeaxis: integer or in ['last', 'first']
        The axis where the time is

    Outputs
    -------
    (Ntmes,) numpy array
        The array with the std of the diff volume time series
    """
#    pass
    ash = np.asarray(a4d.shape)
    if len(a4d.shape) != 4:
        raise ValueError(" Array is not of shape 4")
    
    if timeaxis == 'first':
            timeaxis = 0
    if timeaxis != 'last':
        np.rollaxis(a4d, timeaxis, -1)
        
    print ash, ash[:3].prod(),ash[3]
    diff_arr = (a4d[:,:,:,:-1] - a4d[:,:,:,1:])**2
    diff_arr = diff_arr.reshape(ash[:3].prod(),ash[3]-1)
    return diff_arr.mean(axis=0)

len(arr.shape)

stdts = std_diff_ts(arr)
plot(stdts)

print np.argmax(stdts)
rk = np.argsort(stdts)
print argmax([stdts[rk[-1]-1],stdts[rk[-1]+1]])

bad_volume = np.argmax(stdts) + argmax([stdts[rk[-1]-1],stdts[rk[-1]+1]])