# - compatibility with Python 3
from __future__ import print_function # print('me') instead of print 'me'
from __future__ import division # 1/2 == 0.5, not 0
# - show figures inside the notebook
%matplotlib inline
# - import common modules
import numpy as np # the Python array package
import matplotlib.pyplot as plt # the Python plotting package
# - set defaults for plotting
plt.rcParams['image.cmap'] = 'gray'
plt.rcParams['image.interpolation'] = 'nearest'
Load the image ds114_sub009_t2r1.nii with nibabel. Get the data:
# Load the image 'ds114_sub009_t2r1.nii' with nibabel
# Get the data array from the image
import nibabel as nib
img = nib.load('ds114_sub009_t2r1.nii')
data = img.get_data()
data.shape
As you remember, the first volume in this dataset is a lot different from the rest, and this will mess up our interpolation in time.
So, we need to remove the first volume from the data first, using slicing:
# Remove the first volume by slicing
fixed_data = data[..., 1:]
fixed_data.shape
We start off with example time-courses from the first and second slice.
Use slicing to get a z slice 0 time series for an example voxel at voxel coordinates (23, 19, 0).
Do the same for a z slice 1 time series from (23, 19, 1).
Plot these time series against volume number on the same graph:
# Slice out time series for voxel (23, 19, 0)
# Slice out time series for voxel (23, 19, 1)
# Plot both these time series against volume number, on the same graph
slice_0_ts = fixed_data[23, 19, 0, :]
slice_1_ts = fixed_data[23, 19, 1, :]
plt.plot(slice_0_ts)
plt.plot(slice_1_ts)
The scanner collected slices for these data in an "ascending interleaved" order. That is, the scanner first collected z slice 0, then z slice 2, up to z slice 28. It then went back to collect z slice 1, 3, 5 up to z slice 29.
That means the scanner started collecting slice 0 in each volume, at the beginning of the volume.
The TR (time to collect one volume) is 2.5 seconds.
# - the time between scans
TR = 2.5
Make a time vector, length 172, that corresponds to the start time in seconds of each volume. This also gives the slice 0 start times.
# Make time vector containing start times in second of each volume,
# relative to start of first volume.
# Call this `slice_0_times`
slice_0_times = np.arange(fixed_data.shape[-1]) * 2.5
slice_0_times
The scanner starts to collect z slice 1 exactly half way through the volume (half way through the TR). Make a new vector that is the start time of acquisition of slice 1.
# Make time vector containing start times in seconds of z slice 1,
# relative to start of first volume.
# Call this `slice_1_times`
slice_1_times = slice_0_times + 2.5 / 2
slice_1_times
Now plot the first 10 values for the slice 0 times, against the first 10 values of the slice 0 time series.
To the same plot for the first 10 values of the slice 1 times, against the first 10 values of the slice 1 time series.
Use the :+ line marker for the plots to get the actual position of the points, and dotted lines betweeen them.
# Plot first 10 values of slice 0 times against first 10 of slice 0 time series
# Plot first 10 values of slice 1 times against first 10 of slice 1 time series
# Use ':+' marker
plt.plot(slice_0_times[:10], slice_0_ts[:10], ':+')
plt.plot(slice_1_times[:10], slice_1_ts[:10], ':+')
Import InterpolatedUnivariateSpline from scipy.interpolate. Make a new linear (k=1) interpolation object for slice 1, with the slice 1 times and values.
# Import `InterpolatedUnivariateSpline` from `scipy.interpolate`
# Make a new linear (`k=1`) interpolation object for slice 1, with slice 1 times and values.
from scipy.interpolate import InterpolatedUnivariateSpline
interp = InterpolatedUnivariateSpline(slice_1_times, slice_1_ts, k=1)
Call the object you got with the slice 0 times, to get the estimated time series values for slice 1, if slice 1 had been collected at the same time as slice 0:
# Call interpolator with `slice_0_times` to get estimated values
slice_1_ts_est = interp(slice_0_times)
Repeat the plot of the first 10 values of the time series. This time, on the same plot, add the estimated values for slice 1, if they had been collected at the same time as slice 0. Use a black x for the estimated points (marker kx'):
# Plot first 10 values of slice 0 times against first 10 of slice 0 time series
# Plot first 10 values of slice 1 times against first 10 of slice 1 time series
# Plot first 10 values of slice 0 times against first 10 of interpolated slice 1 time series
plt.plot(slice_0_times[:10], slice_0_ts[:10], ':+')
plt.plot(slice_1_times[:10], slice_1_ts[:10], ':+')
plt.plot(slice_0_times[:10], slice_1_ts_est[:10], 'kx')
Make a new data matrix to contain the slice time corrected values for all voxels. We want to keep the values for z slice 0 unchanged, so make the new data matrix by copying the old one (and therefore getting the slice 0 values):
# Copy old data to a new array
new_data = fixed_data.copy()
Loop over every x voxel coordinate, and then loop over every y voxel coordinate.
For each x, y voxel coordinate:
- extract the time series at this x, y coordinate for slice 1;
- make a linear interpolator object with the slice 1 times and the extracted time series;
- resample this interpolator at the slice 0 times;
- put this new resampled time series into the new data at the same position
# loop over all x coordinate values
# loop over all y coordinate values
# extract the time series at this x, y coordinate for slice 1;
# make a linear interpolator object with the slice 1 times and the extracted time series;
# resample this interpolator at the slice 0 times;
# put this new resampled time series into the new data at the same position
for x in range(fixed_data.shape[0]):
for y in range(fixed_data.shape[1]):
time_series = fixed_data[x, y, 1, :]
interp = InterpolatedUnivariateSpline(slice_1_times, time_series, k=1)
new_series = interp(slice_0_times)
new_data[x, y, 1, :] = new_series
Now we need to do the same thing for all the z slices.
To do this, we want to construct an offset vector (call it time_offset) of length (number of z slices) such that adding the time_offset[z] to the acquisition time of the the first slice will give us the time of acquisition of slice z. The next few steps are to get to that time_offset vector.
First, make a new vector acquisition_order that is length 30, where acquisition_order[i] is the order of acquisition of slice index i. For example, the first 4 elements of acqusition_order should be 0, 15, 1, 16.
# Make acquisition_order vector, length 30, with values:
# 0, 15, 1, 16 ... 14, 29
acquisition_order = np.zeros(30)
acquisition_index = 0
for i in range(0, 30, 2):
acquisition_order[i] = acquisition_index
acquisition_index += 1
for i in range(1, 30, 2):
acquisition_order[i] = acquisition_index
acquisition_index += 1
acquisition_order
Divide the acquisition order vector by number of slices, and multiply by the TR, to get the time offset for each z slice, relative to the start of the scan:
# Divide acquisition_order by number of slices, multiply by TR
time_offsets = acquisition_order / 30 * TR
time_offsets
Now we can do our whole slice time correction, for every slice.
- For each z coordinate (slice index):
- Make a time vector by adding the slice time offset for this slice, to the slice_0 times. Call this the
slice_z_timesvector; - For each x coordinate:
- For each y coordinate:
- extract the time series at this x, y, z coordinate;
- make a linear interpolator object with the
slice_z_timesand the extracted time series; - resample this interpolator at the slice 0 times;
- put this new resampled time series into the new data at the same position
- For each y coordinate:
- Make a time vector by adding the slice time offset for this slice, to the slice_0 times. Call this the
# For each z coordinate (slice index):
## Make `slice_z_times` vector for this slice
### For each x coordinate:
#### For each y coordinate:
##### extract the time series at this x, y, z coordinate;
##### make a linear interpolator object with the `slice_z_times` and the extracted time series;
##### resample this interpolator at the slice 0 times;
##### put this new resampled time series into the new data at the same position
for z in range(fixed_data.shape[2]):
slice_z_times = slice_0_times + time_offsets[z]
for x in range(fixed_data.shape[0]):
for y in range(fixed_data.shape[1]):
time_series = fixed_data[x, y, z, :]
interp = InterpolatedUnivariateSpline(slice_z_times, time_series, k=1)
new_series = interp(slice_0_times)
new_data[x, y, z, :] = new_series
Congratulations - you have just done slice timing correction on this 4D image.