1. pna-notebooks
  2. loading_images.ipynb

Here we are going to have a look at the raw image data using the nibabel package.

Windows installers for nibabel are here: https://pypi.python.org/pypi/nibabel

First we need to set up inline plotting

In [1]:
%pylab inline

Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.

This command did lots of imports for us, but to be explicit, I'm going to import some useful namespaces again

In [2]:
import os # Module with commands to interact with the operating system
In [3]:
import numpy as np # this is the Python package for handling arrays
In [4]:
import matplotlib.pyplot as plt # This is for making 2D plots like MATLAB

Next we are going to load an image into memory and do some math on it.

We first load the anatomical image for the first subject

In [5]:
anat_fname = 'ds107_sub001_highres.nii.gz'
anat_fname
Out[5]:
'ds107_sub001_highres.nii.gz'

What is an image? We have already seen that it is an array, with some metadata. Nibabel does the work of reading the metadata into a nice format, and getting the image data as an array.

We first import the nibabel package to allow us to load the image.

In [6]:
import nibabel as nib

Remember you can check what's in the package by doing tab completion. Try doing tab complete on "nib."

In [7]:
anat_img = nib.load(anat_fname)
anat_img
Out[7]:
<nibabel.nifti1.Nifti1Image at 0x106296a90>

Try tab completing on anat_img. (remember the dot at the end).

An image has a shape

In [8]:
anat_img.shape
Out[8]:
(256, 256, 192)

As we saw before, our image is 256 by 256 by 192. This means that the image consists of an array with dimensions (256, 256, 192), where the elements of the array are the voxel (3D pixel) intensity values.

In [9]:
anat_arr = anat_img.get_data() # get the data array
type(anat_arr)
Out[9]:
numpy.ndarray

An array has a shape. You can check the other attribute it has by tab completing on anat_arr.. Remember the trailing dot.

In [10]:
anat_arr.shape
Out[10]:
(256, 256, 192)

We can pick out a slice using slicing syntax rather similar to MATLAB. Remember the 0-based indexing though.

In [11]:
middle_slice = anat_arr[:, :, 96]
middle_slice.shape
Out[11]:
(256, 256)

This is just a 2D array of numbers, where the numbers are intensity values. We can look at it using matplotlib image display

In [12]:
plt.imshow(middle_slice)
Out[12]:
<matplotlib.image.AxesImage at 0x1063c3310>

This gave us an axial slice. It looks like this image is stored with left to right changing over the first axis of the array, front to back changing over the second axis of the array, and bottom to top changing over the third axis. How would I get a coronal slice? A sagittal slice?

The slice, like the whole of the array, give the intensity values at each location in the image:

In [13]:
middle_slice
Out[13]:
array([[0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       ..., 
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0]], dtype=int16)

We have lots of zeros at the outside, but in the centre... :

In [14]:
middle_slice[122:133, 122:133]
Out[14]:
array([[185, 183, 186, 195, 187, 194, 180, 172, 180, 192, 186],
       [182, 189, 185, 192, 189, 185, 176, 181, 183, 177, 171],
       [179, 180, 173, 183, 187, 179, 182, 182, 170, 163, 147],
       [174, 166, 171, 176, 169, 173, 178, 170, 164, 154, 117],
       [178, 176, 170, 174, 170, 168, 162, 158, 157, 142, 107],
       [187, 189, 175, 177, 173, 171, 170, 171, 166, 150, 135],
       [175, 175, 175, 181, 179, 167, 159, 167, 172, 164, 161],
       [171, 173, 179, 182, 187, 187, 180, 178, 175, 170, 178],
       [171, 175, 182, 187, 183, 182, 182, 182, 185, 182, 187],
       [170, 177, 183, 196, 196, 188, 185, 179, 186, 186, 186],
       [179, 195, 192, 184, 186, 183, 182, 178, 177, 173, 170]], dtype=int16)

What data type does the image have?

Now we have an array, we can do things like find the maximum, and minimum and mean value in the volume:

In [15]:
anat_arr.max()
Out[15]:
794
In [16]:
anat_arr.min()
Out[16]:
0
In [17]:
anat_arr.mean()
Out[17]:
32.011690696080528

We can also do useful things like looking at the histogram of values in the image. We first flatten the array to be a 1D vector, then take the histogram of that:

In [18]:
anat_arr_1d = anat_arr.ravel() # as 1D vector
tmp = plt.hist(anat_arr_1d, 100) # 100 histogram bins for nice detail

It looks like there are some outlier values. We can can knock of the top 5 percent and redo the histogram

In [19]:
percentile_95 = np.percentile(anat_arr_1d, 95)
clipped_vector = anat_arr_1d[anat_arr_1d <= percentile_95]
tmp = plt.hist(clipped_vector, 100)

An image also has meta-data. This is a NIfTI image, and it has NIfTI metadata, in the form of a header:

In [20]:
hdr = anat_img.get_header()
print hdr

<class 'nibabel.nifti1.Nifti1Header'> object, endian='<'
sizeof_hdr      : 348
data_type       : 
db_name         : 
extents         : 0
session_error   : 0
regular         : r
dim_info        : 0
dim             : [  3 256 256 192   1   1   1   1]
intent_p1       : 0.0
intent_p2       : 0.0
intent_p3       : 0.0
intent_code     : none
datatype        : int16
bitpix          : 16
slice_start     : 0
pixdim          : [-1.  1.  1.  1.  1.  0.  0.  0.]
vox_offset      : 352.0
scl_slope       : 1.0
scl_inter       : 0.0
slice_end       : 0
slice_code      : unknown
xyzt_units      : 10
cal_max         : 0.0
cal_min         : 0.0
slice_duration  : 0.0
toffset         : 0.0
glmax           : 0
glmin           : 0
descrip         : FSL4.0
aux_file        : 
qform_code      : aligned
sform_code      : aligned
quatern_b       : 0.0
quatern_c       : 1.0
quatern_d       : 0.0
qoffset_x       : 125.0
qoffset_y       : -144.0
qoffset_z       : -99.0
srow_x          : [  -1.    0.    0.  125.]
srow_y          : [   0.    1.    0. -144.]
srow_z          : [  0.   0.   1. -99.]
intent_name     : 
magic           : n+1

An image also has an affine. The affine is a 4 by 4 matrix that gives the relationship of the array data to real space (in millimeters, relative to the scanner).

In [21]:
aff = anat_img.get_affine()
aff
Out[21]:
array([[  -1.,    0.,    0.,  125.],
       [   0.,    1.,    0., -144.],
       [   0.,    0.,    1.,  -99.],
       [   0.,    0.,    0.,    1.]])

In this case you can see that this affine came from the srow fields in the header.

Specifically the affine specifies the relationship between the array coordinates from our (256, 256, 192) array, and millimeters relative to the center of the scanner. Let's say we want to know the position in millimeter space corresponding to position [0, 0, 0] in array coordinates. To get the mm coordinate, we matrix multiply the affine with the array coordinate.

But wait, the array coordinate has 3 values ([0, 0, 0]) and the affine has 4 columns - that won't work.

This is because the 4 by 4 array operates in homogenous coordinates - see http://en.wikipedia.org/wiki/Homogeneous_coordinates. Homogenous coordinates are a trick to allow the affine matrix to contain translations as well as rotations and zooms and shears. In practice, before we multiply our array coordinate, we have to append a 1 at the end of the 3 coordinate values, to make it work with our 4 by 4 affine.

So, we get the position like this:

In [22]:
array_coords = [0, 0, 0]
homogenous_array_coords = [0, 0, 0, 1]
mm_coordinates = aff.dot(homogenous_array_coords)
mm_coordinates
Out[22]:
array([ 125., -144.,  -99.,    1.])

That means that array position [0, 0, 0] corresponds to millimeter position 125 mm to the right of the scanner center, 144 mm behind the scanner center, and 99 mm below the scanner center.

Another example - the mm location of center of the array:

In [23]:
homogenous_array_coords = [127, 127, 46, 1] # note the extra 1 again for homogenous coordinates
mm_coordinates = aff.dot(homogenous_array_coords)
mm_coordinates
Out[23]:
array([ -2., -17., -53.,   1.])

Images can also be four dimensional, where the fourth dimension is usually time. The EPI images are usually 4D images.

In [24]:
bold_fname = 'ds107_sub001_run01.nii.gz'
bold_fname
Out[24]:
'ds107_sub001_run01.nii.gz'
In [25]:
bold_img = nib.load(bold_fname)
bold_img
Out[25]:
<nibabel.nifti1.Nifti1Image at 0x106589890>
In [26]:
bold_img.shape
Out[26]:
(64, 64, 35, 164)

This EPI run contains 164 volumes, each of size (64, 64, 35). We can look at individual volumes by slicing:

In [27]:
bold_arr = bold_img.get_data()
bold_vol0_arr = bold_arr[:, :, :, 0]
bold_vol0_arr.shape
Out[27]:
(64, 64, 35)
In [28]:
bold_middle_slice = bold_vol0_arr[:, :, 17]
plt.imshow(bold_middle_slice)
Out[28]:
<matplotlib.image.AxesImage at 0x1063ce0d0>

The 4D array allows us to do interesting things like making an image which is the mean across time:

In [29]:
time_course_mean = bold_arr.mean(axis=3) # the last axis (0-based numbering)
time_course_mean.shape
Out[29]:
(64, 64, 35)
In [30]:
plt.imshow(time_course_mean[:, :, 17])
Out[30]:
<matplotlib.image.AxesImage at 0x1065a8510>

Or - mabye of more interest - the standard deviation across time:

In [31]:
time_course_std = bold_arr.std(axis=3) # the last axis (0-based numbering)
middle_slice_std = time_course_std[:, :, 17]
plt.imshow(middle_slice_std)
Out[31]:
<matplotlib.image.AxesImage at 0x10264e350>

Aha - there's lot's of variation at the edge of the brain and in the ventricles.

What's the location of the voxel in this slice with the highest variance?

In [32]:
max_index = np.argmax(middle_slice_std) # flattens the array first
max_coord_2d = np.unravel_index(max_index, middle_slice_std.shape) # get 2D coordinate equivalent
max_coord_2d
Out[32]:
(50, 15)
In [33]:
max_coord_3d = max_coord_2d + (17,) # add the slice number as third coordinate
max_coord_3d
Out[33]:
(50, 15, 17)
In [34]:
time_course_std[max_coord_3d]
Out[34]:
137.59464537780846
In [35]:
middle_slice_std.max()
Out[35]:
137.59464537780846

We can plot voxel time courses:

In [36]:
max_std_over_time = bold_arr[max_coord_3d]
plt.plot(max_std_over_time)
Out[36]:
[<matplotlib.lines.Line2D at 0x107627f50>]

If you have nipy installed, you can try running the diagnostics. They use calculations that are almost identical to the ones we've already run here.

In [37]:
import nipy.algorithms.diagnostics as nads

tdiff = nads.time_slice_diffs(bold_arr, time_axis=-1, slice_axis=-2)
nads.plot_tsdiffs(tdiff)
Out[37]:
[<matplotlib.axes.AxesSubplot at 0x117ae5150>,
 <matplotlib.axes.AxesSubplot at 0x117b36410>,
 <matplotlib.axes.AxesSubplot at 0x117cb2390>,
 <matplotlib.axes.AxesSubplot at 0x117cd2510>]