---
jupyter:
jupytext:
notebook_metadata_filter: all,-language_info
split_at_heading: true
text_representation:
extension: .Rmd
format_name: rmarkdown
format_version: '1.2'
jupytext_version: 1.13.7
kernelspec:
display_name: Python 3
language: python
name: python3
---
# Voxels and time
See also: [Reshaping, 4D to 2D](reshape_and_4d.Rmd).
```{python}
# 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.
```{python}
import nipraxis
bold_fname = nipraxis.fetch_file('ds114_sub009_t2r1.nii')
bold_fname
```
```{python}
import nibabel as nib
img = nib.load(bold_fname)
img.shape
```
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:
```{python}
n_voxels = np.prod(img.shape[:-1])
n_voxels
```
To calculate the mean across all voxels, for a single volume, we can do this:
```{python}
data = img.get_fdata()
first_vol = data[..., 0]
np.mean(first_vol)
```
To calculate the mean across voxels, we could loop across all
volumes, and calculate the mean for each volume:
```{python}
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)
```
We could also flatten the three voxel axes out into one long voxel axis, using
reshape – see: [Reshaping, 4D to 2D](reshape_and_4d.Rmd). 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:
```{python}
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)
```