Tutorial_SimpleElastix.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## These two packages were not included in the requirement sheet.\n",
    "## If you don't have them installed, please run this cell or directly in your command line\n",
    "\n",
    "!pip install matplotlib\n",
    "!pip install ipywidgets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Import Jupyter Notebook specifics\n",
    "%matplotlib notebook\n",
    "from ipywidgets import *\n",
    "\n",
    "## Import numpy\n",
    "import numpy as np\n",
    "\n",
    "## Import SimpleITK - We can access the SimpleElastix from here\n",
    "import SimpleITK as sitk\n",
    "\n",
    "## Import supporting files\n",
    "import support_files.dataManagement as dataManage\n",
    "import support_files.imageUtilities as imageUtilities\n",
    "import support_files.t1MappingPredict as t1Map\n",
    "from support_files.imageBrowse import imageBrowse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to T1 mapping\n",
    "\n",
    "## Why do we need mapping?\n",
    "The image formation in a MRI experiment is a complex process. Underlying parameters of differently tissues respond differently to the external magnetic field produced by the machine. This creates changes in image contrast, allowing for tissue differentiation.\n",
    "\n",
    "However, as most of physical phenomena, this is a noisy procedure, since internal (i.e. variation in tissue susceptibility) and external factors (i.e. inhomogeneity of magnetic fields, imperfect RF pulses) affect the appearance of the final weighted images. This happens because the observed signal amplitudes are not only dependent on the tissue properties, therefore making longitudinal analysis much harder.\n",
    "\n",
    "## Quantitative mapping\n",
    "Quantitative mapping aims to reduce the corrupting factors by estimating the invariant tissue property directly from the observed signal. T1 mapping, the focus of this exercise, can be performed in several ways, but the most common approaches are based on statistical iterative algorithms, such as Bayesian inference or the Maximum Likelihood Estimator (MLE).\n",
    "\n",
    "The optimisation algorithm fits a given signal model to the acquired weighted images to try to recover the tissue properties that best explain the observed signal. In the cases presented here, Inversion Recovery (IR) based sequences were used and the approximated signal model is given by:\n",
    "\n",
    "$$f_{n}\\left(A, B, T_{1}\\right)=\\left|A-B e^{-\\frac{TI_{n}}{T_{1}}}\\right|$$\n",
    "\n",
    "where A is proportional to the proton density of the sample, B is related to the inversion efficiency of the RF pulses, TI<sub>n</sub> is the inversion time for the n<sup>th</sup> image in the series and T<sub>1</sub> is the longitudinal relaxation time in a voxel. \n",
    "\n",
    "By varying the time the RF pulse is applied (the TI), the contrast of the images change according to their A, B and T<sub>1</sub> values and the signal model can be used to fit the time series at a voxel. In the figure below you can see how the signal is dependent on the acquisition TI. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can observe how the signal from different tissues evolve:\n",
    "\n",
    "<img src=\"files/support_images/Synt_T1w.png\">\n",
    "\n",
    "The image above only shows 5 images out of 23 images in total. The plot below shows the signal evolution for gray matter (GM), white matter (WM) and csf. Each point in the curve corresponds to the voxel intensity at a give inversion time.\n",
    "\n",
    "<img src=\"files/support_images/T1w_signalEvolution.png\">\n",
    "\n",
    "The estimation algorithm, through iterative optimisation methods, trys to predict which values for T<sub>1</sub>, A and B that could have originated this signal evolution, based on the signal model. In the plot above, by observing the inflexion points, it is possible to infer that T<sub>1,CSF </sub> > T<sub>1,GM </sub> > T<sub>1,WM </sub>.\n",
    "\n",
    "\n",
    "\n",
    "## Visualising the weighted images\n",
    "Now you can load and visualise our first example, a set of synthetic weighted images simulated using the signal model above. Explore the variation in signal intensity and check for different signal evolution curves for multiple voxels:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Before running the cell below, please unzip the file \"synthetic_t1w.h5.zip\" in the folder \"./support_data/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Loading synthetic data\n",
    "filePath = \"./support_data/synthetic_t1w.h5\"\n",
    "originalSeries, groundTruthMaps, mask = dataManage.loadSyntheticData(filePath, sliceSelect=38)\n",
    "gt_t1Map = groundTruthMaps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## Visualisation\n",
    "imageBrowse(originalSeries[\"weighted_series\"], \n",
    "            im_ax = 0, \n",
    "            rang_min = 0.0,\n",
    "            rang_max = 1.5,\n",
    "            colormap_c = 'gray',\n",
    "            fig_title='T1 weighted series',\n",
    "            ylabel='[a.u.]'\n",
    "           )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How does motion affect the images\n",
    "\n",
    "### Patient motion\n",
    "\n",
    "In an Inversion Recovery experiment, to permit the complete relaxation of the signal, the images at each inversion times are acquired spaced from one another. Due to inherent patient movement during the scan, the images at different TI<sub>s</sub> won't be aligned. This means that throughout the time series, a given voxel in the images won't, necessarily, correspond to the same spatial location, since they are misaligned.\n",
    "\n",
    "In the example below, we introduced simulated motion to the synthetic dataset from before. In the same way as above, you can visualise the different images, and check how the simulated motion affects the signal evolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Loading motion corrupted data\n",
    "filePath = \"./support_data/synthetic_t1w_deformed.nii\"\n",
    "corruptedSeries = {}\n",
    "corruptedSeries[\"weighted_series\"] = dataManage.readNii(filePath)\n",
    "corruptedSeries[\"echo_times\"] = originalSeries[\"echo_times\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Visualisation\n",
    "imageBrowse(corruptedSeries[\"weighted_series\"], \n",
    "            im_ax = 0, \n",
    "            rang_min = 0.0, \n",
    "            rang_max = 1.5, \n",
    "            colormap_c='gray', \n",
    "            fig_title='Motion corr. T1w',\n",
    "            ylabel='[a.u.]'\n",
    "           )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Impact of motion on T1 mapping\n",
    "To demonstrate how this affects the mapping results, you can run the cell below. It evoques a MLE-based optimisation method which estimates the A, B and T<sub>1</sub> maps, however, here, we are only interested in the latter. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the T1Mapping wrapper\n",
    "def t1Mapping(data):\n",
    "    print('-'*40)\n",
    "    print('T1 estimation started')\n",
    "    \n",
    "    configurationFilePath = './support_files/configuration_file_mle.txt'\n",
    "    predictedMaps, simulatedWeightedSequence = t1Map.predictMapping(data, configurationFilePath, True)\n",
    "    \n",
    "    print('T1 estimation finished')\n",
    "    return np.array(predictedMaps), np.array(simulatedWeightedSequence)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Run T1 mapping for original data\n",
    "qMaps, estimatedWeightedSeriesOriginal = t1Mapping(originalSeries)\n",
    "original_t1Map = qMaps[-1].detach().numpy()\n",
    "\n",
    "## Run T1 mapping for motion corrupted data\n",
    "qMaps, estimatedWeightedSeriesCorrupted = t1Mapping(corruptedSeries)\n",
    "corrupted_t1Map = qMaps[-1].detach().numpy()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparing to the known ground-truth (only in simulated experiments)\n",
    "An advantage of using simulated data is that the ground-truth maps (A, B and T<sub>1</sub>) are known. To evaluate the precision of the estimation, we can directly compared the estimated maps to the known ground-truth:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## Comparison plot\n",
    "imageUtilities.compare_predictions(gt_t1Map*mask,\n",
    "                                   original_t1Map*mask, \n",
    "                                   colormap_c='hot', \n",
    "                                   text_a='GT T1 map', \n",
    "                                   text_b='Estimated T1 map'\n",
    "                                  )\n",
    "\n",
    "imageUtilities.compare_predictions(gt_t1Map*mask, \n",
    "                                   corrupted_t1Map*mask,\n",
    "                                   colormap_c='hot', \n",
    "                                   text_a='GT T1 map', \n",
    "                                   text_b='Estimated T1 map'\n",
    "                                  )\n",
    "\n",
    "imageUtilities.compare_predictions(original_t1Map*mask,\n",
    "                                   corrupted_t1Map*mask, \n",
    "                                   colormap_c='hot', \n",
    "                                   text_a='GT T1 map',\n",
    "                                   text_b='Estimated T1 map'\n",
    "                                  )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Checking the residuals\n",
    "An effective way to verify if the mapping was performed correctly is to check the residuals of the fitted data. We can use the signal model to generate simulated weighted images, which can be compared to the original observed data. This shows us if the estimated parameters are, indeed, possible solutions to explain the observed data. The optimal residual would yield only noise (i.e. no brain structure should be visible)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Calculate and plot the residuals for the non-corrupted estimates\n",
    "residualOriginal = estimatedWeightedSeriesOriginal - originalSeries[\"weighted_series\"]\n",
    "imageBrowse(residualOriginal, \n",
    "            im_ax = 0, \n",
    "            rang_min = -0.1, \n",
    "            rang_max = 0.1,\n",
    "            colormap_c='jet', \n",
    "            fig_title='Fitting residual',\n",
    "            ylabel='[a.u.]'\n",
    "           )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Calculate and plot the residuals for the motion corrupted estimates\n",
    "residualCorrupted = estimatedWeightedSeriesCorrupted - originalSeries[\"weighted_series\"]\n",
    "imageBrowse(residualCorrupted, \n",
    "            im_ax = 0, \n",
    "            rang_min = -0.3, \n",
    "            rang_max = 0.3, \n",
    "            colormap_c='jet',\n",
    "            fig_title='Fitting residual',\n",
    "            ylabel='[a.u.]'\n",
    "           )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Motion compensation - Image registration\n",
    "As we observed previously, patient motion during the scans can have a negative effect in the estimated maps. Image registration, as performed by Elastix and SimpleElastix, can help with this by aligning the images to each other, reducing the spatial displacement of voxels.\n",
    "\n",
    "Now, we will use SimpleElastix to perform registration of the corrupted images from the example above. We use SimpleElastix because its integration with Python (and many other languages) is easy. Additionally, the object-oriented structure in which the tool was written gives an intuitive way to use it.\n",
    "\n",
    "#### Parameter File\n",
    "As explained previously, Elastix uses a set of instructions to execute. There are several ways to tell Elastix how to perform the registration. Here, for simplicity, we will use a parameter file available in the extensive community-built database. For your convenience, the file is already located in the working directory.\n",
    "\n",
    "However, if your application demands a different type of data structure, you can go to http://elastix.bigr.nl/wiki/index.php/Parameter_file_database and check what the community made available to us!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the elastix wrapper\n",
    "def elastixRegister(fixedImage, movingImage, parameterMapPath):\n",
    "    print(\"Applying registration\")\n",
    "    \n",
    "    elastixImageFilter = sitk.ElastixImageFilter()\n",
    "    elastixImageFilter.SetFixedImage(fixedImage)\n",
    "    elastixImageFilter.SetMovingImage(movingImage)\n",
    "    \n",
    "    parameterMap = sitk.ReadParameterFile(parameterFilePath)\n",
    "    elastixImageFilter.SetParameterMap(parameterMap)\n",
    "    elastixImageFilter.Execute()\n",
    "    \n",
    "    resultImage = elastixImageFilter.GetResultImage()\n",
    "    transformParameterMap = elastixImageFilter.GetTransformParameterMap()\n",
    "    \n",
    "    procData = sitk.GetArrayFromImage(resultImage)\n",
    "    procData = np.transpose(procData, (0,2,1))\n",
    "    \n",
    "    print(\"Registration finished\")\n",
    "    \n",
    "    return procData, transformParameterMap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Perform registration\n",
    "filePath = \"./support_data/synthetic_t1w_deformed.nii\"\n",
    "parameterFilePath = \"./ParametersFile/par_groupwise/par_real_data/par_groupwise_DTI-BRAIN.txt\"\n",
    "\n",
    "sitk_fixed, sitk_moving = dataManage.prepareDataForElastix(filePath)\n",
    "processedData, transformParameterMap = elastixRegister(sitk_fixed, sitk_moving, parameterFilePath)\n",
    "\n",
    "registeredSeries = {}\n",
    "registeredSeries['weighted_series'] = processedData\n",
    "registeredSeries['echo_times'] = originalSeries['echo_times']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualising registered data\n",
    "You can now visualise the registered data and compare it to the corrupted and original images. Notice how the signal evolution curve changes for the same voxel in each case. Additionally, although much better than the corrupted images, you still might see a slight misalignment between the registered images. Registration is a hard problem and, without fine tunning of the parameters, results are often suboptimal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "imageBrowse(originalSeries[\"weighted_series\"], \n",
    "            im_ax = 0, \n",
    "            rang_min = 0.0, \n",
    "            rang_max = 1.5,\n",
    "            colormap_c='gray', \n",
    "            fig_title='Original T1w',\n",
    "            ylabel='[a.u.]'\n",
    "           )\n",
    "\n",
    "imageBrowse(corruptedSeries[\"weighted_series\"], \n",
    "            im_ax = 0, \n",
    "            rang_min = 0.0, \n",
    "            rang_max = 1.5, \n",
    "            colormap_c='gray', \n",
    "            fig_title='Motion curr. T1w',\n",
    "            ylabel='[a.u.]'\n",
    "           )\n",
    "\n",
    "imageBrowse(registeredSeries[\"weighted_series\"], \n",
    "            im_ax = 0, \n",
    "            rang_min = 0.0,\n",
    "            rang_max = 1.5, \n",
    "            colormap_c='gray', \n",
    "            fig_title='Registered T1w series',\n",
    "            ylabel='[a.u.]'\n",
    "           )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### T1 mapping of registered images\n",
    "After registration, we can now perform the T<sub>1</sub> mapping of the corrected data. We expect improved results compared to the motion corrupted example. However, due to the residual misalignments observed, the T<sub>1</sub> map from the original data has higher accuracy. We can check that by evaluating the residuals of the prediction, as done before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# T1 mapping with registered images\n",
    "qMaps, estimatedWeightedSeriesRegistered = t1Mapping(registeredSeries)\n",
    "registered_t1Map = qMaps[-1].detach().numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Calculate and plot the residuals for the registered estimates\n",
    "residualRegistered = estimatedWeightedSeriesRegistered - originalSeries[\"weighted_series\"]\n",
    "imageBrowse(residualRegistered, \n",
    "            im_ax = 0, \n",
    "            rang_min = 0.0, \n",
    "            rang_max = 0.1, \n",
    "            colormap_c='hot', \n",
    "            fig_title='Motion-free residual',\n",
    "            ylabel='[a.u.]'\n",
    "           )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ---------------------------------------------------------------------------------------\n",
    "# Now it's your turn!\n",
    "\n",
    "In the folder './support_data/' you will find an example of in-vivo data, acquired using an Inversion Recovery sequence. The data structure is similar to what was explained above. Using similar steps as before, you should:\n",
    "\n",
    "1. Load and visualise your data\n",
    "2. Do the T1 mapping with the original data and check the quality of fit (Check the residual)\n",
    "3. Perform the image registration (you can use the same parameter file)\n",
    "4. Again, perform the T1 mapping with the registered data\n",
    "5. Verify results and compare the estimates of step 4 with the original estimates (step 2)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Loading In Vivo data\n",
    "Implemented for you :)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## This part is given, since the loading method is different.\n",
    "## Load In-vivo data\n",
    "filePath = \"./support_data/inVivo.mat\"\n",
    "inVivoSeries = dataManage.loadRealData(filePath, 10)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualise the data\n",
    "Implement me..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imageBrowse(inVivoSeries[\"weighted_series\"], \n",
    "            im_ax = 0,\n",
    "            rang_min = 0.0,\n",
    "            rang_max = 1.5,\n",
    "            ylabel='[a.u.]'\n",
    "           )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Perform T1 mapping on original data\n",
    "Implement me..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Check residuals\n",
    "Implement me..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Perform registration\n",
    "Implement me..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualise the registered data\n",
    "Implement me..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Performing T1 mapping on registered data\n",
    "Implement me..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Check residuals\n",
    "Implement me..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Compare original and registered maps\n",
    "Implement me..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}