# 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>]
```

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)
```