{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Working with MRI data in Python\n", "In this tutorial we will discuss how to interact with *Nifti* files — the file format used most in the MRI community — using the Python package *Nibabel*. We assume that you have experience with basic Python syntax and Numpy and Matplotlib functionality. If not, please consult [this page](https://lukas-snoek.com/NI-edu/other/python_recap.html) for an overview of resources.\n", "\n", "**What you'll learn**: after this week's lab ... \n", "\n", "* you know how to load, manipulate, and analyze Nifti-files\n", "\n", "**Estimated time to complete**: 1-2 hours" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What is \"nifti\"?\n", "Usually, every scanner manufacturer (Philips, Siemens, GE, etc.) have their own proprietary data format. For example, here at the University of Amsterdam, we have a Philips scanner and we often export our data as PAR/REC files. However, to streamline the development and use of neuroimaging software packages, the Neuroimaging InFormatics Technology Initiative came up with a new, standardized format, *nifti*, that most neuroimaging packages should support. As such, usually the first step in any (f)MRI preprocessing pipeline is to convert the scanner-specific files (e.g., PAR/RECs) to nifti. This is beyond the scope of this course, but if at one point you need to do this yourself, we highly recommed the [dcm2niix](https://github.com/rordenlab/dcm2niix) package!\n", "\n", "Throughout these tutorials you'll work with these nifti-files (also called nifti images), which you can recognize by their extension of *.nii* or its compressed version *.nii.gz*. This file-format is supported by most neuroimaging software packages (such as FSL, which you'll use later in this course!). \n", "\n", "However, we'd like to inspect and analyze nifti images in Python as well! *Nibabel* is an awesome Python package that allows us to read and load nifti images, and convert them to numpy arrays in a straightforward manner." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Let's load some other packages we need\n", "import os\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "import nibabel as nib # common way of importing nibabel" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll load in an example anatomical MRI scan (`anat.nii.gz`) from the current directory using nibabel below:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mri_file = 'anat.nii.gz'\n", "img = nib.load(mri_file)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note the type of `img`: the `Nifti1Image` class. This is a custom class just like a Numpy array, with its own attributes and methods!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(type(img))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For example, one attribute of `Nifti1Image` objects is the `shape` (which is similar to the numpy array attribute with the same name):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(img.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, the `shape` attribute is telling us that this is a 3D (anatomical) scan and has 240 voxels in the first dimension, 240 voxels in the second dimension, and 220 voxels in the third dimension." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The three parts of nifti images\n", "Nifti images can be roughly divided into three \"parts\":\n", "1. The header with metadata;\n", "2. The image data;\n", "3. The affine matrix\n", "\n", "All three parts are of course represented in nibabel's `Nifti1Image` class. Let's go through them one by one." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The header\n", "The header of nifti files contain metadata about the scan, such as the units of measurement, the voxel size, etc. In `Nifti1Images`, the header is an attribute:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# here, we're storing the header attribute in a new variable, hdr, for easy of use\n", "hdr = img.header" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Perhaps confusingly, the header is a custom object (a `Nifti1Header` object) as well, with its own methods and attributes. For example, it has a method called `get_zooms()`, which returns the voxel size (and optionally the sampling rate, if it's a fMRI file):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hdr.get_zooms() # it's a 1x1x1 mm MRI file!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another useful method is the `get_xyzt_units` which returns the units of the measurements (here: voxel size in millimeter and time in seconds):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hdr.get_xyzt_units()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's also load in a functional MRI file (`func.nii.gz`) to see the difference with an anatomical MRI file:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fmri_file = 'func.nii.gz'\n", "f_img = nib.load(fmri_file)\n", "print(f_img.shape)\n", "print(f_img.header.get_zooms())\n", "print(f_img.header.get_xyzt_units())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In case of fMRI files, the fourth dimension (almost) always represents the \"time\" dimension. So you can assume that a nifti image of an fMRI file has 4 dimensions, with the first three being the spatial dimensions (similar to the anatomical MRI file: $X \\times Y \\times Z$) and the last (fourth) being the time dimension ($T$).\n", "\n", "So for the above file, you can assume that it has 50 timepoints and has a sampling rate of 0.7 seconds (i.e., a new volume was scanned every 0.7 seconds).\n", "\n", "Moreover, you can infer that this file contains data from $80 \\times 80 \\times 44$ voxels with dimensions $2.7 \\times 2.7 \\times 2.97$ (in millimeters). To be honest, in practice, you won't deal a lot with the header (as you are generally aware of the dimensions/units of your data), so let's look at the *actual* data!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The data\n", "When we loaded in the data and created a `Nifti1Image` object, we actually didn't load in the *actual* data (i.e., voxel intensities) in memory yet! This is postponed because it's quite a memory-intensive operation: remember, loading in an fMRI scan with dimensions $80 \\times 80 \\times 44 \\times 50$ is effectively loading in over 14 million numbers (which is, in this case, more than 50MB in size)! By postponing loading the data, you can, for example, peek at the image dimensions without loading all the data in memory.\n", "\n", "Anyway, to actually load in the data, you can call the `get_fdata()` method, which will return a numpy array with the same dimensions as the image data. We'll take a look at the anatomical MRI data (`anat.nii.gz`):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "img_data = img.get_fdata()\n", "print(type(img_data)) # it's a numpy array!\n", "print(img_data.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Okay, so `img_data` is a 3D numpy array with ... what exactly? Let's check it out:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(img_data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It's just a bunch of numbers! This is important to remember: **like any image (like your standard png, jpeg, or gif), MRI scans are also just (3D or 4D) arrays of numbers**! The higher the number, the more signal the scanner recorded at that voxel. (It's actually a little more complex than this, but that's beyond the scope of this course!) \n", "\n", "Often, the absolute value of the signal intensities is not necessarily most or only important thing. The *relative differences* between different voxels in space or within voxels across time is also important. For example, for anatomical scans, apart from a high overall signal intensity, you often want a good *contrast* between white and gray matter (i.e., different signal intensities between voxels in these two tissue types). For functional MRI, apart from a clear tissue contrast and a high overall signal intensity, you often also want little (non-experimentally related) differences in signal intensity within a voxel across time.\n", "\n", "Anyway, when we printed this 3D array (of $240 \\times 240 \\times 220$), the notebook chose not to display all 12,672,000 numbers but instead truncated it. The reason it only printed zeros is because the first and last couple of numbers in each dimension likely doesn't contain any (signal related to) brain, just empty space!\n", "\n", "So, let's index a small $3 \\times 3 \\times 3$ patch of voxels in the middle of the brain to check what values these voxels contain:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mid_vox = img_data[118:121, 118:121, 108:111]\n", "print(mid_vox)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That's better! The exact values are not necessarily directly interpretable, i.e., you cannot say that the value 61978.46 is good (or bad), because the exact *scale* of the signal intensities (whether it goes from 0-100 or from 0-100,000) depends on the specific scanner hardware and specific nifti conversion software.\n", "\n", "Like any image, (f)MRI scans can be visualized by plotting the numbers and assigning colors to the numbers. Often we visualize brain images in brain and white, with higher signal intensities being brighter and lower signal intensities being darker. Importantly, remember that our data here cannot directly be plotted as a (2D) image, because our data (an anatomical scan) is 3D! However, we *can* just plot a single *slice* of the 3D volume, for example, the middle slice of our first voxel axis:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mid_slice_x = img_data[119, :, :]\n", "print(mid_slice_x.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can use matplotlib to plot this slice as an image using the `imshow` function you're seen before:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Note that the transpose the slice (using the .T attribute).\n", "# This is because imshow plots the first dimension on the y-axis and the\n", "# second on the x-axis, but we'd like to plot the first on the x-axis and the\n", "# second on the y-axis. Also, the origin to \"lower\", as the data was saved in\n", "# \"cartesian\" coordinates.\n", "plt.imshow(mid_slice_x.T, cmap='gray', origin='lower')\n", "plt.xlabel('First axis')\n", "plt.ylabel('Second axis')\n", "plt.colorbar(label='Signal intensity')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Alright, time to practice." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "