# Working with four dimensional images, with some revision¶

## Loading data with nibabel¶

First let’s load a three-dimensional file, like the one we were looking at last week.

```
import nibabel as nib
```

First we `load`

the image, to give an “image” object:

```
img = nib.load('ds107_sub001_highres.nii')
```

You can explore the returned image object with `img.`

followed by TAB in the
IPython console.

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

method.

```
data = img.get_data()
type(data)
```

This gives:

```
numpy.core.memmap.memmap
```

The `memmap`

is a special type of array that saves memory, but
otherwise behaves the same as any normal array.

The data is a three dimensional array. We can show slices from this data with matplotlib.

## Four dimensional arrays - space + time¶

Images can also be four-dimensional. It’s easiest to think of four dimensional images as a stack of 3-dimensional images (volumes).

Here we load a four dimensional image. It is an functional MRI image, where the volumes are collected in sequence a few seconds apart.

```
img = nib.load('ds114_sub009_t2r1.nii')
data = img.get_data()
```

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

```
data.shape
```

```
(64, 64, 30, 173)
```

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:

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

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

```
first_vol.shape
```

```
(64, 64, 30)
```

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:

```
first_vol = data[:, :, :, 0]
first_vol_again = data[..., 0] # Using the ellipsis
```

`first_vol`

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

```
# A slice over the third dimension of a 3D image
plt.imshow(first_vol[:, :, 14], cmap='gray')
```

## Numpy operations work on the whole array by default¶

Numpy operations like `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:

```
np.max(data)
```

```
6793
```

This is exactly the same as:

```
# 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 variance over the 173 values at that 3D position (the variance across time):

```
# variance across the final (time) axis
var_vol = np.var(data, axis=3)
plt.imshow(var_vol[:, :, 14], cmap='gray')
```

## Indexing with boolean arrays¶

We covered this briefly a few classes back.

Let’s say we have an array like this:

```
arr = np.array([[0, 1, 3, 0], [5, 2, 0, 1]])
```

We can get a True / False (boolean) array to tell us whether these values are above some threshold:

```
tf_array = arr > 2
tf_array
```

```
array([[False, False, True, False],
[ True, False, False, False]], dtype=bool)
```

You can flip the True / False values with `~`

(bitwise not):

```
~tf_array
```

```
array([[ True, True, False, True],
[False, True, True, True]], dtype=bool)
```

We can use boolean arrays to *index* into the original array (or any array
with a suitable shape). This will select only the elements where the boolean
array is `True`

. The returned array may well have selected only a few
members from any particular row or column or (in general) higher axis, so the
returned array is always one-dimensional to reflect the loss of shape:

```
selected_elements = arr[tf_array]
selected_elements
```

This gives:

```
array([3, 5])
```

The returned array is shape `(2,)`

(one-dimensional).

We can use this to select values in our image as well. For example, if we
wanted to select only values less than 10 in `first_vol`

:

```
tf_lt_10 = first_vol < 10
vals_lt_10 = first_vol[tf_lt_10]
np.unique(vals_lt_10)
```

This gives:

```
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9], dtype=int16)
```