{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 10. Transformix: Spatial Jacobian calculations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the transformix algorithm the spatial jacobian and the determinant of the spatial jacobian of the transformation can be calculated.\n",
    "Especially the determinant of the spatial Jacobian, which identifies the amount of local\n",
    "compression or expansion and can be useful, for example in lung ventilation studies."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Elastix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import itk\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Images\n",
    "fixed_image = itk.imread('data/CT_2D_head_fixed.mha', itk.F)\n",
    "moving_image = itk.imread('data/CT_2D_head_moving.mha', itk.F)\n",
    "\n",
    "# Import Default Parameter Map\n",
    "parameter_object = itk.ParameterObject.New()\n",
    "parameter_map_rigid = parameter_object.GetDefaultParameterMap('rigid')\n",
    "parameter_object.AddParameterMap(parameter_map_rigid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Registration with the registration function.\n",
    "The output directory has to be specified, \n",
    "otherwise elastix will not save the transformparameter file as .txt file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Call registration function and specify output directory\n",
    "result_image, result_transform_parameters = itk.elastix_registration_method(\n",
    "    fixed_image, moving_image,\n",
    "    parameter_object=parameter_object,\n",
    "    output_directory='exampleoutput/')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transformix Jacobian Calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Image to transform, transformix is transforming from moving -> fixed;\n",
    "# for this example the exact same moving image is used, this however is normally not \n",
    "# very usefull since the elastix algorithm already transformed this image.\n",
    "moving_image_transformix = itk.imread('data/CT_2D_head_moving.mha', itk.F)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The calculation of the Jacobian matrix and it's determinant can be done in one line..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate Jacobian matrix and it's determinant in a tuple\n",
    "jacobians = itk.transformix_jacobian(moving_image_transformix, result_transform_parameters)\n",
    "\n",
    "# Casting tuple to two numpy matrices for further calculations.\n",
    "spatial_jacobian = np.asarray(jacobians[0]).astype(np.float32)\n",
    "det_spatial_jacobian = np.asarray(jacobians[1]).astype(np.float32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ".. or by initiating an transformix image filter object, calculating the jacobian and reading the jacobian from Disk IO."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load Transformix Object\n",
    "transformix_object = itk.TransformixFilter.New(moving_image_transformix)\n",
    "transformix_object.SetTransformParameterObject(result_transform_parameters)\n",
    "\n",
    "# Set advanced options\n",
    "transformix_object.SetComputeSpatialJacobian(True)\n",
    "transformix_object.SetComputeDeterminantOfSpatialJacobian(True)\n",
    "\n",
    "# Set output directory for spatial jacobian and its determinant,\n",
    "# default directory is current directory.\n",
    "transformix_object.SetOutputDirectory('exampleoutput/')\n",
    "\n",
    "# Update object (required)\n",
    "transformix_object.UpdateLargestPossibleRegion()\n",
    "\n",
    "# Results of Transformation\n",
    "result_image_transformix = transformix_object.GetOutput()\n",
    "\n",
    "# Load Jacobian from Disk IO and cast to numpy matrices\n",
    "spatial_jacobian = itk.imread('exampleoutput/fullSpatialJacobian.nii', itk.F)\n",
    "spatial_jacobian = np.asarray(spatial_jacobian).astype(np.float32)\n",
    "\n",
    "det_spatial_jacobian = itk.imread('exampleoutput/spatialJacobian.nii', itk.F)\n",
    "det_spatial_jacobian = np.asarray(det_spatial_jacobian).astype(np.float32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inspect the deformation field by looking at the determinant of the Jacobian of Tµ(x). Values smaller\n",
    "than 1 indicate local compression, values larger than 1 indicate local expansion, and 1 means volume\n",
    "preservation. The measure is quantitative: a value of 1.1 means a 10% increase in volume. If this\n",
    "value deviates substantially from 1, you may be worried (but maybe not if this is what you expect for\n",
    "your application). In case it is negative you have “foldings” in your transformation, and you definitely\n",
    "should be worried. For more information see [elastix manual](https://elastix.lumc.nl/download/elastix-5.0.1-manual.pdf)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of foldings in transformation: 0\n"
     ]
    }
   ],
   "source": [
    "print(\"Number of foldings in transformation:\",np.sum(det_spatial_jacobian < 0))"
   ]
  }
 ],
 "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.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}