Combined slice timing and motion correction¶
Nipy has a routine that combines the calculation of the motion between scans, and the slice time interpolation.
You can imagine that, as the motion correction moves the slices, we might want to take into account the fact that the slices in the images now correspond to different times of acquisition.
Nipy has a routine that does that by combining slice time correction with simultaneous motion correction.
See:
Roche, Alexis (2011) A four-dimensional registration algorithm with application to joint correction of motion and slice timing in fMRI. Medical Imaging, IEEE Transactions on; 30:1546--1554
import numpy as np
# - make numpy print numbers to 4dp only
np.set_printoptions(precision=4, suppress=True)
import matplotlib.pyplot as plt
%matplotlib inline
# - set gray colormap and nearest neighbor interpolation by default
plt.rcParams['image.cmap'] = 'gray'
plt.rcParams['image.interpolation'] = 'nearest'
We load the 4D image to which we will apply the algorithm:
import nipy
img = nipy.load_image('ds114_sub009_t2r1.nii')
To calculate the slice timing, we need the image TR. In our case, we can see the TR recorded in the db_name and pixdim fields of the header, but this is often not the case, and we will have to get the TR from some other source such as the original DICOM files or our lab notebook written when we acquired the data.
print(img.metadata['header'])
In fact nipy has taken the TR information out of the header, and put it into the 4th row, 4th column of the 5D affine matrix for this image:
img.affine
TR = img.affine[3, 3]
TR
Before we get going, we drop the first four scans because of the T1 effect (see the diagnostics notebook):
img = img[..., 4:]
Now we are ready to do the combined motion and slice-time correction.
import nipy.algorithms.registration as reg
We have to specify the times at which the slices were acquired, relative to the beginning of the volume. We can do this by passing in the times as an array, or by passing in a function that returns the times, given the TR and number of slices, but we can also select common schemes by passing in a string such as "ascending" (slices acquired sequentially from bottom to top). Here we've specified that this acquisition acquired the second slice in space first in time, and then alternate slices to the top slice, followed by the first slice in space, and alternate slices to the top slice. We call this asc_alt_2_1 - but see nipy.algorithms.slicetiming.timefuncs for a collection of these functions with their standard names.
slice_info=2 tells nipy that the third image axes is the slice axis - meaning that the first slice collected in time was img[:, :, 1, 0], the second slice collected was img[:, :, 3, 0] and so on.
realigner = reg.SpaceTimeRealign(img, TR,
slice_times='asc_alt_2_1',
slice_info=2)
The estimation step takes a fairly long time:
realigner.estimate()
The estimation step results in a list of transforms, one list for each run. We only have one run, so we have a list with a single entry:
transforms = realigner._transforms
len(transforms)
Each list has a transform object corresponding to every volume in the run. These transform objects give the transformation to match each scan to the first:
run_0_xforms = transforms[0]
len(run_0_xforms)
run_0_xforms[0]
affines = [t.as_affine() for t in run_0_xforms]
print(affines[0])
print(affines[1])
We can get the translations and rotations from the affines, like this:
translations = np.array([aff[:3, 3] for aff in affines])
translations.shape
plt.plot(translations[:, 0], label='x trans')
plt.plot(translations[:, 1], label='y trans')
plt.plot(translations[:, 2], label='z trans')
plt.legend()
It takes a little more math to get the rotation (Euler) angles out of the affines, but not much:
rotations = np.array([reg.aff2euler(aff) for aff in affines])
plt.plot(rotations[:, 0], label='pitch')
plt.plot(rotations[:, 1], label='roll')
plt.plot(rotations[:, 2], label='yaw')
plt.legend()
Finally, we can generate a new 4D image with the estimated slice time and motion effects removed, using the resample method of the realigner object:
resampled = realigner.resample()
resampled