{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# NIfTI images and python\n",
    "\n",
    "The [`nibabel`](http://nipy.org/nibabel/) module is used to read and write NIfTI\n",
    "images and also some other medical imaging formats (e.g., ANALYZE, GIFTI,\n",
    "MINC, MGH).  `nibabel` is included within the FSL python environment.\n",
    "\n",
    "\n",
    "Building upon `nibabel`, the\n",
    "[`fslpy`](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/) library\n",
    "contains a number of FSL-specific classes and functions which you may find\n",
    "useful. But let's start with `nibabel` - `fslpy` is introduced in a different\n",
    "practical (`advanced_topics/08_fslpy.ipynb`).\n",
    "\n",
    "\n",
    "## Contents\n",
    "\n",
    "* [Reading images](#reading-images)\n",
    "* [Header info](#header-info)\n",
    " * [Voxel sizes](#voxel-sizes)\n",
    " * [Coordinate orientations and mappings](#orientation-info)\n",
    "* [Writing images](#writing-images)\n",
    "* [Exercise](#exercise)\n",
    "\n",
    "---\n",
    "\n",
    "<a class=\"anchor\" id=\"reading-images\"></a>\n",
    "## Reading images\n",
    "\n",
    "It is easy to read an image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import nibabel as nib\n",
    "import os.path as op\n",
    "filename =  op.expandvars('${FSLDIR}/data/standard/MNI152_T1_1mm.nii.gz')\n",
    "imobj = nib.load(filename, mmap=False)\n",
    "\n",
    "# display header object\n",
    "imhdr = imobj.header\n",
    "\n",
    "# extract data (as a numpy array)\n",
    "imdat = imobj.get_fdata()\n",
    "print(imdat.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Make sure you use the full filename, including the `.nii.gz` extension.\n",
    "> `fslpy` provides FSL-like automatic file suffix detection though.\n",
    "\n",
    "> We use the `expandvars()` function above to insert the FSLDIR\n",
    "> environmental variable into our string. This function is\n",
    "> discussed more fully in the file management practical.\n",
    "\n",
    "Reading the data off the disk is not done until `get_fdata()` is called.\n",
    "\n",
    "> Pitfall:\n",
    ">\n",
    "> The option `mmap=False` disables memory mapping, which would otherwise be\n",
    "> invoked for uncompressed NIfTI files but not for compressed files. Since\n",
    "> some functionality behaves differently on memory mapped objects, it is\n",
    "> advisable to turn this off unless you specifically want it.\n",
    "\n",
    "Once the data is read into a numpy array then it is easily manipulated.\n",
    "\n",
    "> The `get_fdata` method will return floating point data, regardless of the\n",
    "> underlying image data type. If you want the image data in the type that it\n",
    "> is stored (e.g. integer ROI labels), then use\n",
    "> `imdat = np.asanyarray(imobj.dataobj)` instead.\n",
    "\n",
    "---\n",
    "\n",
    "<a class=\"anchor\" id=\"header-info\"></a>\n",
    "## Header info\n",
    "\n",
    "There are many methods available on the header object - for example, look at\n",
    "`dir(imhdr)` or `help(imhdr)` or the [nibabel webpage about NIfTI\n",
    "images](http://nipy.org/nibabel/nifti_images.html)\n",
    "\n",
    "<a class=\"anchor\" id=\"voxel-sizes\"></a>\n",
    "### Voxel sizes\n",
    "\n",
    "Dimensions of the voxels, in mm, can be found from:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "voxsize = imhdr.get_zooms()\n",
    "print(voxsize)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a class=\"anchor\" id=\"orientation-info\"></a>\n",
    "### Coordinate orientations and mappings\n",
    "\n",
    "Information about the NIfTI qform and sform matrices can be extracted like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sform = imhdr.get_sform()\n",
    "sformcode = imhdr['sform_code']\n",
    "qform = imhdr.get_qform()\n",
    "qformcode = imhdr['qform_code']\n",
    "print(qformcode)\n",
    "print(qform)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also get both code and matrix together like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "affine, code = imhdr.get_qform(coded=True)\n",
    "print(affine, code)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "<a class=\"anchor\" id=\"writing-images\"></a>\n",
    "## Writing images\n",
    "\n",
    "\n",
    "If you have created a modified image by making or modifying a numpy array then\n",
    "you need to put this into a NIfTI image object in order to save it to a file.\n",
    "The easiest way to do this is to copy all the header info from an existing\n",
    "image like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newdata = imdat * imdat\n",
    "newhdr = imhdr.copy()\n",
    "newobj = nib.nifti1.Nifti1Image(newdata, None, header=newhdr)\n",
    "nib.save(newobj, \"mynewname.nii.gz\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where `newdata` is the numpy array (the above is a random example only) and\n",
    "`imhdr` is the existing image header (as above).\n",
    "\n",
    "> It is possible to also just pass in an affine matrix rather than a\n",
    "> copied header, but we *strongly* recommend against this if you are\n",
    "> processing an existing image as otherwise you run the risk of\n",
    "> swapping the left-right orientation.  Those that have used\n",
    "> `save_avw` in matlab may well have been bitten in this way in the\n",
    "> past.  Therefore, copy a header from one of the input images\n",
    "> whenever possible, and just use the affine matrix option if you are\n",
    "> creating an entirely separate image, like a simulation.\n",
    "\n",
    "If the voxel size of the image is different, then extra modifications will be\n",
    "required.  Take a look at the `fslpy` practical for some extra image\n",
    "manipulation options, including cropping and resampling\n",
    "(`advanced_topics/08_fslpy.ipynb`).\n",
    "\n",
    "---\n",
    "\n",
    "\n",
    "<a class=\"anchor\" id=\"exercises\"></a>\n",
    "## Exercise\n",
    "\n",
    "\n",
    "Write some code to read in a 4D fMRI image (you can find one\n",
    "[here](http://www.fmrib.ox.ac.uk/~mark/files/av.nii.gz) if you don't have one\n",
    "handy), calculate the tSNR and then save the 3D result.\n",
    "\n",
    "> The tSNR of a time series signal is simply its mean divided by its standard\n",
    "> deviation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate tSNR"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 2
}