For this notebook, you'll need the image from the course website: http://nipy.bic.berkeley.edu/practical_neuroimaging/bold.nii.gz
Download it into the same directory as the notebook.
%pylab inline
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline]. For more information, type 'help(pylab)'.
Remember to be explicit about our imports
import numpy as np
import matplotlib.pylab as plt
Image loading library
import nibabel as nib
fname = 'bold.nii.gz'
img = nib.load(fname)
img.shape
(64, 64, 35, 165)
data = img.get_data()
data.dtype
dtype('int16')
data = data.astype(np.float32) # We need higher precision for our calculations
difference_volumes = np.diff(data, axis=3)
difference_volumes.shape
(64, 64, 35, 164)
sq_diff = difference_volumes ** 2
We want to take the mean of the squared difference across voxels within a volume. So we need to reshape to an array that is voxels by scans.
n_voxels = np.prod(sq_diff.shape[:-1])
n_scans = sq_diff.shape[-1]
n_voxels, n_scans
(143360, 164)
vox_by_scans = np.reshape(sq_diff, (n_voxels, n_scans))
vox_by_scans.shape
(143360, 164)
scan_means = np.mean(vox_by_scans, axis=0)
rms = np.sqrt(scan_means)
plt.plot(rms)
[<matplotlib.lines.Line2D at 0x7f7108969810>]
rms
array([ 21.37829399, 21.55786514, 22.90981674, 21.0092392 , 19.14093018, 19.34042931, 20.33923531, 19.72624779, 19.20626831, 19.1655426 , 19.31467247, 19.29933929, 19.64847565, 20.38537979, 21.09734154, 20.79466629, 21.05896378, 21.22260666, 21.51582146, 22.25740433, 22.14149094, 20.5478611 , 20.3484745 , 19.33082008, 19.5223484 , 19.62931633, 20.07916451, 19.83069229, 20.37201118, 20.68777466, 20.72026062, 20.72160721, 22.06051064, 21.7769928 , 20.5906868 , 20.39357185, 20.0847187 , 20.22652435, 20.22488976, 20.64841461, 20.31625938, 20.32613564, 20.82729912, 20.82919312, 20.49733734, 20.63652802, 20.71328354, 20.60631943, 20.63438606, 22.7435894 , 25.08402443, 22.15480232, 20.31813622, 21.49239731, 22.93136215, 21.85350037, 20.71713257, 20.75389862, 20.6148243 , 20.80734253, 21.40749741, 21.45947647, 20.65410042, 20.53489494, 26.78358841, 30.37088585, 21.19126511, 20.53188705, 20.66870308, 20.78804207, 20.68810272, 20.48601151, 20.27190208, 20.46998215, 21.73822212, 21.66768646, 22.29109764, 21.82066727, 20.9644146 , 20.64793777, 21.54236794, 21.66319084, 21.15878105, 20.92150307, 20.55589294, 21.12077332, 20.64767265, 20.47261047, 21.01286888, 20.71836472, 20.30692291, 20.72637939, 20.55831718, 20.69306755, 20.46246147, 20.15977669, 20.35247231, 20.33913803, 21.17718506, 21.22943687, 20.6655426 , 20.59701347, 20.62110519, 20.99123573, 23.87256622, 22.76289368, 20.32838821, 19.5613575 , 19.70575905, 19.7425251 , 21.26521111, 21.3801899 , 22.28580284, 21.09300232, 20.63578606, 20.55285263, 20.69499397, 20.53114891, 20.72611427, 20.62341309, 20.91604996, 20.96081543, 20.42824745, 20.47571564, 20.88507462, 20.90059662, 20.65011215, 20.88697243, 20.72612953, 20.63994026, 20.66607475, 20.59076118, 22.50296974, 21.29590034, 22.22260857, 21.04150963, 20.23121071, 19.65003395, 20.13300323, 20.39213181, 20.44101524, 20.34570312, 20.69584274, 20.68807602, 20.78573799, 20.53284645, 20.23099709, 20.24572372, 20.52599716, 20.39331245, 20.2400856 , 20.39010429, 20.2218895 , 20.51807022, 20.66083527, 20.5424633 , 21.45739174, 20.94854736, 20.21707535, 22.44255829, 24.59318733, 22.98219299, 22.61777115, 26.26057816], dtype=float32)
assert True
assert 1 == 1
We may need to get this RMS trace for more than one array, or image.
Make a function that takes an array, and returns the RMS time course.
def difference_rms(img_arr):
""" Your code goes here """
# You want to return the RMS from this function
result = difference_rms(data)
# assert np.all(result == rms) # This should work
worst = np.argmax(rms)
worst
65
worst_vol = difference_volumes[:, :, :, worst]
plt.gray()
plt.imshow(worst_vol[:, :, 17], interpolation='nearest')
<matplotlib.image.AxesImage at 0x7f7108bce350>
def montage(vol3d):
"""Returns a 2d image montage as an array, given a 3d volume array
The code here is a little fancy, no need to worry about the details for now. Just pretend it's FSL
"""
n_slices = vol3d.shape[-1]
n_cols = int(np.ceil(np.sqrt(n_slices)))
# Split and distrubute the 3D array over the last (slice) axis
rows = np.array_split(vol3d, n_cols, axis=2)
# ensure the last row is the same size as the others
rows[-1] = np.dstack((rows[-1], np.zeros(rows[-1].shape[0:2] + (rows[0].shape[2]-rows[-1].shape[2],))))
im = np.vstack([np.squeeze(np.hstack(np.dsplit(row, n_cols))) for row in rows])
return(im)
plt.imshow(montage(worst_vol), interpolation='nearest')
<matplotlib.image.AxesImage at 0x7f7108e204d0>
in_bad_order = np.argsort(rms)
in_bad_order
array([ 4, 9, 8, 11, 10, 23, 5, 24, 107, 25, 12, 137, 108, 7, 109, 27, 26, 36, 138, 95, 158, 152, 38, 37, 146, 136, 150, 147, 72, 90, 40, 52, 41, 106, 97, 6, 141, 22, 96, 28, 13, 151, 139, 149, 35, 122, 140, 94, 73, 87, 123, 71, 44, 153, 148, 117, 67, 145, 63, 155, 21, 115, 84, 92, 34, 131, 101, 47, 58, 102, 119, 48, 114, 45, 129, 86, 79, 39, 126, 62, 154, 100, 130, 68, 29, 143, 70, 93, 116, 142, 46, 56, 89, 30, 31, 118, 128, 91, 57, 144, 69, 15, 59, 42, 43, 124, 127, 125, 120, 83, 157, 121, 78, 103, 3, 88, 135, 16, 113, 14, 85, 82, 98, 66, 17, 99, 110, 133, 0, 111, 60, 156, 61, 53, 18, 80, 1, 81, 75, 74, 33, 77, 55, 32, 20, 51, 134, 19, 112, 76, 159, 132, 162, 49, 105, 2, 54, 161, 104, 160, 50, 163, 64, 65])
Maybe we want to remove the worst volume. How would I do that? How would I check I had done the right thing?
img_arr
def replace_vol(img_arr, vol_no):
""" Replace volume number `vol_no` with mean of vols either side """
# Copy the original data, ``img_arr``
# Take the mean of volumes 64, 66
# Replace volume 65 with the mean of 64, 66
fixed = replace_vol(data, 65) # Call our new function
assert not np.all(fixed == data) # We changed the array
left = data[:, :, :, 64]
right = data[:, :, :, 66]
mean_either_side = (left + right) / 2.0
# assert np.all(fixed[:, :, :, 65] == mean_either_side) # This should work
import awesome # Yes, I wrote it
fixed = awesome.replace_vol(data, 65)
new_rms = awesome.difference_rms(fixed)
plt.plot(rms)
plt.plot(new_rms)
[<matplotlib.lines.Line2D at 0x7f7108986190>]
What happens if I try and fix volume 0? Volume 164? What should happen?