Voxels and time#
See also: Reshaping, 4D to 2D.
# Import common modules
import numpy as np # the Python array package
import matplotlib.pyplot as plt # the Python plotting package
# Display array values to 6 digits of precision
np.set_printoptions(precision=4, suppress=True)
In this example, we calculate the mean across all voxels at each time point.
We’re working on ds114_sub009_t2r1.nii
. This is a 4D FMRI image.
import nipraxis
bold_fname = nipraxis.fetch_file('ds114_sub009_t2r1.nii')
bold_fname
'/home/runner/.cache/nipraxis/0.5/ds114_sub009_t2r1.nii'
import nibabel as nib
img = nib.load(bold_fname)
img.shape
(64, 64, 30, 173)
We want to calculate the mean across all voxels. Remember that a voxel is a pixel with volume, and refers to a position in space. Therefore we have this number of voxels in each volume:
n_voxels = np.prod(img.shape[:-1])
n_voxels
122880
To calculate the mean across all voxels, for a single volume, we can do this:
data = img.get_fdata()
first_vol = data[..., 0]
np.mean(first_vol)
414.40107421875
To calculate the mean across voxels, we could loop across all volumes, and calculate the mean for each volume:
n_trs = img.shape[-1]
means = []
for vol_no in range(n_trs):
vol = data[..., vol_no]
means.append(np.mean(vol))
plt.plot(means)
[<matplotlib.lines.Line2D at 0x7fef28150a60>]
![_images/34fdff0f1f341c167e0648659db33248010f219817c6aa577300f175deb3eae4.png](_images/34fdff0f1f341c167e0648659db33248010f219817c6aa577300f175deb3eae4.png)
We could also flatten the three voxel axes out into one long voxel axis, using
reshape – see: Reshaping, 4D to 2D. Then we can use the axis
parameter to
the np.mean
function to calculate the mean across voxels, in one
shot. This is “vectorizing”, where we take an operation that needed a loop,
and use array operations to do the work instead:
voxels_by_time = np.reshape(data, (n_voxels, n_trs))
means_vectorized = np.mean(voxels_by_time, axis=0)
# The answer is the same, allowing for tiny variations.
assert np.allclose(means_vectorized, means)