1 - read_viz.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9c889c80-839d-4825-88d3-46d82dd3c248",
   "metadata": {},
   "source": [
    "# Reading and Visualizing Structural MRI Data\n",
    "\n",
    "In this lesson we will learn how to read structural MRI images of a brain into Python, visualize them as slices, and perform a couple of simple image processing operations on images.\n",
    "\n",
    "Structural MRI is a generic term for any MRI scan intended to image the structure of the body. These are sometimes called anatomical scans. This is in contrast to functional MRI, which are scans designed to measure some aspect of physiological function (typically neural activation). We will work with functional MRI data later in this chapter."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e9a31aa-2012-4fbe-92fc-cb74286e9a2e",
   "metadata": {},
   "source": [
    "## Load in required packages\n",
    "\n",
    "New to us are two packages for working with images:\n",
    "- [imageio](https://pypi.org/project/imageio/) is a Python library that provides an easy interface to read and write a wide range of image data, including animated images, volumetric data, and scientific formats\n",
    "- [ndimage](https://docs.scipy.org/doc/scipy/reference/ndimage.html) is part of the SciPy pacakge, which contains functions for multidimensional image processing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa8b701b",
   "metadata": {},
   "source": [
    "~~~python\n",
    "import imageio as iio\n",
    "import scipy.ndimage as ndi\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "bbab236b-7564-4562-b9cb-d22e3303f275",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "bdc0a407-9584-4dda-8f38-e5c7b6987edf",
   "metadata": {},
   "source": [
    "## Read a DICOM image\n",
    "\n",
    "[DICOM](https://www.dicomstandard.org/) is the international standard for medical images and related information. It defines the formats for medical images that can be exchanged with the data and quality necessary for clinical use. DICOM stands for *Digital Imaging and Communications in Medicine* and is recognized by the International Organization for Standardization as the ISO 12052 standard. DICOM is the standard format for most medical imaging devices used in hospitals and clinical care settings, including MRI, PET, CT, Xray, and ultrasound. For this reason, it is common that raw MRI data comes off the scanner in DICOM format. \n",
    "\n",
    "The imageio package provides an `imread()` function that is capable of reading DICOM images. The first argument is the file name, and the second is the format. Below we read in a single slice from a structural MRI scan of a head, and show they type of the resulting object:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64d37180",
   "metadata": {},
   "source": [
    "~~~python\n",
    "brain_slice = iio.imread('data/DICOM/IM-0004-0096.dcm', 'DICOM')\n",
    "\n",
    "type(brain_slice)\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5518c2f7-4c42-4705-8746-88b21bb97d0e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "a3be39ff-b0ce-48d0-b5b3-2a6832e08e46",
   "metadata": {},
   "source": [
    "## View metadata\n",
    "We can view the image meta-data which is stored in the `.meta` property of the imageio object. This contains lots of information, such as the date and location of the scan, the manufacturer of the scanner, the type of scanner (modality; i.e., that it is an MRI scan), and details of the person who was scanned. Since this scan was acquired for research purposes, the \"patient\" information does not contain personally-identifying information like the name or birth date of the person who was scanned."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22ebfafd",
   "metadata": {},
   "source": [
    "~~~python\n",
    "brain_slice.meta\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f12e59c0-c2dc-4401-a861-8e69417fa96f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "8298e19d-cede-4d28-ac26-d494fff5f28f",
   "metadata": {},
   "source": [
    "## View image dimensions\n",
    "\n",
    "Because we loaded in a single image *slice*, our image is 2-dimensional, like a photograph. We can view its *x* and *y* dimensions with the `.shape` property. Although the image is stored as an object of type `imageio.core.util.Array`, the object is designed in a way that allows us to treat it in many ways like a simple NumPy array, so standard properties like `.shape` are available:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de06efd5",
   "metadata": {},
   "source": [
    "~~~python\n",
    "brain_slice.shape\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9caa897a-2a91-4d4f-9125-9dfb2ac9a6c1",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "baac33d1-3c27-48d2-bd87-1663735e46c2",
   "metadata": {},
   "source": [
    "The output above tells us that the image comprises 256 x 256 pixels. In medical imaging, the term **voxels** is typically used instead of pixels, because each element in the image actually represents a cube with a volume. In the metadata above, the `'SliceSpacing'` field tells us that the thickness of this slice is 1 mm, and the `'sampling'` field tells us that each voxel in the slice is 1 x 1 mm. Since these are pieces of information stored in the `.meta` prperty, we can access them individually using the following syntax:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "679fc559",
   "metadata": {},
   "source": [
    "~~~python\n",
    "brain_slice.meta['sampling']\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bab582e8-18aa-42ef-85ca-cdf32a4c1d2b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "50aa98c6-5f33-454b-a473-4cc55ccbdb54",
   "metadata": {},
   "source": [
    "## View contents of image\n",
    "\n",
    "We can use NumPy indexing to view the data values of the first row of the image:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e292a1b3",
   "metadata": {},
   "source": [
    "~~~python\n",
    "brain_slice[0]\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11a43218-09a1-4827-bbe7-3745c4897b3b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "f4c2b18b-e758-47bc-93ac-1efb8ccce36e",
   "metadata": {},
   "source": [
    "Or the first column of the image:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "814380c5",
   "metadata": {},
   "source": [
    "~~~python\n",
    "brain_slice[:, 0]\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b29305ea-b71d-404f-81dc-ce9ff39f7b17",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "546285b2-7bfd-40e6-b280-2b7eace7226d",
   "metadata": {},
   "source": [
    "## Visualize the image\n",
    "\n",
    "Since the MRI data is a 2D NumPy array of integers, we can plot it using Matplotlib's `plt.imshow()`:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13e31a11",
   "metadata": {},
   "source": [
    "~~~pythong\n",
    "plt.imshow(brain_slice)\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c690474-5a48-4ce9-9483-c06dc53dc917",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "78180bf8-9452-4448-9b4d-f0a3d5e6367c",
   "metadata": {},
   "source": [
    "We can see that this image was acquired in the **sagittal** plane, which is like looking at the person from the side.\n",
    "\n",
    "We can clean the image up by using a greyscale colour map, and removing the numerical axis labels:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0da7af50",
   "metadata": {},
   "source": [
    "~~~python\n",
    "plt.imshow(brain_slice, cmap='gray')\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c49b897-0bea-480b-a1dd-800d1e484397",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "e53306b5-cc20-44e2-a0c1-8c53a946ce07",
   "metadata": {},
   "source": [
    "## Load a brain volume\n",
    "\n",
    "Above we loaded a single slice through the head, but in fact this was from a scan that covered the entire head. As noted, each slice was 1 mm thick, and it took 184 slices to cover the entire head in the sagittal plane. Each slice is saved in a separate DICOM file. All of the DICOM images (slices) from this scan are stored in the `data/DICOM` folder. Above we used `iio.imread()` and specified a specific filename as the first argument, to load a single 2D slice/image. To load the entire 3D brain volume, we use imageio's `volread()` function, and pass as the first argument the name of the folder, rather than a list of files:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2aca403",
   "metadata": {},
   "source": [
    "~~~python\n",
    "brain_vol = iio.volread('data/DICOM', 'DICOM')\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "befff9ed-42ea-4e18-a507-68c7dac8b705",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "9c7cced8-02e9-4591-8d12-74352a720472",
   "metadata": {},
   "source": [
    "Now we can see from the `.shape` property that this is a 3D volume, with slices as the first dimension in the image array:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b48ec5ee",
   "metadata": {},
   "source": [
    "~~~python\n",
    "brain_vol.shape\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1842a04-db2e-4be2-b55f-bd931fb4ab9b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "8fe82638-9ef3-4250-a94d-85b439b2c712",
   "metadata": {},
   "source": [
    "## Visualize one slice of the volume\n",
    "\n",
    "We can plot a slice of the 3D volume using the same `plt.imshow()` command as before; the one difference is we need to specify which slice number to plot. Since slices are the first dimension of the array, we only need to supply one index, even though it is a 3D array (the other dimensions are treated as if we specified `brain_vol[96, :, :]`). We pick a slice in the middle of the volume, because if we picked one near the edges (e.g., slice 0) we would likely see little or no interesting anatomy. This time we use the `bone` colourmap, which is another monochrome palette but with a slight blueish hue:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "360edab9",
   "metadata": {},
   "source": [
    "~~~python\n",
    "plt.imshow(brain_vol[96], cmap='bone')\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a917976-1b68-400d-8885-c5cb3053f618",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "6e3212d3-4366-40a3-836b-90145915092f",
   "metadata": {},
   "source": [
    "## Visualize a slice through different planes\n",
    "\n",
    "Since the data are a 3D NumPy array, it is very easy to \"reslice\" the image and visualize the head from one of the other two orientations, axial and coronal. Below we pass `:` to select all sagittal slices, and `128` for the second dimension to get the slice midway through the volume, in the axial plane. Again, we omit the third dimension and so `:` is assumed:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7a8a18f",
   "metadata": {},
   "source": [
    "~~~python\n",
    "plt.imshow(brain_vol[:, 128], cmap='bone')\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3eb4215-d16c-44a2-bdd0-4981b3de81c8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "80680187-920e-4bbd-8dc1-d54500690e67",
   "metadata": {},
   "source": [
    "### Rotate the image\n",
    "Typically in showing axial slices, we orient them so that the nose and eyes are at the top of the image. The `scipy.ndimage` package (which we imported with the alias `ndi` provides a tool to rotate images, and we can embed it inside our `.imshow()` command to apply a rotation. The first argument to `ndi.rotate()` is the image, and the second is the amount of rotation, in degrees. The rotation is counter-clockwise, so here we need to use 270 deg rotation:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aaa0ba9f",
   "metadata": {},
   "source": [
    "~~~python\n",
    "plt.imshow(ndi.rotate(brain_vol[:, 128], 270), cmap='bone')\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74e32663-c544-43ce-a44a-56906115344b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "c1409caa-b468-4c18-88fd-9e57a68f4899",
   "metadata": {},
   "source": [
    "### Coronal plane\n",
    "Finally, we can do this in the coronal place as well; again rotation is necessary. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e898d9d",
   "metadata": {},
   "source": [
    "~~~python\n",
    "plt.imshow(ndi.rotate(brain_vol[:, :, 128], 270), cmap='bone')\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd004512-54de-4a9c-84f6-091cfd72ece4",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "9ea79f3f-dbac-49c5-afa1-7c539c0811c9",
   "metadata": {},
   "source": [
    "## Plotting a series of slices through a volume\n",
    "\n",
    "We can use our skills with Matplotlib subplots to plot a series of slices through the brain, which is a more comprehensive way of visualizing the data. The biggest trick with this is deciding on the number of subplots (slices) we want, and then doing the necessary math to select the appropriate slices from the 3D volume such that the slices are evenly-spaced through the volume, and centered in the middle of the volume. For instance, below we will generate a 4 x 4 array of 16 subplots. Our number of slices — 184 — does not divide evenly by 16 (`184 / 16 = 11.5`). For this reason, we can't simply run a `for` loop over a range of slice numbers that starts at 0 and goes up to the number of slices, in steps of `n_slices / n_subplots`. Instead, we use floor division (`//`) to generate the integer result of dividing the number of slices by subplots (`11`) so that we get a step size that ensures we have 16 evenly-spaced slices through the volume. Then, we determine the number of slices that will be covered by 16 subplots, spaced 11 slices apart from each other (which will be < 184). Finally, we compute a `start_stop` value which tells us which slice to start from (i.e., the first slice to plot), such that the set of slices we generate will be centered on the volume. This is determined by computing how many slices in the volume are *not* covered by our 16 subplots (i.e., the slices at either edge of the volume), and dividing this by 2, since half of those slices should be at one edge of the volume, at the other at the other end of the volume."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f683e049",
   "metadata": {},
   "source": [
    "~~~python\n",
    "fig_rows = 4\n",
    "fig_cols = 4\n",
    "n_subplots = fig_rows * fig_cols\n",
    "n_slice = brain_vol.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",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "5ff3027b-6067-4b5b-b2b2-1d3c44a79378",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "d16a4147-28f1-4101-85fb-44782c6e1aca",
   "metadata": {},
   "source": [
    "Having done this math, we now loop through the slices and plot them, starting with the slice defined by `start_stop`:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdf3b64f",
   "metadata": {},
   "source": [
    "~~~python\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(brain_vol[img, :, :], cmap='gray')\n",
    "    axs.flat[idx].axis('off')\n",
    "        \n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56c7e19e-50f5-4e12-a1b2-4cd7d05549f0",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "783739c9-c3cb-49a4-91f9-2b51cc9975c2",
   "metadata": {},
   "source": [
    "## Slice through a different axis\n",
    "\n",
    "We can use the same approach to plot the data through other image planes. The only things we need to change are: \n",
    "- which dimension of the image to use to derive `n_slice`\n",
    "- the dimension that we specify the slice number in, inside the `.imshow()` command\n",
    "- adding rotation for axial and coronal slices, as we did above when plotting a single slice"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41547b63",
   "metadata": {},
   "source": [
    "~~~python\n",
    "n_slice = brain_vol.shape[1]\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[:, img, :], 270), \n",
    "                         cmap='gray')\n",
    "    axs.flat[idx].axis('off')\n",
    "        \n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7556e54f-66a4-4d16-bcaa-cc79d1c437e0",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "271a64e4-b0fc-4f34-841e-9675009e2f9d",
   "metadata": {},
   "source": [
    "### Coronal plane"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "433efdea",
   "metadata": {},
   "source": [
    "~~~python\n",
    "fig_rows = 4\n",
    "fig_cols = 4\n",
    "n_subplots = fig_rows * fig_cols\n",
    "n_slice = brain_vol.shape[2]\n",
    "step_size = n_slice // n_subplots\n",
    "plot_range = n_subplots * step_size\n",
    "\n",
    "start_stop = int((n_slice - plot_range) / 2)\n",
    "\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[:, :, img], 270), \n",
    "                         cmap='gray')\n",
    "    axs.flat[idx].axis('off')\n",
    "        \n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77eedb2c-e2dd-4cbb-8128-22deb82f4e44",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "88940a11-aa08-48b0-9521-f1b80f580861",
   "metadata": {},
   "source": [
    "## Plotting an image histogram\n",
    "\n",
    "As we saw above, the image is stored as a NumPy array, in which each voxel in the image is represented as a number, which is mapped to an intensity value in the colourmap when plotting. Larger values appear as brighter (whiter), and lower values appear as darker.\n",
    "\n",
    "Histograms of the anatomical images show the number of voxels of a given intensity value. These can be informative because the distribution of intensity values in an anatomical image is not uniform. Instead, as we can see above, there are many very dark voxels (outside of the head, and in some of the fluid-filled spaces inside the head), and then clusters of voxels that are darker grey (the grey matter, largely in the cerebral cortex that forms the outer layer of the brain), lighter grey (the white matter that comprises much of the inside of the brain), and also some very bright areas that are primarily due to areas of fat concentration.\n",
    "\n",
    "We can use ndimage's `.histogram()` function to plot a histogram of our brain volume. We use this rather than the NumPy `histogram` function, because ndimage's function is designed to work with 3D images. This function requires several arguments, including the minimum and maximum intensity values that define the range of the *x* axis of the histogram, and the number of bins:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "963436ea",
   "metadata": {},
   "source": [
    "~~~python\n",
    "plt.plot(ndi.histogram(brain_vol, min=0, max=np.max(brain_vol), bins=50))\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7511b0df-f6c1-4cfc-9785-34b96f868413",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "329313fe-ef97-4f08-86ee-fa5df0c681a8",
   "metadata": {},
   "source": [
    "In the histogram above, there is a large peak close to zero which represents the fact that a large number of voxels in the image don't contain the head at all, and therefore have values at or close to zero. We can see a peak just above 10 on the *x* axis (note that the numbers on the *x* axis are bin numbers, not intensity values), with a slight decrease and then a second small peak just before 20, followed by a flat area. The peaks just above 10 and just below 20 reflect the concentration of similar intensity values corresponding to grey and white matter respectively. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "244f589f-a3d9-4702-ae39-f992b18a2b6c",
   "metadata": {},
   "source": [
    "## Mask an image\n",
    "\n",
    "We can use an image histogram like this to create a **mask** that isolates a particular range of intensity values in the image, while setting all other intensity values to zero. This can be useful in certain types of analysis, such as if we want to isolate only grey matter, or only white matter, from the rest of the brain. This is often used in studies that measure the volume of grey and/or white matter separately, as well as in functional MRI analyses where the interest is often only in activation in the cerebral cortex (grey matter). \n",
    "\n",
    "In the code below, we will attempt to isolate grey matter based on intensity values from the histogram above. Note that this is an overly-simplistic approach to separating grey and white matter in MRIs, compared to more sophisticated, automated, and accurate approaches that would be used in research or clinical care. However, it is an important step in those more sophisticated approaches.\n",
    "\n",
    "First, we manually define the range of intensity values that we consider grey matter. Based on the histogram above, we will choose a range from bins 10 – 15. Since the *x* axis of the histogram is bin numbers, not intensity values, we need to determine what intensities those bin numbers correspond to. Since we know the histogram divided the range from 0 to the maximum intensity value in the image, into 50 bins, we can divide the max intensity value by 50 to get the width of each bin, then multiply by the values we observed on the *x* axis:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbe8c0cb",
   "metadata": {},
   "source": [
    "~~~python\n",
    "gm_min = ((np.max(brain_vol)) / 50) * 10\n",
    "gm_max = ((np.max(brain_vol)) / 50) * 15\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "c237b28f-bc75-424d-a23d-358acc2217f0",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "d6db519a-a798-404a-820f-9856b10b39c0",
   "metadata": {},
   "source": [
    "Next we use `np.where()` to create two *binary masks* of the image. A binary mask is an image (in this case, a NumPy array) in which each voxel's value is either 1 or 0. For the first mask, any voxel whose intensity is greater than `gm_min` is set to 1, and the others (lower values) are set to 0. In the second mask, any voxel less than `gm_max` is set to 1, and larger values are set to zero."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5c73d24",
   "metadata": {},
   "source": [
    "~~~python\n",
    "brain_mask1 = np.where(brain_vol > gm_min, 1, 0)\n",
    "brain_mask2 = np.where(brain_vol < gm_max, 1, 0)\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "5ed605bc-5d87-4943-8ac6-1ff0a3598248",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "a427ee8f-35f0-4f24-a8b3-8005789b144c",
   "metadata": {},
   "source": [
    "Having done this, we combine the two masks by adding them. Now, any voxel that is in the range between `gm_min` and `gm_max` will have a value of 2. Finally, we use `np.where()` again to create a final mask in which any voxel from the combined mask with a value of 2 (grey matter) is set to 1, and all other voxels (those both above and below our grey matter range) are set to zero."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2081743b",
   "metadata": {},
   "source": [
    "~~~python\n",
    "brain_mask = brain_mask1 + brain_mask2\n",
    "brain_mask = np.where(brain_mask == 2, 1, 0)\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "9f256ee4-1785-4866-af2c-4f74778770d5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "63ad888e-9e49-4514-9079-b88db9d0fecc",
   "metadata": {},
   "source": [
    "When we plot this final mask, we can see that it does a pretty good job of isolating the grey matter of the cortex. There is also grey matter around the ventricles (those dark areas in the centre of the brain in the images we plotted above) which shows up here. Unfortunately, a lot of the tissue in the scalp is also in the same intensity range as grey matter, so we are not able to separate these based on intensity values. Often prior to performing this type of image segmentation, a preprocessing step called *skull stripping* is applied, which (as it sounds) isolates and removes the skull and other non-brain tissues from the image. Skull stripping is more advanced than we will cover here. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "397955fc",
   "metadata": {},
   "source": [
    "~~~python\n",
    "plt.imshow(brain_mask[:, 96, :], cmap='gray')\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc782600-987e-4542-9cef-4cbeb44d8d61",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "4984e582-3d5b-48c9-b309-5f1915c2cc01",
   "metadata": {},
   "source": [
    "## Image smoothing\n",
    "\n",
    "One preprocessing operation that can be useful in working with images is **smoothing**. Smoothing is really spatial filtering, and it is commonly applied by averaging together the intensity values of nearby voxels. This averaging is done in a weighted fashion, based on distance between voxels. For example, for each voxel we could average together the intensity of that voxel, with half the intensity of the voxels adjacent to it on all sides, and perhaps 10% of intensities that were 2 voxels away from it. In fact, the most common **smoothing kernel** (the mathematical function used to perform this weighted averaging) is called a *Gaussian* kernel, which is shaped like a normal (bell) curve. In 2 dimensions, it looks like this:\n",
    "\n",
    "![](./images/gaussian_kernel_2D.jpg)\n",
    "\n",
    "\n",
    "When applied to an image slice, a Gaussian smoothing kernel reduces noise, that is, the amount of voxel-to-voxel variation in intensities. So, adjacent voxels will have more similar intensity values after smoothing, making the image appear more blurred. The image below shows an example of the effects of smoothing. In this image, the intensity at each point in the image slice is shown by the height of the mesh at that location, as well as the colour. \n",
    "\n",
    "![](./images/smoothed_vs_unsmoothed_mesh.jpg)\n",
    "\n",
    "Below we will apply Gaussian smoothing to our MRI volume using ndimage's `gaussian_filter()` function. The one argument we need to supply is `sigma`, which is the width of the smoothing kernel, expressed as units of standard deviation. A larger sigma value will result in a smoother (blurrier) image, because we average over a larger number of voxels.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b69673b",
   "metadata": {},
   "source": [
    "~~~python\n",
    "sigma = 2\n",
    "smoothed = ndi.gaussian_filter(brain_vol, sigma)\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "id": "2fe3c036-98c3-49f0-be2b-cb551aae0228",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "37cf431c-80b1-475c-9b35-1605db84613c",
   "metadata": {},
   "source": [
    "To plot the image we just re-use the code we used above:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b86592ba",
   "metadata": {},
   "source": [
    "~~~python\n",
    "fig_rows = 4\n",
    "fig_cols = 4\n",
    "n_subplots = fig_rows * fig_cols\n",
    "n_slice = brain_vol.shape[0]\n",
    "step_size = n_slice // n_subplots\n",
    "plot_range = n_subplots * step_size\n",
    "\n",
    "start_stop = int((n_slice - plot_range) / 2)\n",
    "\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(smoothed[img, :, :], cmap='gray')\n",
    "    axs.flat[idx].axis('off')\n",
    "        \n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a150354-c5c6-40a7-bd77-d0ff18ef64ed",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "81d3eef4-6948-402a-acc6-3770a20b61ca",
   "metadata": {},
   "source": [
    "Below we increase `sigma` to 4, resulting in a more smoothed image:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbb14bc3",
   "metadata": {},
   "source": [
    "~~~python\n",
    "sigma = 4\n",
    "smoothed = ndi.gaussian_filter(brain_vol, sigma)\n",
    "\n",
    "fig_rows = 4\n",
    "fig_cols = 4\n",
    "n_subplots = fig_rows * fig_cols\n",
    "n_slice = brain_vol.shape[0]\n",
    "step_size = n_slice // n_subplots\n",
    "plot_range = n_subplots * step_size\n",
    "\n",
    "start_stop = int((n_slice - plot_range) / 2)\n",
    "\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(smoothed[img, :, :], cmap='gray')\n",
    "    axs.flat[idx].axis('off')\n",
    "        \n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2e79394-a352-4140-85d6-00c4cae559b8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "1d89ed14-7f09-4de3-b398-aafd90c4a6b7",
   "metadata": {},
   "source": [
    "## Segmenting smoothed images\n",
    "\n",
    "One benefit of smoothing, and the corresponding reduction in noise in the image, is that it can make segmentation a bit cleaner, because adjacent voxels will have more similar values — so they are less extreme overall, and less variable from voxel to voxel. Below we plot the histogram of the smoothed image, and use it (as we did earlier) to determine cutoffs for grey matter and then segment the image. Note that the histogram is different, and more smooth, due to the smoothing applied to the intensity values."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b555fc9",
   "metadata": {},
   "source": [
    "~~~python\n",
    "filt = ndi.gaussian_filter(brain_vol, sigma=2)\n",
    "plt.plot(ndi.histogram(filt, min=0, max=np.max(filt), bins=50))\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0bf30267-bce9-4714-97cf-0f21524fe201",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "370706e7-916d-4a40-8c67-795bf8a49706",
   "metadata": {},
   "source": [
    "Based on the above histogram we select 16 – 25 as our range of bins for grey matter."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3371425",
   "metadata": {},
   "source": [
    "~~~python\n",
    "gm_min = ((np.max(filt)) / 50) * 16\n",
    "gm_max = ((np.max(filt)) / 50) * 25\n",
    "\n",
    "brain_mask1 = np.where(filt > gm_min, 1, 0)\n",
    "brain_mask2 = np.where(filt < gm_max, 1, 0)\n",
    "\n",
    "brain_mask = brain_mask1 + brain_mask2\n",
    "brain_mask = np.where(brain_mask == 2, 1, 0)\n",
    "\n",
    "plt.imshow(brain_mask[:, 96, :], cmap='gray')\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eea15f3c-3c1e-4b77-8528-f65ffca15078",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "1a54f729-5d55-44ac-83ad-adb526a153ec",
   "metadata": {},
   "source": [
    "We can see that this results in a 'cleaner' depiction of the grey matter, and in general the grey matter is separated by black (0) voxels from the skull."
   ]
  }
 ],
 "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
}