%pylab inline
import nibabel as nib
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline]. For more information, type 'help(pylab)'.
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__)
-c:3: UserWarning: The DICOM readers are highly experimental, unstable, and only work for Siemens time-series at the moment Please use with caution. We would be grateful for your help in improving them
'/home/jb/.local/lib/python2.7/site-packages/nibabel/__init__.pyc'
pwd
u'/home/jb/data/qcpna'
CURDIR = os.path.abspath(os.path.curdir)
ddata = os.path.join(CURDIR,"bolddata")
print os.listdir(ddata)
print os.listdir(CURDIR)
['smallbold.nii.gz', 'epi_01', 'dimg.nii.gz', 'bold.nii.gz'] ['bold_QC_teacher.ipynb', 'bolddata', 'bold_QC.ipynb']
%alias
# install mricron ...
# change to directory with epi_01 : os.chdir( pjoin(CURDIR,bolddata) )
#!dcm2nii epi_01
# move back : os.chdir( CURDIR )
Total number of aliases: 12
[('cat', 'cat'), ('cp', 'cp -i'), ('ldir', 'ls -F -o --color %l | grep /$'), ('lf', 'ls -F -o --color %l | grep ^-'), ('lk', 'ls -F -o --color %l | grep ^l'), ('ll', 'ls -F -o --color'), ('ls', 'ls -F --color'), ('lx', 'ls -F -o --color %l | grep ^-..x'), ('mkdir', 'mkdir'), ('mv', 'mv -i'), ('rm', 'rm -i'), ('rmdir', 'rmdir')]
ls
bolddata/ bold_QC.ipynb bold_QC_teacher.ipynb
# !mv epi_01/20090421_155120EPI9613TRSAGs005a001.nii.gz dimg.nii.gz
fimg = pjoin(ddata, "bold.nii.gz")
print fimg
img = nib.load(fimg)
/home/jb/data/qcpna/bolddata/bold.nii.gz
print img.shape
arr = img.get_data()
print arr.shape
(64, 64, 35, 165) (64, 64, 35, 165)
Let's reduce the size of this time series : take the 42 first scans
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)
(64, 64, 35, 42) [[ 3. 0. 0. -93. ] [ 0. 3. 0. -103.41815186] [ 0. 0. 3. -45.47966003] [ 0. 0. 0. 1. ]]
plt.gray()
print fname_small_img
simg = nib.load(fname_small_img)
imshow(simg.get_data()[:,:,17,5])
/home/jb/data/qcpna/bolddata/smallbold.nii.gz
<matplotlib.image.AxesImage at 0x4475fd0>
plot(simg.get_data()[24,12,17,:])
[<matplotlib.lines.Line2D at 0x46f6090>]
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
(64, 64, 35) 273.428762921 273.428762921 temporal_mean[:,:,1].shape : (64, 64)
t = 12
arr12 = temporal_mean - arr[:,:,:,t]
print arr12.shape
arr12sq = arr12**2
print arr12sq.shape
print np.sqrt(arr12sq.mean())
(64, 64, 35) (64, 64, 35) 15.4273905419
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)
[<matplotlib.lines.Line2D at 0x47293d0>]
ashape = np.asarray(arr.shape)
print ashape, 64*64*35, ashape[:3].prod()
print ashape[:3].prod(),ashape[3]
[64 64 35 42] 143360 143360 143360 42
# 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))
(143360, 42)
arr_demeaned = arr - arr.mean(axis=3)
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-21-bb07f6068054> in <module>() ----> 1 arr_demeaned = arr - arr.mean(axis=3) ValueError: operands could not be broadcast together with shapes (64,64,35,42) (64,64,35)
So - that didn't work. Any idea why ? - how to correct ?
tmp = np.arange(24)
tmp
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23])
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
[[ 0 1 2 3 4 5 6 7 8 9 10 11] [12 13 14 15 16 17 18 19 20 21 22 23]] [ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23] [[ 0 1] [ 2 3] [ 4 5] [ 6 7] [ 8 9] [10 11] [12 13] [14 15] [16 17] [18 19] [20 21] [22 23]] [[ 0 1 2 3 4 5 6 7 8 9 10 11] [12 13 14 15 16 17 18 19 20 21 22 23]] [[ 0 12] [ 1 13] [ 2 14] [ 3 15] [ 4 16] [ 5 17] [ 6 18] [ 7 19] [ 8 20] [ 9 21] [10 22] [11 23]] [ 0 12 1 13 2 14 3 15 4 16 5 17 6 18 7 19 8 20 9 21 10 22 11 23] [[ 0 1 2 3 4 5 6 7 8 9 10 11] [12 13 14 15 16 17 18 19 20 21 22 23]]
# has our array changed ?
tmp.shape = 2, 3, 4
print tmp
print tmp.ravel()
[[[ 0 1 2 3] [ 4 5 6 7] [ 8 9 10 11]] [[12 13 14 15] [16 17 18 19] [20 21 22 23]]] [ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23]
print tmp[0].shape, "\n", tmp[0], "\n", tmp[-1]
print tmp[:,:,0].shape, tmp[:,:,-1].shape, "\n", tmp[:,:,0]
(3, 4) [[ 0 1 2 3] [ 4 5 6 7] [ 8 9 10 11]] [[12 13 14 15] [16 17 18 19] [20 21 22 23]] (2, 3) (2, 3) [[ 0 4 8] [12 16 20]]
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()
[[[ 0 1 2 3] [ 4 5 6 7] [ 8 9 10 11]] [[12 13 14 15] [16 17 18 19] [20 21 22 23]]] [[[ 0 1] [ 2 3] [ 4 5]] [[ 6 7] [ 8 9] [10 11]] [[12 13] [14 15] [16 17]] [[18 19] [20 21] [22 23]]] [ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23] [[[ 0 1 2 3] [ 4 5 6 7] [ 8 9 10 11]] [[12 13 14 15] [16 17 18 19] [20 21 22 23]]] [ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23]
# what does transpose do ?
tmp.shape = 2, 3, 4
print tmp
print tmp.T.shape
print tmp.T
print tmp.T.ravel()
[[[ 0 1 2 3] [ 4 5 6 7] [ 8 9 10 11]] [[12 13 14 15] [16 17 18 19] [20 21 22 23]]] (4, 3, 2) [[[ 0 12] [ 4 16] [ 8 20]] [[ 1 13] [ 5 17] [ 9 21]] [[ 2 14] [ 6 18] [10 22]] [[ 3 15] [ 7 19] [11 23]]] [ 0 12 4 16 8 20 1 13 5 17 9 21 2 14 6 18 10 22 3 15 7 19 11 23]
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()
(2, 3, 4) [ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23] (4, 2, 3) [ 0 4 8 12 16 20 1 5 9 13 17 21 2 6 10 14 18 22 3 7 11 15 19 23] (4, 2, 3) [ 0 4 8 12 16 20 1 5 9 13 17 21 2 6 10 14 18 22 3 7 11 15 19 23] (4, 3, 2) [ 0 12 4 16 8 20 1 13 5 17 9 21 2 14 6 18 10 22 3 15 7 19 11 23]
#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]
[[[ 0 1 2 3 4] [ 5 6 7 8 9] [10 11 12 13 14] [15 16 17 18 19]] [[20 21 22 23 24] [25 26 27 28 29] [30 31 32 33 34] [35 36 37 38 39]]] 120 (5, 4, 3, 2) [[[ 0 60] [ 20 80] [ 40 100]] [[ 5 65] [ 25 85] [ 45 105]]] [ 0 60 20 80 40 100 5 65 25 85 45 105 10 70 30 90 50 110 15 75 35]
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')
(64, 64, 35) [64 64 35] (42, 143360) (42, 143360) (143360,)
<matplotlib.text.Text at 0x4b591d0>
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)
0.0 0.0 True True False
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)
(64, 64, 35, 42) (64, 64, 35, 41) 143360 (143360, 41)
figure(num=None, figsize=(12,2))
plot(m2diff)
[<matplotlib.lines.Line2D at 0x76cd990>]
Select the biggest difference:
print np.argmax(m2diff)
18
Create a vector that you would put in the design matrix
nuisance = np.zeros(n_scans)
nuisance[np.argmax(m2diff)] = 1
print nuisance
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
Exercice : create a volume that will show a big difference - check that the covariate is indeed taking out this volume. This is the concept of a unit test !
# 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
(64, 64, 35, 42) (35, 42) (42, 143360) (42, 64, 64, 35) arr_demean_reshaped.shape (42, 64, 64, 35) 5.41365893912e-15
[<matplotlib.lines.Line2D at 0x8515b50>]
import nipy.algorithms.diagnostics as nads
tdiff = nads.time_slice_diffs(arr, time_axis=-1, slice_axis=-2)
nads.plot_tsdiffs(tdiff)
[<matplotlib.axes.AxesSubplot at 0xb611250>, <matplotlib.axes.AxesSubplot at 0xb945a10>, <matplotlib.axes.AxesSubplot at 0xb961ed0>, <matplotlib.axes.AxesSubplot at 0xbacf410>]
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)
(64, 64, 35, 42) (42, 143360)
plot(s)
[<matplotlib.lines.Line2D at 0x7f8f0c023250>]
plot(u[:,0])
[<matplotlib.lines.Line2D at 0x68b7410>]
tximg_centred = tximg - tximg.mean(axis=0)
print tximg_centred.shape
aTa = tximg_centred.dot(tximg_centred.T)
v,s,vh = lin.svd(aTa)
(42, 143360)
plot(s)
[<matplotlib.lines.Line2D at 0x6746690>]
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])
(143360,) (42, 143360) (64, 64, 35)
<matplotlib.image.AxesImage at 0x7f8f0d08a0d0>
First - write a test !
# 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)
/home/jb/data/qcpna/bolddata/smallbold.nii.gz
len(arr.shape)
stdts = std_diff_ts(arr)
plot(stdts)
[64 64 35 42] 143360 42
[<matplotlib.lines.Line2D at 0x7f8f04b15590>]
print np.argmax(stdts)
rk = np.argsort(stdts)
print argmax([stdts[rk[-1]-1],stdts[rk[-1]+1]])
17 1
bad_volume = np.argmax(stdts) + argmax([stdts[rk[-1]-1],stdts[rk[-1]+1]])