2 - nifti.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7f248ff0-3fa5-4f85-8180-8c6e9cbbbe63",
   "metadata": {},
   "source": [
    "# Working with NIfTI images\n",
    "\n",
    "[NIfTI](https://nifti.nimh.nih.gov/) stands for *Neuroimaging Informatics Technology Initiative*, which is jointly sponsored by the US  National Institute of Mental Health and the National Institute of Neurological Disorders and Stroke. NIfTI defines a file format for neuroimaging data that is meant to meet the needs of the fMRI research community. In particular, NIfTI was developed to support inter-operability of tools and software through a common file format. Prior to NIfTI there were a few major fMRI analysis software packages, and each used a different file format. NIfTI was designed to serve as a common file format for all of these (and future) neuroimaging software packages. \n",
    "\n",
    "NIfTI was derived from an existing medical image format, called ANALYZE. ANALYZE was originally developed by the Mayo Clinic in the US, and was adopted by several neuroimaging analysis software packages in the 1990s. The ANALYZE header (where meta-data are stored) had extra fields that were not used, and NIfTI format basically expands on ANALYZE by using some of those empty fields to store information relevant to neuroimaging data. In particular the header stores information about the position and orientation of the images. This was a huge issue prior to NIfTI. In particular, there were different standards for how to store the order of the image data. For example, some software packages stored the data in an array that started from the most right, posterior, and inferior voxel, with the three spatial dimensions ordered right-to-left, posterior-to-anterior, and then inferior-to-superior. This is referred to as RPI orientation. Other packages that also used ANALYZE data stored the voxels in RAI format (with the second dimension going anterior-to-posterior) or LPI format (reversing left and right). This caused a lot of problems for researchers, especially if they wanted to try different analysis software, or use a pipeline that involved tools from different software packages. In some cases, this was just annoying (e.g., having to reverse the anterior-posterior dimension of an image). In other cases, it was confounding and potentially created erroneous results. This was especially true of the right-left (*x*) dimension. While it is immediately obvious when viewing an image which the front and back, and top and bottom, of the brain are, the left and right hemispheres are typically indistinguishable from each other, so a left-right swap could easily go undetected, potentially leading researchers to make completely incorrect conclusions about which side of the brain activation occurred on! The NIfTI format was designed to help prevent this by more explicitly storing orientation information in the header.\n",
    "\n",
    "Another improvement with the NIfTI format was to allow a single file. ANALYZE format requires two files, a header (with a `.hdr` extension) and the image data itself (`.img`). These files had to have the same name prior to the extension (e.g., `brain_image.hdr` and `brain_image.img`), and doubled the number of files in a directory of images, which created more clutter. NIfTI defines a single image file ending in a `.nii` extension. As well, NIfTI images can be compressed using a standard, open-source algorithm known as Gzip, which can significantly reduce file sizes and thus the amount of storage required for imaging data. Since neuroimaging data files tend to be large, this compression was an important feature. \n",
    "\n",
    "Although other file formats are still used by some software, NIfTI has become the most widely used standard for fMRI and other MRI research data file storage. Here we will learn how to convert a DICOM file to NIfTI format, which is typically the first step in an MRI research analysis pipeline, since most MRI scanners produce DICOM files, but the software researchers use to process their data reads NIFTI and not DICOM format."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "388b0765-126f-4ba6-b6b7-0d86f773d8b3",
   "metadata": {},
   "source": [
    "## Import packages\n",
    "\n",
    "Here we load in three new Python packages designed to work with NIfTI data:\n",
    "- `dicom2nifti` converst DICOM images to NIfTI format\n",
    "- `NiBabel` reads and converts between NIfTI and several other common neuroimaging file formats, including ANALYZE\n",
    "- `NiLearn` is primarily designed to provide statistical analysis and machine learning tools for neuroimaging data. However, it also provides a number of utilities for reading and writing NIfTI images, and working with and visualizing data\n",
    "\n",
    "As well we'll load SciPy's ndimage package, and Matplotlib"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bd0bae4",
   "metadata": {},
   "source": [
    "~~~python\n",
    "import dicom2nifti\n",
    "import nibabel as nib\n",
    "import nilearn as nil\n",
    "import scipy.ndimage as ndi\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6dc25a97-910e-42c5-a962-892cabba37cd",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "f2636ee4-0ee0-4427-a9f9-b2645d6cf576",
   "metadata": {},
   "source": [
    "We will use `dicom2nifti`'s `convert_directory()` function to convert the structural MRI images we worked with in the previous lesson from DICOM to NIfTI. We pass it the name of the folder in which the DICOM images are saved, and also instruct it to compress the resulting NIfTI file (to save space). We also use the `reorient=True` kwarg to force the image to be written in LAI orientation (i.e., starting with the most left, anterior, and inferior voxel), which ensures there is no ambiguity about the resulting NIfTI image.\n",
    "\n",
    "`convert_directory` does not take an argument for the output file name. Instead, it uses the name of the scan that was used when it was acquired on the MRI scanner. This might seem like a frustrating lack of control, however it does ensure that there are no user errors in the conversion process, that could result in mis-identified files. Here we will first list the contents of the `data` folder, then run `convert_directory`, then list the contents again to see the new NIfTI file and what it is named:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0c674f3",
   "metadata": {},
   "source": [
    "~~~python\n",
    "# List the contents of the data folder\n",
    "os.listdir('data')\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15fb156f-0910-4c90-98f9-2ce0e0738293",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "00a404ef",
   "metadata": {},
   "source": [
    "~~~python\n",
    "# Run 'convert_directory'\n",
    "dicom2nifti.convert_directory('data/DICOM', 'data', compression=True, reorient=True)\n",
    "\n",
    "# List the contents of the folder again\n",
    "os.listdir('data')\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc8bca8b-8614-4eaf-b8bf-b8551db063d0",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "1f270762-fa39-41e3-8c12-9d78f6a0ba90",
   "metadata": {},
   "source": [
    "So the new converted file is `4_sag_3d_t1_spgr.nii.gz`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62f8edff-dc64-4536-8c4c-099131d3f74c",
   "metadata": {},
   "source": [
    "## Load NIfTI file\n",
    "\n",
    "We use `NiBabel` to read in the NIfTI file we just created:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02ab6cc3",
   "metadata": {},
   "source": [
    "~~~python\n",
    "brain_vol = nib.load('data/4_sag_3d_t1_spgr.nii.gz')\n",
    "\n",
    "# What is the type of this object?\n",
    "type(brain_vol)\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56bac565-fa1d-4467-a1da-947fd9a37b08",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "f0c5ebef-feee-44fd-b2a2-dd1ccd55f4e2",
   "metadata": {},
   "source": [
    "## View metadata\n",
    "\n",
    "We can view the image's header by printing it (note that due to how the NiBabel `Nifti1Image` object is coded, we need to `print()` the header rather than just asking for it as a property):"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bd4d747",
   "metadata": {},
   "source": [
    "~~~python\n",
    "print(brain_vol.header)\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecc91a3e-aa9f-49d0-bbd6-bd66614a89de",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "39f82d91-ca0a-48fd-82d1-b162ad0bfee5",
   "metadata": {},
   "source": [
    "## Access data in the NIfTI object\n",
    "\n",
    "NiBabel's handling of the NIfTI format data is not quite as elegant as what we saw in the previous lesson. Rather than being able to access the data directly by referencing the name of the object (in this case, `brain_vol`), we need to use the method `get_fdata()` to do this (the \"f\" in this method name stands for \"floating point\", as this is the type of data it returns). We will assign the result of this to a new variable so that it's easy to work with. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "517b884b",
   "metadata": {},
   "source": [
    "~~~python\n",
    "# Assigning the data to a new variable\n",
    "brain_vol_data = brain_vol.get_fdata()\n",
    "\n",
    "# Getting the data type\n",
    "type(brain_vol_data)\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b26e91e-75e5-44a0-8e69-99112ed7d6df",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "f5f48fa5-7d8c-42d5-aaf1-79f1ab6bc626",
   "metadata": {},
   "source": [
    "We see that the data is a familiar NumPy array, and below we see the dimensions are identical to what we saw for this image in the previous lesson:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b222b260",
   "metadata": {},
   "source": [
    "~~~python\n",
    "brain_vol_data.shape\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f72dc94-f45e-423d-9028-5066a93aceaa",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "d1174a1c-fc2f-4df0-903d-6e6e8b411553",
   "metadata": {},
   "source": [
    "## Visualize a slice\n",
    "\n",
    "We can use `.plt.imshow()` as in the previous lesson:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46eaf858",
   "metadata": {},
   "source": [
    "~~~python\n",
    "plt.imshow(brain_vol_data[96], cmap='bone')\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4090ff3-cf95-4b1a-b612-7a69a18ee978",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "25d5bc0e-7088-4784-ae56-620a08e6d048",
   "metadata": {},
   "source": [
    "Note that our image is rotated, so use can use `ndi.rotate` to fix this:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed1fd579",
   "metadata": {},
   "source": [
    "~~~python\n",
    "plt.imshow(ndi.rotate(brain_vol_data[96], 90), cmap='bone')\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75aab96d-38fa-4fc1-83c3-ddfb77fb6235",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "1aa4df76-2df7-4d75-9cb6-e96c1b20bf44",
   "metadata": {},
   "source": [
    "## Plot a series of slices"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d94fbda6",
   "metadata": {},
   "source": [
    "~~~python\n",
    "fig_rows = 4\n",
    "fig_cols = 4\n",
    "n_subplots = fig_rows * fig_cols\n",
    "n_slice = brain_vol_data.shape[0]\n",
    "step_size = n_slice // n_subplots\n",
    "plot_range = n_subplots * step_size\n",
    "start_stop = int((n_slice - plot_range) / 2)\n",
    "\n",
    "fig, axs = plt.subplots(fig_rows, fig_cols, figsize=[10, 10])\n",
    "\n",
    "for idx, img in enumerate(range(start_stop, plot_range, step_size)):\n",
    "    axs.flat[idx].imshow(ndi.rotate(brain_vol_data[img, :, :], 90), cmap='gray')\n",
    "    axs.flat[idx].axis('off')\n",
    "        \n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d58cc45-4645-4742-a056-9a3322977d6c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "8be54919-b388-45fd-9352-3b149aea9039",
   "metadata": {},
   "source": [
    "## Plot with NiLearn\n",
    "\n",
    "While SciPy's ndimage module was designed for working with a wide variety of image types, NiLearn was designed to work with neuroimaging data specifically. As such, it's tools are a bit easier to use and more purpose-built for tasks that neuroimaging data scientists might want to perform. For example, we can plot the NiBabel NIfTI image object directly without first having to extract the data, using the `plot_img()` function from NiLearn's `plotting` module:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c10f8e6c",
   "metadata": {},
   "source": [
    "~~~python\n",
    "from nilearn import plotting\n",
    "\n",
    "plotting.plot_img(brain_vol)\n",
    "plt.show()\n",
    "~~~\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9858856-35d6-4f46-b545-08869e55dc3b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "e8c83f71-891e-4764-bdc2-00044731873d",
   "metadata": {},
   "source": [
    "One nice thing that we see is that since NiLearn is neuroimaging-aware, it explicitly adds labels to our plot showing us clearly which the left and right hemispheres are.\n",
    "\n",
    "NiLearn's plotting library uses Matplotlib, so we can use familiar tricks to do things like adjust the image size and colourmap:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55cbaac3",
   "metadata": {},
   "source": [
    "~~~python\n",
    "from nilearn import plotting\n",
    "\n",
    "fig, ax = plt.subplots(figsize=[10, 5])\n",
    "plotting.plot_img(brain_vol, cmap='gray', axes=ax)\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43d4c387-27ae-48a8-95af-eccf449b1db5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "1c3ca9bf-da00-4bcc-bf87-740e8b32e682",
   "metadata": {},
   "source": [
    "The `plot_img()` function also provides a variety of ways to display the brain, with much less code than we had to use when working with raw NumPy arrays and Matplotlib functions:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f96d771b",
   "metadata": {},
   "source": [
    "~~~python\n",
    "plotting.plot_img(brain_vol, display_mode='tiled', cmap='gray')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d13fc207-0623-430b-a417-8c4ca9742f26",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "7446ac35-77ec-40c8-82b3-b67a3d045e64",
   "metadata": {},
   "source": [
    "We can use the `cut_coords` kwarg to specify there to centre the crosshairs and \"cuts\" through the image that we visualize. In this image, the coordinates are relative to the *isocentre* of the MRI scanner — the centre of the magnetic field inside the scanner. The position of a person's head relative to this isocentre will vary from individual to individual, and scan to scan, due to variations in head size and the optimizations used by the MRI technician and scanner. But we can use the coorindates printed in the above image (which defaulted to the centre of the image volume) and some trial-and-error to get a different view through the brain:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e31f4ab4",
   "metadata": {},
   "source": [
    "~~~python\n",
    "plotting.plot_img(brain_vol, cmap='gray', cut_coords=(-45, 40, 0))\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d34702bb-8d13-45ce-9124-ec62292f1e17",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "30dc9066-a531-4318-90f5-4944d3dafbf1",
   "metadata": {},
   "source": [
    "`plot_img()` also has a few other ways to see multiple slices at once:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e47501d0",
   "metadata": {},
   "source": [
    "~~~python\n",
    "plotting.plot_img(brain_vol, display_mode='x', cmap='gray')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cedd7187-628d-4dc7-a4c2-b6b9987d1df3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "96d2cbf7",
   "metadata": {},
   "source": [
    "~~~python\n",
    "plotting.plot_img(brain_vol, display_mode='mosaic', cmap='gray')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2081e6f-cde5-48b6-a9ac-e743356fb4b6",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "86c55338-9a00-49a9-8853-87b6314054b5",
   "metadata": {},
   "source": [
    "## Smoothing\n",
    "\n",
    "NiLearn has its own function for applying Gaussian spatial smoothing to images as well. The only real difference from scipy.ndimage's `gaussian_filter()` function is that instead of specifying the smoothing kernel in standard deviations, we specify it in units of **full width half-maximum (FWHM)**. This is the standard way that most neuroimaging analysis packages specify smoothing kernel size, so it is preferable to SciPy's approach. As the term implies, FWHM is the width of the smoothing kernel, in millimetres, at the point in the kernel where it is half of its maximum height. Thus a larger FWHM value applies more smoothing."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db390b17",
   "metadata": {},
   "source": [
    "~~~python\n",
    "from nilearn import image\n",
    "\n",
    "fwhm = 4\n",
    "\n",
    "brain_vol_smth = image.smooth_img(brain_vol, fwhm)\n",
    "plotting.plot_img(brain_vol_smth, cmap='gray', cut_coords=(-45, 40, 0))\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "551a836f-cd44-4351-ab08-fd6b126271cc",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}