{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Tissue Classification of a T1-weighted Structural Image\n\nThis example explains how to segment a T1-weighted structural image by using\nBayesian formulation. The observation model (likelihood term) is defined as a\nGaussian distribution and a Markov Random Field (MRF) is used to model the\na priori probability of context-dependent patterns of different tissue\ntypes of the brain. Expectation Maximization and Iterated Conditional\nModes are used to find the optimal solution. Similar algorithms have been\nproposed by Zhang et al. [Zhang2001]_ and Avants et al. [Avants2011]_ available\nin FAST-FSL and ANTS-atropos, respectively.\n\nHere we will use a T1-weighted image, that has been previously skull-stripped\nand bias field corrected.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport matplotlib.pyplot as plt\nfrom dipy.data import get_fnames\nfrom dipy.io.image import load_nifti_data\nfrom dipy.segment.tissue import TissueClassifierHMRF\nimport time"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First we fetch the T1 volume from the Syn dataset and determine its shape.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "t1_fname, _, _ = get_fnames('tissue_data')\nt1 = load_nifti_data(t1_fname)\nprint('t1.shape (%d, %d, %d)' % t1.shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We have fetched the T1 volume. Now we will look at the axial and coronal\nslices of the image.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure()\na = fig.add_subplot(1, 2, 1)\nimg_ax = np.rot90(t1[..., 89])\nimgplot = plt.imshow(img_ax, cmap=\"gray\")\na.axis('off')\na.set_title('Axial')\na = fig.add_subplot(1, 2, 2)\nimg_cor = np.rot90(t1[:, 128, :])\nimgplot = plt.imshow(img_cor, cmap=\"gray\")\na.axis('off')\na.set_title('Coronal')\nplt.savefig('t1_image.png', bbox_inches='tight', pad_inches=0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: t1_image.png\n   :align: center\n\n   T1-weighted image of healthy adult.\n\nNow we will define the other two parameters for the segmentation algorithm.\nWe will segment three classes, namely corticospinal fluid (CSF), white matter\n(WM) and gray matter (GM).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nclass = 3"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Then, the smoothness factor of the segmentation. Good performance is achieved\nwith values between 0 and 0.5.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "beta = 0.1"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We could also set the number of iterations. By default this parameter is set to\n100 iterations, but most of the time the ICM (Iterated Conditional Modes)\nloop will converge before reaching the 100th iteration.\nAfter setting the necessary parameters we can now call an instance of the class\n\"TissueClassifierHMRF\" and its method called \"classify\" and input the\nparameters defined above to perform the segmentation task.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "t0 = time.time()\n\nhmrf = TissueClassifierHMRF()\ninitial_segmentation, final_segmentation, PVE = hmrf.classify(t1, nclass, beta)\n\nt1 = time.time()\ntotal_time = t1-t0\nprint('Total time:' + str(total_time))\n\nfig = plt.figure()\na = fig.add_subplot(1, 2, 1)\nimg_ax = np.rot90(final_segmentation[..., 89])\nimgplot = plt.imshow(img_ax)\na.axis('off')\na.set_title('Axial')\na = fig.add_subplot(1, 2, 2)\nimg_cor = np.rot90(final_segmentation[:, 128, :])\nimgplot = plt.imshow(img_cor)\na.axis('off')\na.set_title('Coronal')\nplt.savefig('final_seg.png', bbox_inches='tight', pad_inches=0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we plot the resulting segmentation.\n\n.. figure:: final_seg.png\n   :align: center\n\n   Each tissue class is color coded separately, red for the WM, yellow for\n   the GM and light blue for the CSF.\n\nAnd we will also have a look at the probability maps for each tissue class.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure()\na = fig.add_subplot(1, 3, 1)\nimg_ax = np.rot90(PVE[..., 89, 0])\nimgplot = plt.imshow(img_ax, cmap=\"gray\")\na.axis('off')\na.set_title('CSF')\na = fig.add_subplot(1, 3, 2)\nimg_cor = np.rot90(PVE[:, :, 89, 1])\nimgplot = plt.imshow(img_cor, cmap=\"gray\")\na.axis('off')\na.set_title('Gray Matter')\na = fig.add_subplot(1, 3, 3)\nimg_cor = np.rot90(PVE[:, :, 89, 2])\nimgplot = plt.imshow(img_cor, cmap=\"gray\")\na.axis('off')\na.set_title('White Matter')\nplt.savefig('probabilities.png', bbox_inches='tight', pad_inches=0)\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: probabilities.png\n   :align: center\n   :scale: 120\n\n   These are the probability maps of each of the three tissue classes.\n\n.. [Zhang2001] Zhang, Y., Brady, M. and Smith, S. Segmentation of Brain MR\n   Images Through a Hidden Markov Random Field Model and the\n   Expectation-Maximization Algorithm IEEE Transactions on Medical Imaging,\n   20(1): 45-56, 2001\n\n.. [Avants2011] Avants, B. B., Tustison, N. J., Wu, J., Cook, P. A. and Gee,\n   J. C. An open source multivariate framework for n-tissue segmentation with\n   evaluation on public data. Neuroinformatics, 9(4): 381-400, 2011.\n\n\n"
      ]
    }
  ],
  "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.9.13"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}