NiBabel¶

Neuroimaging data and file structures in Python¶

Christopher J Markiewicz¶

NeuroHackademy 2020¶

The goal of this presentation is to familiarize you with some broad classes of neuroimaging data that can be interacted with in Python, using the NiBabel library.

This document is intended to be viewed as a RISE presentation. It works fine as a notebook, but blocks with images may look strange because they are formatted to work as slides.

NiBabel¶

NiBabel is a low-level Python library that gives access to a variety of imaging formats, with a particular focus on providing a common interface to the various volumetric formats produced by scanners and used in common neuroimaging toolkits:

It also supports surface file formats:

Tractography files:

As well as the CIFTI-2 format for composite volume/surface data.

Installation¶

NiBabel is available on PyPI:

pip install nibabel

And conda-forge:

conda install -c conda-forge nibabel

Note: This notebook assumes NiBabel 3+, which requires a minimum Python version of 3.5.

Learning objectives¶

1. Be able to load and save different types of files in NiBabel
2. Become familiar with the SpatialImage API and identify its components
3. Understand the differences between array and proxy images
4. Acquire a passing familiarity with the structures of surface images, CIFTI-2 files, and tractograms

Basic I/O¶

Loading¶

Saving¶

All NiBabel images have a .to_filename() method:

nibabel.save will attempt to convert to a reasonable image type, if the extension doesn't match:

Some image types separate header and data into two separate images. Saving to either filename will generate both.

Serialization¶

Some NiBabel images can be serialized to and deserialized from byte strings, for performing stream-based I/O.

Images that save to single files can generally be serialized. NIfTI-1/2, GIFTI and MGH are currently supported.

Spatial Images¶

For MRI studies, neuroimaging data is typically acquired as one or more volumes, a regular grid of values. NiBabel represents these data as spatial images, a structure containing three things:

1. The image data array: a 3D or 4D array of image data
2. An affine matrix: 4x4 array relating voxel coordinates and "world" coordinates
3. Image metadata: usually a format-specific header

Many file types encode this basic structure. NiBabel will read any of ANALYZE (plain, SPM99, SPM2 and later), NIfTI-1, NIfTI-2, MINC 1.0, MINC 2.0, AFNI BRIK/HEAD, MGH, ECAT, DICOM and Philips PAR/REC, and expose a simple interface:

Spatial images have some properties that should be familiar from NumPy arrays:

The data array¶

The data is a simple NumPy array. It can be accessed, sliced and generally manipulated as you would any array:

Each location in the image data array is a voxel (pixel with a volume), and can be referred to with voxel coordinates (array indices).

This is a natural way to describe a block of data, but is practically meaningless with regard to anatomy.

NiBabel has a basic viewer that scales voxels to reflect their size and labels orientations.

The crosshair is focused at the origin $(0, 0, 0)$.

All of this information is encoded in the affine:

Affine transforms¶

The affine is a 4 x 4 numpy array. This describes the transformation from the voxel space (indices $(i, j, k)$) to a reference space (coordinates $(x, y, z)$). These coordinates are, by convention, distance in mm right, anterior and superior of an origin*.

For an unmodified image, this reference space typically refers to an origin in the isocenter of the imaging magnet, and the directions right, anterior and superior are defined assuming a subject is lying in the scanner, face up.

* If a file format uses an alternative convention, NiBabel converts on read/write, so affines are always RAS+.

The affine as a series of transformations¶

An affine transformation can be decomposed into translation, rotation, scaling (zooms) and shear transformations, applied right-to-left.

NiBabel provides functions for extracting information from affines:

• Orientation (or axis codes) indicates the direction an axis encodes. If increasing index along an axis moves to the right or left, the axis is coded "R" or "L", respectively.
• Voxel sizes (or zooms)
• Obliquity measures the amount of rotation from "canonical" axes.

You can also use it to answer specific questions. For instance, the inverse affine allows you to calculate indices from RAS coordinates and look up the image intensity value at that location.

Don't Panic¶

If you didn't follow all of the above, that's okay. Here are the important points:

1. Affines provide the spatial interpretation of the data.
2. NiBabel has some useful methods for working with them.

You'll go over this again with Noah Benson in Introduction to the Geometry and Structure of the Human Brain.

Matthew Brett's Coordinate systems and affines tutorial is an excellent resource.

Image headers¶

The image header is specific to each file format. Minimally, it contains information to construct the data and affine arrays, but some methods are useful beyond that.

• get_zooms() method returns voxel sizes. If the data is 4D, then repetition time is included as a temporal zoom.
• get_data_dtype() returns the numpy data type in which the image data is stored (or will be stored when the image is saved).

For images with fixed header structures, header fields are exposed through a dict-like interface.

The MGH header is similarly accessible, but its structure is quite different:

Other formats use text key-value pairs (PAR/REC, BRIK/HEAD) or keys in NetCDF (MINC 1.0) or HDF5 (MINC 2.0) containers, and their values are accessible by other means.

Often, we're not particularly interested in the header, or even the affine. But it's important to know they're there and, especially, to remember to copy them when making new images, so that derivatives stay aligned with the original image.

Creating new images¶

Reading data will only take you so far. In many cases, you will want to use the image to compute a new image in the same space. Such a function might take the form:

def my_function(image):
    orig_data = image.get_fdata()
    new_data = some_operation_on(orig_data)
    return image.__class__(new_data, affine=image.affine, header=image.header)

Note the image.__class__ ensures that the output image is the same type as the input. If your operation requires specific format features, you might use a specific class like nb.Nifti1Image.

For example, perhaps we want to save space and rescale our T1w image to an unsigned byte:

Data objects - Arrays and Array Proxies¶

Recall that the spatial image contains data, affine and header objects. When creating a new image, the data array is typically an array.

When loading a file, an ArrayProxy is used for most image types.

An array proxy is an object that knows how to access data, but does not read it until requested. The usual way to request data is get_fdata(), which returns all data as floating point:

get_fdata() provides a very predictable interface. When you need more control, you'll want to work with the ArrayProxy directly.

Converting array proxies to arrays¶

Array proxies are designed to be one step away from a numpy array:

Memory maps are arrays that remain on disk, rather than in RAM. This is only possible with uncompressed images.

We can also cast to any type we please, however unwisely.

Indexing and slicing array proxies¶

One of the primary motivations for an array proxy is to avoid loading data unnecessarily. Accessing the array proxy with array indices or slices will return the requested values without loading other data:

For example, this is useful for fetching a single volume from a BOLD series:

Slicing works with compressed data, as well. Install the indexed-gzip package for significant speedups with gzipped files.

Scaling array proxies¶

Several image types, including NIfTI, have a slope and intercept (scaling factors) in the header that allows them to extend the data range and/or precision beyond those of the on-disk type.

For example uint16 can take the values 0-65535. If a BOLD series includes values from 500-2000, we can calculate a slope that will allow us to utilize the 16 bits of precision to encode our desired values.

Let's create an image from some random values in our plausible range:

Save as uint16 and check its values:

We clearly lost some precision...

But what's going on?

The ArrayProxy keeps track of scaling factors:

The scaling is done automatically when the data is accessed, by slice or whole.

The ArrayProxy guarantees that the data has the intended value, but the type can vary based on the on-disk type and the values of scaling factors.

get_fdata() sweeps these details under the rug and always gives you the same type.

Don't Panic¶

If you didn't follow all of the above, that's okay. Here are the important points:

1. When in doubt, use img.get_fdata() will fetch all of the data, and it will always be a float
2. img.dataobj exists if you want to load only some data or control the data type
3. Both methods transparently scale data when needed

In the NiBabel docs, The image data array gives you an overview of both methods, and Images and memory has even more details.

Slicing images¶

Slicing array proxies is nice. Wouldn't it be nicer to keep track of the affine and header?

The slicer attribute provides an interface that allows you to apply slices to an image, and updates the affine to ensure that the spatial information matches.

Consider the T1-weighted image from earlier:

We can use the slicer to crop unneeded voxels in the left-right and inferior-superior directions:

Note the origin crosshair points to the same structure. The affines now differ in translation:

You can even downsample an image, and the zooms will reflect the increased distance between voxels.

Note that this is a bad idea in most circumstances because it induces aliasing.

The better approach would be to anti-alias and then slice:

For non-spatial dimensions, slices or indices may be used to select one or more volumes.

Aliasing considerations apply to time series as well, so be careful with down-sampling here, too.

Surface Images¶

Although the scanner samples data in three dimensions, some brain structures are better represented as a convoluted sheet than a volume. Data may be usefully resampled onto a cortical sheet, but in the process, it loses the intrinsic geometry of a 3D array.

To represent data on a surface, you need the following structures:

1. The surface mesh
   1. Vertices: a list of coordinates in world space
   2. Faces: a list of 3-tuples of indices into the coordinate list
2. The data array: a 1D or 2D array of values or vectors at each vertex

Unlike spatial images, these components are frequently kept in separate files.

Terminological note: You are likely to encounter multiple names for each of these. Vertices may be called coordinates or a point set. Faces are often called triangles. When a data array is plotted on the surface, it might be called a texture.

Image-specific interfaces¶

The two supported surface formats at present are GIFTI and FreeSurfer geometry files. Unlike spatial images, there is not yet a common interface for working with surface-based data.

FreeSurfer encodes a great deal of information in its directory structure and file names, allowing the necessary data arrays to have relatively simple formats.

GIFTI is more of an interchange format, and so has a rich metadata structure and can store different kinds of data.

Surface meshes¶

To load a surface mesh in a FreeSurfer directory, use nb.freesurfer.read_geometry:

NiBabel does not have any viewer for surfaces, but Nilearn has plotting utilities:

Let's do the same thing with GIFTI. Here, the file-level metadata is minimal:

The vertices and faces are stored as separate data arrays, each with their own metadata. These can be queried by NIfTI intent code:

The actual values are stored as the .data attribute:

This API can be cumbersome, if all you want is the NumPy arrays. Thanks to a project started at Neurohackademy 2019, the agg_data() (aggregate data) method will find the requested data array(s) and return just the data:

Surface-sampled data¶

For data sampled to the surface, FreeSurfer has a few data types:

GIFTIs will be GIFTIs. This one is a BOLD series sampled to the fsaverage5 surface:

Each time point is an individual data array with intent NIFTI_INTENT_TIME_SERIES. agg_data() will aggregate these into a single array.

We can plot the mean BOLD signal. This time we will use a FreeSurfer surface, which Nilearn knows what to do with:

Don't Panic¶

If you didn't follow all of the above, that's okay. Here are the important points:

1. Practical considerations make dealing with surfaces a multi-file affair
2. The data structures are common (vertices, faces, per-vertex data arrays), but the arrangement can vary
3. For working with FreeSurfer files, learn to traverse the directory and identify file types
4. For working with GIFTI, learn to query the file for its contents

Hopefully next year, we'll be able to talk about a common interface that will simplify things.

CIFTI-2¶

From Glasser, et al., 2016. doi:10.1038/nn.4361

CIFTI-2 is a file format intended to cover many use cases for connectivity analysis.

Files have 2-3 dimensions and each dimension is described by one of 5 types of axis.

• Brain models: each row/column is a voxel or vertex
• Parcels: each row/column is a group of voxels and/or vertices
• Scalars: each row/column has a unique name
• Labels: each row/column has a unique name and label table
• Series: each row/column is a point in a series (typically time series), which increases monotonically

For example, a "parcellated dense connectivity" CIFTI-2 file has two dimensions, indexed by a brain models axis and a parcels axis, respectively. The interpretation is "connectivity from parcels to vertices and/or voxels".

On disk, the file is a NIfTI-2 file with an alternative XML header as an extension, schematized here.

NiBabel loads a header that closely mirrors this structure, and makes the NIfTI-2 header accessible as a nifti_header attribute.

The Cifti2Header is useful if you're familiar with the XML structure and need to fetch an exact value or have fine control over the header that is written.

Most of the time, the Axis format will be more useful:

The simplest way to get a handle on CIFTI-2 data is to use it. Let's take an axis and a data block and repackage the voxels as a regular NIfTI-1 image:

Now we can extract the values on a surface vertex. This time, as a simple numpy array:

Finally, combine into a function that will break a CIFTI-2 matrix into a volume and two surface components:

Tractography¶

From DIPY documentation

Diffusion MRI uses estimates of the directional flow of water to detect white matter tracts. These tracts are represented as a series of points in world space.

NiBabel represents tractogram files as a collection of these things:

1. A tractogram, which is:
   1. A set of streamlines: each streamline is a series of coordinates describing a path
   2. Streamline-level data: a set of data arrays associated with each streamline
   3. Point-level data: a set of data arrays associated with each point of each streamline
2. An affine matrix: 4x4 array relating voxel coordinates of a reference image and world coordinates
3. File-level metadata: a format-specific header

Files of this type are considered TractogramFiles. Like SpatialImage, the goal is to provide a uniform interface where possible.

The two currently-supported file types are TCK (MRtrix) and TRK (TrackVis). We will load some very small ones that were created for testing:

Like a spatial image, printing the tractogram file gives an overview:

The streamlines API¶

The overall API is as follows:

A single streamline¶

When looking at an individual streamline, the objects are familiar: NumPy arrays and Python dictionaries.

All streamlines¶

When looking at all streamlines, the objects are custom. These deal with issues arising from the fact that each streamline may have a different number of coordinates, so N-dimensional arrays are not appropriate.

The dictionaries operate like normal dictionaries, providing NumPy array and ArraySequence values.

Each component is also indexable by integer:

As we can see, there are a few ways to get at the same data, based on whether you prefer to operate by streamline or by data element (such as fa).

Summary¶

You made it!

Learning objectives (reprise)¶

1. Be able to load and save different types of files in NiBabel
   • nb.load, nb.save, img.to_filename(), img.to_bytes(), img.from_bytes()
2. Become familiar with the SpatialImage API and identify its components
   • img.get_fdata()/img.dataobj
   • img.affine
   • img.header
3. Understand the differences between array and proxy images
   • Arrays are what it says on the tin.
   • Proxies conserve memory and automatically scale values
   • get_fdata() and dataobj provide tradeoffs between consistency and control
4. Acquire a passing familiarity with:
   • Surfaces: Meshes and per-vertex data, idiosyncratic
   • CIFTI-2: 2-3D NIfTI with some interesting axes; brainmodels can be decomposed into volume and surface components
   • Tractograms: Streamlines are paths in space, and data can be attached to the point or the path