1. pna2015
  2. day12
In [1]:
# - compatibility with Python 3
from __future__ import print_function  # print('me') instead of print 'me'
from __future__ import division  # 1/2 == 0.5, not 0
In [2]:
# - show figures inside the notebook
%matplotlib inline
In [3]:
# - import common modules
import numpy as np  # the Python array package
import matplotlib.pyplot as plt  # the Python plotting package
In [4]:
# - set gray colormap and nearest neighbor interpolation by default
plt.rcParams['image.cmap'] = 'gray'
plt.rcParams['image.interpolation'] = 'nearest'
In [5]:
# - import numpy.linalg with a shorter name
import numpy.linalg as npl

Comparing arrays

A comparison between two arrays returns the elementwise result of the comparison:

In [6]:
arr1 = np.array([[0, 1, 2],
                [3, 4, 5]])
arr2 = np.array([[1, 1, 2],
                 [3, 4, 6]])
arr1 == arr2
Out[6]:
array([[False,  True,  True],
       [ True,  True, False]], dtype=bool)

Sometimes we want to know if two arrays are equal, in the sense that all the elements of the two arrays are equal to each other. For this we use np.all. np.all accepts an array as input, and returns True if all the elements are non-zero (evaluate to True if forced to a True / False value):

In [7]:
np.all([True, True, True])
Out[7]:
True
In [8]:
np.all([True, False, True])
Out[8]:
False

To ask whether all the elements in two arrays are equal, we can pass the result of the elementwise comparison to np.all:

In [9]:
np.all(arr1 == arr2)
Out[9]:
False
In [10]:
arr3 = arr1.copy()
np.all(arr1 == arr3)
Out[10]:
True

Sometimes we want to know if any of the values in an array are non-zero (evaluate to True if forced to a True / False value). Enter np.any:

In [11]:
np.any([False, False, False])
Out[11]:
False
In [12]:
np.any([False, False, True])
Out[12]:
True
In [13]:
np.any(arr1 == arr2)
Out[13]:
True
In [14]:
np.any(arr1 != arr3)
Out[14]:
False

Formatting output when printing

We sometimes want to print numbers or strings with particular formatting. For example, we might want to print a floating point number to only 2 digits of precision. One good way to do this is using the Python curly brace string formatting. It works like this:

In [15]:
'Here is a {} floating point number'.format(3.33333)
Out[15]:
'Here is a 3.33333 floating point number'

The curly braces indicate that the number should go at that position in the string. Calling the format method with the number as an argument makes the number into a string and replaces the curly braces with the number represented as a string.

You can use multiple values in the formatting, and mix strings and numbers:

In [16]:
'Here is a number - {}; here is a string - {}'.format(2.0, 'Hello')
Out[16]:
'Here is a number - 2.0; here is a string - Hello'

This system allows us to give formatting instructions for things like numbers, by using a : inside the curly braces, followed by the formatting instructions. Here is an example of printing a fixed precision floating point value (f) with exactly 4 digits after the decimal point:

In [17]:
'A formatted number - {:.4f}'.format(.2)
Out[17]:
'A formatted number - 0.2000'

See the Python string formatting docs for more details and examples.

Comparing floats and floating point error

There are an infinite number of real - er - numbers, for example, there are an infinite number of real numbers between 0 and 1.

The computer stores real numbers as floating point numbers. The typical format the computer uses for floating point numbers uses 8 bytes (64 bits of memory) per number. That means the computer can only represent some of the numbers between 0 and 1 (for example).

One number that standard floating point numbers cannot represent exactly is 0.1. Here we print the value with 20 digits after the decimal point to show that:

In [18]:
print('{:.20f}'.format(0.1))

0.10000000000000000555

This means that when we do mathematical operations on floating point numbers, the exact result of the operation may be a number that the computer cannot represent exactly. In that case the software implementing the calculation can at best give us the nearest floating point number that it can represent:

In [19]:
print('{:1.20f}'.format(1 / 10))

0.10000000000000000555

For this reason, sometimes the results of doing floating point calculations are not exactly accurate. This is called floating point error. In the case above the floating point error is very small - 0.00000000000000000555.

Floating point error means that we have to be careful comparing results of floating point calculations, because the result may be off by some small amount, due to floating point error:

In [20]:
c = np.cos(1)
# We hope that the inverse cos should give back the original number, but
np.arccos(c) == 1
Out[20]:
False
In [21]:
np.arccos(c)
Out[21]:
0.99999999999999989

In this case, we probably want to check that the result is equal within some reasonable error due to the limited precision of floating point. For example, np.allclose checks whether two numbers are almost equal within some default tolerance levels:

In [22]:
np.allclose(np.arccos(c), 1)
Out[22]:
True

Making and saving new images in nibabel

We often want to do some processing on an image, then save the processed image back to an image file on disk.

When we load an image from disk, we get back an image object. When we load a NIfTI .nii image, we get an image object of type Nifti1Image.

In [23]:
import nibabel as nib
img = nib.load('ds114_sub009_highres.nii')
type(img)
Out[23]:
nibabel.nifti1.Nifti1Image

Maybe we were worried about some very high values in the image, and we wanted to clip them down to a more reasonable number:

In [24]:
data = img.get_data()
data.max()
Out[24]:
memmap(3237.0, dtype=float32)

We might consider clipping the top 5 percent of voxel values:

In [25]:
data = img.get_data()
top_95_thresh = np.percentile(data, 95)
top_95_thresh
Out[25]:
722.0
In [26]:
new_data = data.copy()
new_data[new_data > top_95_thresh] = top_95_thresh
new_data.max()
Out[26]:
memmap(722.0, dtype=float32)

We can make a new Nifti1Image by constructing directly, passing the new data, the image affine, and (optionally) a template header for the image:

In [27]:
clipped_img = nib.Nifti1Image(new_data, img.affine, img.header)
type(clipped_img)
Out[27]:
nibabel.nifti1.Nifti1Image

The nib.Nifti1Image routine copies and adapts the passed header to the new image data shape, and affine.

In [28]:
# Show the original data array shape from the original header
img.header.get_data_shape()
Out[28]:
(256, 156, 256)
In [29]:
# Here we construct a new empty header
empty_header = nib.Nifti1Header()
empty_header.get_data_shape()
Out[29]:
(0,)

If we make a new image with this header, the constructor routine fixes the header to have the correct shape for the data array:

In [30]:
another_img = nib.Nifti1Image(new_data, img.affine, empty_header)
another_img.header.get_data_shape()
Out[30]:
(256, 156, 256)

We can save our new image with nib.save:

In [31]:
nib.save(clipped_img, 'clipped_image.nii')

This image has the clipped data:

In [32]:
nib.load('clipped_image.nii').get_data().max()
Out[32]:
memmap(722.0, dtype=float32)

Sometimes, the TR in a NIfTI image is in the header

The NIfTI standard suggests putting the TR of a functional image, into the voxel dimension field of the header.

You can get the voxel (plus TR) dimensions with the get_zooms method of the header object:

In [33]:
func_img = nib.load('ds114_sub009_t2r1.nii')
header = func_img.header
header.get_zooms()
Out[33]:
(4.0, 4.0, 4.0000162, 2.5)

In this case, the image spatial voxel sizes are (4 by 4 by 4) millimeters, and the TR is 2.5 seconds.

In fact these values come from the NIfTI header field called pixdim:

In [34]:
print(header)

<class 'nibabel.nifti1.Nifti1Header'> object, endian='<'
sizeof_hdr      : 348
data_type       : 
db_name         : ?TR:2500.000 TE:50
extents         : 0
session_error   : 0
regular         : r
dim_info        : 0
dim             : [  4  64  64  30 173   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.          4.          4.          4.00001621  2.5         1.          1.
  1.        ]
vox_offset      : 0.0
scl_slope       : nan
scl_inter       : nan
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           : 255
glmin           : 0
descrip         : 
aux_file        : 
qform_code      : scanner
sform_code      : scanner
quatern_b       : 0.0
quatern_c       : 0.9958999753
quatern_d       : -0.0904608666897
qoffset_x       : 124.244003296
qoffset_y       : -103.449684143
qoffset_z       : -33.4928588867
srow_x          : [  -4.           0.           0.         124.2440033]
srow_y          : [   0.            3.93453455    0.72072268 -103.44968414]
srow_z          : [  0.          -0.72071975   3.93455052 -33.49285889]
intent_name     : 
magic           : n+1

Unfortunately, it is common for people writing NIfTI images not to write this information correctly into the header, so we have to be careful, and very suspicious, if the TR value appears to be 0 or 1.