image_processing/image_processing_FSL.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "00c84ee5",
      "metadata": {},
      "source": [
        "# Image Processing using FSL\n",
        "<hr/>\n",
        "***As-Is Software Disclaimer***\n",
        "\n",
        "This content in this repository is delivered “As-Is”. Notwithstanding anything to the contrary, DNAnexus will have no warranty, support, liability or other obligations with respect to Materials provided hereunder.\n",
        "\n",
        "<hr/>\n",
        "\n",
        "This notebook demonstrates usage of the Image Processing package <a href=\"https://fsl.fmrib.ox.ac.uk/fsl/fslwiki\">FSL</a> on the DNAnexus platform, to perform Brain Extraction and Isotropic Smooth as individual analyses and also in the form of a NiPype workflow. \n",
        "\n",
        "<a href=\"https://github.com/dnanexus/OpenBio/blob/master/LICENSE.md\">MIT License</a> applies to this notebook."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "f6b94e08",
      "metadata": {},
      "source": [
        "## 1. Preparing your environment\n",
        "### Launch spec:\n",
        "\n",
        "* App name: JupyterLab with Python, R, Stata, ML, Image Processing\n",
        "* App flavor: IMAGE_PROCESSING\n",
        "* Kernel: Python\n",
        "* Instance type: mem1_ssd1_v2_x16\n",
        "* Runtime: =~ 5 min\n",
        "* Data description: Input for this notebook is a T1-weighted MRI image. We are downloading the primer example data from <a href=\"https://www.fmrib.ox.ac.uk/primers/intro_primer/ExBox13/IntroBox13.html\">Oxford NeuroImaging Primer</a>. The downloaded zipped folder (size ~30 MB) contains the following images: \n",
        "    * Original (non-brain-extracted) structural image: `T1.nii.gz`\n",
        "    * Two brain extracted images: `T1_v1_brain.nii.gz` and `T1_v2_brain.nii.gz`\n",
        "    \n",
        "### Package and tools dependency:\n",
        "\n",
        "| Package | License | \n",
        "| --- | --- |\n",
        "| <a href=\"https://graphviz.org/license/\">graphviz</a> | <a href=\"https://opensource.org/licenses/CPL-1.0\">Common Public License</a> |"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "bb5775d8",
      "metadata": {},
      "source": [
        "## 2. Download and Unzip an MRI"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "1f0c414c",
      "metadata": {},
      "outputs": [],
      "source": [
        "! wget https://www.fmrib.ox.ac.uk/primers/intro_primer/ExBox13/ExBox13.zip"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "78ed1a2d",
      "metadata": {},
      "outputs": [],
      "source": [
        "! apt-get update\n",
        "! apt-get install unzip\n",
        "! unzip ExBox13.zip\n",
        "! ls ExBox13/"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "18a3b57a",
      "metadata": {},
      "source": [
        "## 3. Install and Import Packages"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "d8263228",
      "metadata": {},
      "outputs": [],
      "source": [
        "! yes | apt-get install graphviz"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "e676434b",
      "metadata": {},
      "outputs": [],
      "source": [
        "from nipype.interfaces import fsl\n",
        "from nipype import Node, Workflow\n",
        "import nibabel as nib\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import pandas as pd\n",
        "import os\n",
        "\n",
        "plt.style.use(\"seaborn-dark\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "fb502e10",
      "metadata": {},
      "source": [
        "## 4. Load MRI and Get Metadata"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a11c7b81",
      "metadata": {},
      "source": [
        "### Declare input file"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "58bb7042",
      "metadata": {},
      "outputs": [],
      "source": [
        "mri_dir = \"/opt/notebooks/ExBox13/\"\n",
        "t1_filename = \"T1.nii.gz\"\n",
        "t1_in_file = os.path.join(mri_dir, t1_filename)\n",
        "t1_in_file"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "346864ac",
      "metadata": {},
      "source": [
        "The `nib.load()` function returns an nibabel object and allows us to get metadata about the MRI image without loading the actual image into memory."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "635f0063",
      "metadata": {},
      "outputs": [],
      "source": [
        "t1 = nib.load(t1_in_file)\n",
        "print(\"MRI shape:\", t1.shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a378b629",
      "metadata": {},
      "source": [
        "The header attribute contains important metadata for an MRI (e.g. image dimensions, data type, etc.)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "32c37ce7",
      "metadata": {},
      "outputs": [],
      "source": [
        "print(t1.header)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d222c170",
      "metadata": {},
      "source": [
        "## 5. Visualize MRI Slices"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "f541c602",
      "metadata": {},
      "source": [
        "Set the default output type of images as `*.nifti.gz`"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "6ac78da0",
      "metadata": {},
      "outputs": [],
      "source": [
        "fsl.FSLCommand.set_default_output_type(\"NIFTI_GZ\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6216a3c5",
      "metadata": {},
      "source": [
        "The nibabel object generated above does not load the actual image data. We use the `get_fdata()` function to load the image array data."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "a146b36f",
      "metadata": {},
      "outputs": [],
      "source": [
        "t1_img_data = t1.get_fdata()\n",
        "\n",
        "fig, ax = plt.subplots(1, 3, figsize=(10, 4))\n",
        "fig.suptitle(\"T1.nii.gz\", fontsize=14)\n",
        "\n",
        "ax[0].imshow(t1_img_data[100, :, :], cmap=\"gray\")\n",
        "ax[0].set_title(\"Sagittal Slice\")\n",
        "\n",
        "ax[1].imshow(t1_img_data[:, 100, :], cmap=\"gray\")\n",
        "ax[1].set_title(\"Frontal Slice\")\n",
        "\n",
        "ax[2].imshow(t1_img_data[:, :, 100], cmap=\"gray\")\n",
        "ax[2].set_title(\"Axial Slice\")\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "7ddfdea9",
      "metadata": {},
      "source": [
        "## 6. Extract Brain image from the whole head\n",
        "We use the `BET` (Brain Extraction Tool) function to extract only the brain tissue and delete non-brain parts from an image of the whole head."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "19e81f94",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Skullstrip process\n",
        "skullstrip = fsl.BET(\n",
        "    in_file=t1_in_file, out_file=mri_dir + \"/T1_skullstrip.nii.gz\", mask=True\n",
        ")\n",
        "skullstrip.run()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "6c25454a",
      "metadata": {},
      "outputs": [],
      "source": [
        "skullstrip_img_data = nib.load(\n",
        "    os.path.join(mri_dir, \"T1_skullstrip.nii.gz\")\n",
        ").get_fdata()\n",
        "\n",
        "fig, ax = plt.subplots(1, 3, figsize=(10, 4))\n",
        "fig.suptitle(\"T1_skullstrip.nii.gz\", fontsize=14)\n",
        "\n",
        "ax[0].imshow(skullstrip_img_data[100, :, :], cmap=\"gray\")\n",
        "ax[0].set_title(\"Sagittal Slice\")\n",
        "\n",
        "ax[1].imshow(skullstrip_img_data[:, 100, :], cmap=\"gray\")\n",
        "ax[1].set_title(\"Frontal Slice\")\n",
        "\n",
        "ax[2].imshow(skullstrip_img_data[:, :, 100], cmap=\"gray\")\n",
        "ax[2].set_title(\"Axial Slice\")\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ca7a93c0",
      "metadata": {},
      "source": [
        "## 7. Apply Isotropic Smooth\n",
        "We use the `IsotropicSmooth()` function to spatially smooth an image by applying a Gaussian filter."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "64699ad6",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Smoothing process\n",
        "smooth = fsl.IsotropicSmooth(\n",
        "    in_file=t1_in_file, out_file=mri_dir + \"/T1_smooth.nii.gz\", fwhm=4\n",
        ")\n",
        "smooth.run()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "e09c0500",
      "metadata": {},
      "outputs": [],
      "source": [
        "smooth_img_data = nib.load(os.path.join(mri_dir, \"T1_smooth.nii.gz\")).get_fdata()\n",
        "\n",
        "fig, ax = plt.subplots(1, 3, figsize=(10, 4))\n",
        "fig.suptitle(\"T1_smooth.nii.gz\", fontsize=14)\n",
        "ax[0].imshow(smooth_img_data[100, :, :], cmap=\"gray\")\n",
        "ax[0].set_title(\"Sagittal Slice\")\n",
        "\n",
        "ax[1].imshow(smooth_img_data[:, 100, :], cmap=\"gray\")\n",
        "ax[1].set_title(\"Frontal Slice\")\n",
        "\n",
        "ax[2].imshow(smooth_img_data[:, :, 100], cmap=\"gray\")\n",
        "ax[2].set_title(\"Axial Slice\")\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6c37e34f",
      "metadata": {},
      "source": [
        "## 8. Running a Workflow\n",
        "So far, we ran analysis on the MRI in the form of individual interfaces. If we want to run multiple interfaces in a specific order, we can build and run a workflow. The workflow is a Directed Acyclic Graph consisting of individual interfaces as nodes and connections between the interfaces as edges. \n",
        "\n",
        "For more information information on workflow, visit: https://github.com/miykael/nipype_tutorial/blob/master/notebooks/basic_workflow.ipynb\n",
        "\n",
        "### Create nodes"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "f82354f9",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Skullstrip process\n",
        "skullstrip = Node(fsl.BET(in_file=t1_in_file, mask=True), name=\"skullstrip\")\n",
        "\n",
        "# Smooth process\n",
        "smooth = Node(fsl.IsotropicSmooth(in_file=t1_in_file, fwhm=4), name=\"smooth\")\n",
        "\n",
        "# Mask process\n",
        "mask = Node(fsl.ApplyMask(), name=\"mask\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "26b7e231",
      "metadata": {},
      "source": [
        "### Initiate workflow"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "36617ab4",
      "metadata": {},
      "outputs": [],
      "source": [
        "wf = Workflow(name=\"smoothflow\", base_dir=mri_dir)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "13d84327",
      "metadata": {},
      "source": [
        "### Connect nodes to form edges"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "13d037f2",
      "metadata": {},
      "outputs": [],
      "source": [
        "wf.connect(\n",
        "    [\n",
        "        (skullstrip, mask, [(\"mask_file\", \"mask_file\")]),\n",
        "        (smooth, mask, [(\"out_file\", \"in_file\")]),\n",
        "    ]\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b880f27d",
      "metadata": {},
      "source": [
        "### Visualize the workflow"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "e3da4f02",
      "metadata": {},
      "outputs": [],
      "source": [
        "wf.write_graph(graph2use=\"flat\")\n",
        "from IPython.display import Image\n",
        "\n",
        "Image(filename=wf.base_dir + \"/smoothflow/graph_detailed.png\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6554bf40",
      "metadata": {},
      "source": [
        "### Run the workflow"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "dd2c6003",
      "metadata": {},
      "outputs": [],
      "source": [
        "wf.run()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2b6de8d1",
      "metadata": {},
      "source": [
        "### Visualize the results"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "39b4ed5b",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Helper function to plot images\n",
        "def plot_slice(file_name):\n",
        "\n",
        "    # Load the image\n",
        "    img = nib.load(file_name)\n",
        "    data = img.get_fdata()\n",
        "\n",
        "    # Cut in the middle of the brain\n",
        "    cut = int(data.shape[-1] / 2) + 10\n",
        "\n",
        "    # Plot the data\n",
        "    plt.title(os.path.basename(file_name))\n",
        "    plt.imshow(data[..., cut], cmap=\"gray\")\n",
        "    plt.gca().set_axis_off()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "9868a101",
      "metadata": {},
      "outputs": [],
      "source": [
        "fig = plt.figure(figsize=(12, 4))\n",
        "for i, img in enumerate(\n",
        "    [\n",
        "        t1_in_file,\n",
        "        wf.base_dir + \"/smoothflow/smooth/T1_smooth.nii.gz\",\n",
        "        wf.base_dir + \"/smoothflow/skullstrip/T1_brain_mask.nii.gz\",\n",
        "        wf.base_dir + \"/smoothflow/mask/T1_smooth_masked.nii.gz\",\n",
        "    ]\n",
        "):\n",
        "    fig.add_subplot(1, 4, i + 1)\n",
        "    plot_slice(img)"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "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.6.5"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}