{
 "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": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First two import are currently necessary to run ITKElastix on MacOs\n",
    "from itk import itkElastixRegistrationMethodPython\n",
    "from itk import itkTransformixFilterPython\n",
    "import itk\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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": 3,
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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": [
    "Transformation can either be done in one line with the transformix function..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "transformed_image = itk.transformix_filter(\n",
    "    moving_image=moving_image_transformix,\n",
    "    transform_parameter_object=result_transform_parameters,\n",
    "    compute_spatial_jacobian=True,\n",
    "    compute_determinant_of_spatial_jacobian=True,\n",
    "    output_directory='exampleoutput/')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ".. or by initiating an transformix image filter object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load Transformix Object\n",
    "transformix_object = itk.TransformixFilter.New()\n",
    "transformix_object.SetMovingImage(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()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Jacobian Matrix and it's determinant"
   ]
  },
  {
   "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](http://elastix.isi.uu.nl/download/elastix-5.0.0-manual.pdf)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of foldings in transformation: 0\n"
     ]
    }
   ],
   "source": [
    "det_spatial_jacobian = itk.imread('exampleoutput/spatialJacobian.nii', itk.F)\n",
    "det_spatial_jacobian = np.asarray(det_spatial_jacobian).astype(np.float32)\n",
    "\n",
    "print(\"Number of foldings in transformation:\",np.sum(det_spatial_jacobian < 0))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ElastixEnv",
   "language": "python",
   "name": "elastixenv"
  },
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}