intro_to_4d.Rmd

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

# Working with four dimensional images, masks and functions

```{python}
#: Our usual imports
import numpy as np  # the Python array package
import matplotlib.pyplot as plt  # the Python plotting package
```

## Loading data with Nibabel

Nibabel allows us to `load` images. First we need to `import` the nibabel
library:

```{python}
import nibabel as nib
```

We are going to load the image with filename `ds114_sub009_t2r1.nii`.

First we use the `nipraxis` utility to find the file on the web and download
it to the computer we are working on.

```{python}
# Import the utility.
import nipraxis
# Use the fetch_file function to find the file and download
# it to this computer.
image_fname = nipraxis.fetch_file('ds114_sub009_t2r1.nii')
# Show the filename of the downloaded file.
image_fname
```

Next we `load` the image, to give me an “image” object:

```{python}
img = nib.load(image_fname)
type(img)
```

You can explore the image object with `img.` followed by the Tab key at the
Jupyter / IPython prompt.

## Saving memory

Because images can have large arrays, Nibabel does not load the image array
when you `load` the image, in order to save time and memory. The best way to
get the image array data is with the `get_fdata()` method.

```{python}
data = img.get_fdata()
type(data)
```

See [Nibabel images and
memory](https://nipy.org/nibabel/images_and_memory.html) for more detail.

## Four dimensional arrays: space + time

The image we have just loaded is a four-dimensional image, with a
four-dimensional array:

```{python}
data.shape
```

The first three axes represent three dimensional space. The last axis
represents time. Here the last (time) axis has length 173. This means that,
for each of these 173 elements, there is one whole three dimensional image. We
often call the three-dimensional images *volumes*. So we could say that this
4D image contains 173 volumes.

We have previously been taking slices over the third axis of a
three-dimensional image. We can now take slices over the fourth axis of this
four-dimensional image:

```{python}
# A slice over the final (time) axis
first_vol = data[:, :, :, 0]
```

This slice selects the first three-dimensional volume (3D image) from the 4D
array:

```{python}
first_vol.shape
```

You can use the ellipsis `...` when slicing an array. The ellipsis is a
short cut for “everything on the previous axes”. For example, these two have
exactly the same meaning:

```{python}
first_vol = data[:, :, :, 0]  # All rows, all columns, all planes, vol 0
first_vol_again = data[..., 0]  # The same thing, using the ellipsis
```

`first_vol` is a 3D image just like the 3D images you have already seen:

```{python}
# A slice over the third axis (dimension).
plt.imshow(first_vol[:, :, 14], cmap='gray')
```

## The voxel time course

We can think of this 4D data as a series of 3D volumes.  That is the way we
have been thinking of the 4D data so far:

```{python}
# This is slicing over the last (time) axis
first_vol = data[..., 0]
first_vol.shape
```

We can also index over the first three axes.  The first three axes in this
array represent space.

* The first axis goes from right to left (0 index value means right, 63 means
  left);
* The second axis goes from back to front (0 index value means back, 63 means
  front);
* The third axes goes from bottom to top (0 means bottom, 29 means top).

If you give me index values for these first three axes, you have given me a
*coordinate* in the first three axes of the array.

For example, you could give me an index tuple for the first three axes like
this: `(42, 32, 19)`.  The first index of 42 refers to a position towards
the left of the brain (> 31).  The second index of 32 refers to a position
almost in the center front to back.  The last index of 19 refers to a position
a little further towards the top of the brain – in this image.

This coordinate therefore refers to a particular part of the image:

```{python}
# Where is this in the brain?
# Don't worry about this line, we come to that soon.
mean_data = np.mean(data, axis=-1)
# Make a nice bright dot in the right place
mean_data[42, 32, 19] = np.max(mean_data)
plt.imshow(mean_data[:, :, 19], cmap='gray', interpolation='nearest')
```

If we slice into the data array with these coordinates, we will get a vector,
with the image value at that position (43, 32, 19), for every point in time:

```{python}
# This is slicing over all three of the space axes
voxel_time_course = data[42, 32, 19, :]
voxel_time_course.shape
```

![](images/volumes_and_voxel.png)

```{python}
# Do a plot of these values.
plt.plot(voxel_time_course)
plt.xlabel('Volume index')
plt.ylabel('Signal')
plt.title('Voxel time course');
```

We could call this a “voxel time course”.

## Numpy operations work on the whole array by default

Numpy operations like `mean` and `min`, and `max` and `std` operate on the
whole numpy array by default, ignoring any array shape. For example, here is
the maximum value for the **whole 4D array**:

```{python}
np.max(data)
```

This is exactly the same as:

```{python}
# Maximum when flattening the array to 1 dimension
np.max(data.ravel())
```

You can ask Numpy to operate over particular axes instead of operating over
the whole array. For example, this will generate a 3D image, where each array
value is the mean over the 173 values at the corresponding 3D position (the
mean across time):

```{python}
# Mean across the final (time) axis
mean_vol = np.mean(data, axis=3)
# Show the middle plane of this mean volume.
plt.imshow(mean_vol[:, :, 14], cmap='gray')
```