The scanner collected each volume slice by slice. That means that each slice corresponds to a different time.
We are going to create a new BOLD time series where we have corrected the data for this effect, so the data for each slice corresponds to our best guess at what the data would have looked like, if we had acquired all the data at the beginning of the volume.
This will involve interpolating in time.
%matplotlib inline
import os
import numpy as np
import matplotlib.pyplot as plt
# Set some plot defaults
plt.rcParams['image.cmap'] = 'gray'
plt.rcParams['image.interpolation'] = 'nearest'
import nibabel as nib
bold_fname = 'fixed_bold.nii.gz'
bold_img = nib.load(bold_fname)
The fixed image has one less volume than the original dataset:
bold_img.shape
(64, 64, 30, 172)
The scanner acquired these data with interleaved slices in time. That means that the scanner first acquired slice 0, then slice 2 up to slice 29, and then returned to acquire slice 1, slice 3 up to slice 29. This type of acquisition is often called "ascending interleaved".
Get the data as an array:
bold_data = bold_img.get_data()
bold_data.shape
(64, 64, 30, 172)
Let's get the data for the voxel at coordinate 43, 38, 0. The last 0 tells us this data is from the first (bottom) slice. Here is the voxel value at that coordinate, for the first volume (volume 0):
bold_data[43, 38, 0, 0]
196
Here is the voxel value for the same coordinate in the second volume:
bold_data[43, 38, 0, 1]
210
This is the voxel value for all volumes:
time_course_slice_0 = bold_data[43, 38, 0, :]
time_course_slice_0.shape
(172,)
Of course the first two values are the values we have already seen:
time_course_slice_0[:2]
array([196, 210], dtype=int16)
# Plot the time course from the voxel from slice 0
plt.plot(time_course_slice_0)
[<matplotlib.lines.Line2D at 0x3e60e50>]
Note that the x axis is in scans. To convert to seconds, we need to get the times for the start of each TR:
TR = 2.5 # time to acquire one scan
T = len(time_course_slice_0) # number of scans
times_slice_0 = np.arange(T) * TR
plt.plot(times_slice_0, time_course_slice_0)
plt.xlabel('Time (seconds)')
<matplotlib.text.Text at 0x3e745d0>
We can also get a time course for the equivalent voxel, but one slice up:
time_course_slice_1 = bold_data[43, 38, 1, :] # same voxel, slice 1
plt.plot(time_course_slice_1)
[<matplotlib.lines.Line2D at 0x432e850>]
Notice our x axis is still in TRs (scans). What times were these slice 1 values acquired? They were acquired half way through the TR, after the scanner had finished acquiring all the even slices (slice 0, slice 2, ....). So the times in seconds for these values are:
times_slice_1 = np.arange(T) * TR + TR / 2. # Add time offset for this slice
plt.plot(times_slice_1, time_course_slice_1)
plt.xlabel('Time (seconds)')
<matplotlib.text.Text at 0x40ccb10>
At the moment we can't see the offset in time between these slices; let's plot only the first eight values:
plt.plot(times_slice_0[:8], time_course_slice_0[:8], 'x:')
plt.plot(times_slice_1[:8], time_course_slice_1[:8], 'rx:')
plt.xlabel('Time (seconds)')
<matplotlib.text.Text at 0x40f6790>
The job of slice timing, is to find new values for the red line, which correspond to the times of acquisition of the blue line (the blue crosses). In order to do this, we can use linear interpolation. You can think of this as drawing a vertical line at $x=t$ where $t$ is the time of acqusition of the first slice. We then take the value from the red line where this $x=t$ line crosses the red line:
plt.plot(times_slice_0[:8], time_course_slice_0[:8], 'x:')
plt.plot(times_slice_1[:8], time_course_slice_1[:8], 'rx:')
ax = plt.gca()
min_y, max_y= ax.get_ylim()
for i in range(1, 8):
x = times_slice_0[i]
plt.plot([x, x], [min_y, max_y], 'k:')
plt.xlabel('Time (seconds)')
<matplotlib.text.Text at 0x47dbc50>
The value at $x$ for the straight line between any two of the points $(x_0, y_0)$ and $(x_1, y_1)$ on the red line is given by:
$$ y = y_0 + (y_1 - y_0) \frac{x_1 - x}{x_1 - x_0} $$You can probably convince yourself of this by drawing the line on a piece of paper, and working out the slope and intercept.
plt.plot(times_slice_0[:8], time_course_slice_0[:8], 'x:')
plt.plot(times_slice_1[:8], time_course_slice_1[:8], 'rx:')
ax = plt.gca()
min_y, max_y= ax.get_ylim()
for i in range(1, 8):
x = times_slice_0[i]
plt.plot([x, x], [min_y, max_y], 'k:')
x0 = times_slice_1[i-1]
x1 = times_slice_1[i]
y0 = time_course_slice_1[i-1]
y1 = time_course_slice_1[i]
y = y0 + (y1 - y0) * (x1 - x) / (x1 - x0)
plt.plot(x, y, 'kx')
plt.xlabel('Time (seconds)')
<matplotlib.text.Text at 0x433fb90>
So now we can just replace the original values from the red line (values for slice 1) with our best guess values if the slice had been taken at the same times as slice 0 (the beginning of the volume):
est_time_course_slice_1 = np.zeros(T)
for i in range(1, T): # Don't interpolate for first value
x = times_slice_0[i]
x0 = times_slice_1[i-1]
x1 = times_slice_1[i]
y0 = time_course_slice_1[i-1]
y1 = time_course_slice_1[i]
y = y0 + (y1 - y0) * (x1 - x) / (x1 - x0)
est_time_course_slice_1[i] = y
# Fill in the first value with something
est_time_course_slice_1[0] = est_time_course_slice_1[1]
plt.plot(times_slice_0[:8], time_course_slice_0[:8], 'x:')
plt.plot(times_slice_0[1:8], est_time_course_slice_1[1:8], 'rx:')
plt.xlabel('Time (seconds)')
<matplotlib.text.Text at 0x49ff710>
We do this for every slice except the first, and then make a new dataset containing the interpolated time courses for all slices. Now we have done slice timing correction.