\(\newcommand{L}[1]{\| #1 \|}\newcommand{VL}[1]{\L{ \vec{#1} }}\newcommand{R}[1]{\operatorname{Re}\,(#1)}\newcommand{I}[1]{\operatorname{Im}\, (#1)}\)
Four dimensions exercise¶
- For code template see:
four_dimensions_code.py
; - For solution see: Four dimensions solution.
>>> #: Our usual imports
>>> import numpy as np # the Python array package
>>> import matplotlib.pyplot as plt # the Python plotting package
We are going to load a four-dimensional (X, Y, Z, t) BOLD image called
ds107_sub012_t1r2.nii
. Download the file using the link. Import
the nibabel
module, and load the image with nibabel to create an image
object.
>>> #- Load image object using nibabel
>>> #- Turn off nibabel memory mapping.
In the usual way get the array data from this image and show the image shape.
>>> #- Get image array data from image object
Select the 10th volume (time index 9) from 4D image data array, by slicing over the last dimension. What shape is it?
>>> #- Get the 10th volume and show shape
Get the standard deviation across all voxels in this 3D volume:
>>> #- Standard deviation across all voxels for 10th volume
Loop over all volumes in the 4D image and store the standard deviation for each volume in a list. Plot these to see if there are any volumes with particularly unusual standard deviation.
>>> #- Get standard deviation for each volume; then plot the values
SPM uses a measure called “global signal” to give a very rough estimate of the average signal value within the brain. The idea is that you need a threshold to select voxels with signal high enough to be inside the brain, and not in the air around the brain, which has very low (near 0) signal.
The algorithm goes like this.
- Get a single 3D volume V;
- Calculate the mean signal of the voxels in V; call that M;
- Make a threshold T where T = M / 8.
- Select all voxel values in V that are greater than T; call these values W;
- Return the mean of all values in W.
In the SPM code, the algorithm is implemented in a MATLAB function called
spm_global
.
I used the MATLAB script get_global_signals.m
to run the
spm_global
MATLAB function on the volumes of ds107_sub012_t1r2.nii
.
The script saved the SPM global values to a text file
global_signals.txt
. The first four lines of the
global_signals.txt
file look like this:
376.53
375.75
375.26
376.01
Read these global values calculated by SPM into a list, and plot the values.
>>> #- Read global signal values calculated by SPM, and plot
Now implement the algorithm above to recalculate the SPM global signal for the first volume (volume index zero). Hint: you will likely need to index using a boolean (mask) array. Remember, the steps are:
- Get a single 3D volume V;
- Calculate the mean signal of the voxels in V; call that M;
- Make a threshold T where T = M / 8.
- Select all voxel values in V that are greater than T; call these values W;
- Get the mean of all values in W.
You should get the same value as SPM - the first value you read from
global_signals.txt
.
>>> #- Apply algorithm for SPM global calculation to first volume
Make a function called spm_global
that accepts a 3D array as input, and
returns the global signal using the SPM algorithm. Call that function on the
first volume to show that it is working (as in print(spm_global(data[:, :,
:, 0]))
).
>>> #- Make an `spm_global` function that accepts a 3D array as input,
>>> #- and returns the global mean for the volume according to the SPM
>>> #- algorithm
Make a function called get_spm_globals
that accepts an image filename as
an argument. The function will load the image, get the array data for the
image, use your new spm_global
function calculate the global value for
each volume, and return these values as a list. Finally, show this is working
by plotting the values for the ds107_sub012_t1r2.nii
image with something
like:
all_globals = get_spm_globals('ds107_sub012_t1r2.nii')
plt.plot(all_globals)
>>> #- Write a function `get_spm_globals` that returns the global values
>>> #- for each volume