Feat: Implement a visualisation feature to NiTransforms

.circleci/config.yml

@@ -142,7 +142,8 @@ jobs:
               -v /tmp/data/nitransforms-tests:/data -e TEST_DATA_HOME=/data \
               -e COVERAGE_FILE=/tmp/summaries/.pytest.coverage \
               -v /tmp/fslicense/license.txt:/opt/freesurfer/license.txt:ro \
-              -v /tmp/tests:/tmp nitransforms:latest \
+              -v /tmp/tests:/tmp -e TEST_OUTPUT_DIR=/tmp/artifacts \
+              nitransforms:latest \
               pytest --junit-xml=/tmp/summaries/pytest.xml \
                      --cov nitransforms --cov-report xml:/tmp/summaries/unittests.xml \
                      nitransforms/

docs/examples.rst

@@ -9,3 +9,4 @@ A collection of Jupyter Notebooks to serve as interactive tutorials.
 
     notebooks/isbi2020
     notebooks/Reading and Writing transforms.ipynb
+    notebooks/Visualizing transforms.ipynb

docs/notebooks/Visualizing transforms.ipynb

@@ -0,0 +1,379 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Vis - visualizing transforms\n",
+    "This notebook showcases the `nitransforms.vis` module, which implements the functionality to illustrate and view transforms."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Preamble\n",
+    "Prepare a Python environment and use a temporal directory for the outputs. After that, fetch the actual file from NiBabel documentation."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import os\n",
+    "from pathlib import Path\n",
+    "from tempfile import TemporaryDirectory\n",
+    "\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "import nibabel as nb\n",
+    "from nitransforms.vis import PlotDenseField\n",
+    "\n",
+    "cwd = TemporaryDirectory()\n",
+    "os.chdir(cwd.name)\n",
+    "print(f\"This notebook is being executed under <{os.getcwd()}>.\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "test_dir = \"Desktop/nitransforms-tests\"\n",
+    "\n",
+    "anat_file = Path(os.getenv(\"TEST_DATA_HOME\", str(Path.home() / test_dir))) / \"someones_anatomy.nii.gz\"\n",
+    "transform_file = Path(os.getenv(\"TEST_DATA_HOME\", str(Path.home() / test_dir))) / \"someones_displacement_field.nii.gz\"\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Load a displacement field\n",
+    "Info about the transform file here..."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Load the example\n",
+    "nii = nb.load(anat_file)\n",
+    "hdr = nii.header.copy()\n",
+    "aff = nii.affine.copy()\n",
+    "data = np.asanyarray(nii.dataobj)\n",
+    "nii.orthoview()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Displacement fields\n",
+    "About displacement fields here..."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Using `PlotDenseField()` from the `nitransforms.vis` module, we can prepare the ground for the illustration of a dense field transform object. `PlotDenseField()` takes three aruments: the transform file, an indication regarding the nature of the transform (default: deltas field) and a reference (default: None). We must also define slices to select which axial, coronal and sagittal planes we wish to visualise.\n",
+    "\n",
+    "Looking at the corresponding field for the example nifti file showcased above:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "xslice, yslice, zslice = 30, 30, 30\n",
+    "\n",
+    "pdf = PlotDenseField(\n",
+    "    transform=transform_file,\n",
+    "    is_deltas=True,\n",
+    ")\n",
+    "print(pdf, pdf._xfm._field.shape)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Plot distortion\n",
+    "To sbegin, we can use the `DenseFieldTransform().plot_distortion` module to visualise the deformed grid and superimpose it onto the density map of the deformation. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axes = plt.subplots(1, 3, figsize=(12, 4), tight_layout=True)\n",
+    "pdf.plot_distortion(\n",
+    "    axes=axes,\n",
+    "    xslice=xslice,\n",
+    "    yslice=yslice,\n",
+    "    zslice=zslice,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In this example, we notice that the coronal slice contains no deformations, hinting that the transform is only dependent on the y-dimension, ie that the deltas field only contains non-zero $V_y$ vector-components. The distortions of both Axial and Sagittal planes confirm this, with the grid being unaffected in $x$ and $z$ (ie $V_{x, i} = V_{y, i} = 0$ for all $i$). \n",
+    "\n",
+    "## Quiver deformation scalar map\n",
+    "To verify this, we can plot the field using `DenseFieldTransform().plot_quiverdsm` and view the transformation as a quiver plot. Here, the colors represent the dominant vector components of the transform field accroding to the standard convention for deformation scalar maps: red ($V_x$-dominant), green ($V_y$-dominant) and blue ($V_z$-dominant)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axes = plt.subplots(1, 3, figsize=(12, 4), tight_layout=True)\n",
+    "pdf.plot_quiverdsm(\n",
+    "    axes=axes,\n",
+    "    xslice=xslice,\n",
+    "    yslice=yslice,\n",
+    "    zslice=zslice,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "As expected, the field is restrained to a single spatial dimension ($V_y$) and is most intense in the areas where gridlines are subject to the most distortions. The magnitude of the field at any point is represented by the length of the correspodning arrow, while the colour represence the dominance of the dimensional component. In this example, given that the transformation is confined to one dimension, these measures are equivalent and, indeed, the color intensity is also a proxy for arrow length. As we will see further down, this is not always true when dealing with non-zero, three-dimensional transformation fields."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Transformation fields in three dimensions\n",
+    "Let's now look at a more complicated dense field transform, this time containing deformations in all three spatial dimensions. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "transform_3d = Path(os.getenv(\"TEST_DATA_HOME\", str(Path.home() / \"workspace/nitransforms/nitransforms/tests/data\"))) / \"ds-005_sub-01_from-OASIS_to-T1_warp_fsl.nii.gz\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "xslice, yslice, zslice = 50, 50, 50\n",
+    "\n",
+    "pdf_3d = PlotDenseField(\n",
+    "    transform=transform_3d,\n",
+    "    is_deltas=True,\n",
+    ")\n",
+    "print(pdf_3d, pdf_3d._xfm._field.shape)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can start by taking a look at the deformation grid."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axes = plt.subplots(1, 3, figsize=(12, 4), tight_layout=True)\n",
+    "pdf_3d.plot_distortion(\n",
+    "    axes=axes,\n",
+    "    xslice=xslice,\n",
+    "    yslice=yslice,\n",
+    "    zslice=zslice,\n",
+    "    lw=0.25,\n",
+    "    show_brain=False,\n",
+    "    show_grid=True,\n",
+    ")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The three-dimensional deformations can easily be identified through the distortion of the coordinate grid alone. \n",
+    "\n",
+    "As in 1D, we again want to see how the field looks as a quiver plot."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Quiver and deformation scalar map\n",
+    "Once again, it is straightfowrad to highlight the dominance of each dimension using `PlotDensefield().plot_quiverdsm`... "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axes = plt.subplots(1, 3, figsize=(12, 4), tight_layout=True)\n",
+    "pdf_3d.plot_quiverdsm(\n",
+    "    axes = axes,\n",
+    "    xslice=xslice,\n",
+    "    yslice=yslice,\n",
+    "    zslice=zslice\n",
+    ")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "... which now allows us to identify which regions of the transform field have dominant $x$- (red), $y$- (green) or $z$- (blue) vector components.\n",
+    "\n",
+    "The color of each arrow therefore highlights which of the three directional components ($V_x$, $V_y$ or $V_z$) is dominant. The intensity of the color maps the strength of the dominant dimension (eg pale red -> dark red represents low $V_x$ -> high $V_x$, and similarly for $y$- and $z$- components with green and blue, respectively). The magnitude of the vector is still represented by the arrow length. \n",
+    "\n",
+    "For example, a vector $\\textbf{V} = (0.1, 0.2, 0.5)$ with magnitude $V \\approx 0.5$ has a slightly dominant $V_z$ and a relatively low magnitude, so will appear pale blue with a relatively short arrow. On the other hand, although a different vector $\\textbf{V} = (7.2, 7.4, 7.8)$ with $V \\approx 13$ still has a very small $z$-dominace, the $V_z$ component is significantly larger causing the arrow to take a dark blue color. Independently of the color, this vector will be long due to it's large magnitude."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Linear, planar and spherical coefficients\n",
+    "We have shown that the transformation field can be mapped according to the dominance of its individual vector components. However, in the case where the difference between the most dominant vector component and the second most dominant component is very small, the transformation is in fact dictated by two dimensions. \n",
+    "\n",
+    "We therefore compute the Linear, Planar and Spherical coefficients of each point in the vector field to visualise whether the transformation is dictated by one, two or three dimensions, respectively. \n",
+    "\n",
+    "For example:\n",
+    "- A vector $V = (3, 0.2, 0.4)$ dominant in one dimension ($V_x$) will have a linear coefficient $c_L > c_p, c_s$\n",
+    "- A vector $V = (3, 2.8, 0.4)$ dominant in two dimensions ($V_x$ and $V_y$) will have a planar coefficient $c_p > c_l, c_s$\n",
+    "- A vector $V = (3, 2.8, 3.1)$ dominant in three dimensions will have a spherical coefficient $c_s > c_l, c_p$\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axes = plt.subplots(3, 3, figsize=(10, 9), layout='constrained')\n",
+    "pdf_3d.plot_coeffs(\n",
+    "    fig=fig,\n",
+    "    axes=axes,\n",
+    "    xslice=xslice,\n",
+    "    yslice=yslice,\n",
+    "    zslice=zslice,\n",
+    ")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Jacobians\n",
+    "Now that our field has components in more than one dimension, we can also use `PlotDenseField().plot_jacobian` to map out the jacobians of the vector field and highlight regions in which the transformation represents an expansion (red, $J>0$) or a contraction (blue, $J<0$) of the brain."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axes = plt.subplots(1, 3, figsize=(12, 4.5), tight_layout=True)\n",
+    "pdf_3d.plot_jacobian(\n",
+    "    axes=axes,\n",
+    "    xslice=xslice,\n",
+    "    yslice=yslice,\n",
+    "    zslice=zslice,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Transform overview\n",
+    "Now that we have unpacked the dense field transform, we can now bring this all together with `PlotDenseField().show_transform`. This creates a 3x3 grid of plots to provide the user with an overview of the transformation field contained in the nifti file."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pdf_3d.show_transform(\n",
+    "    xslice=xslice,\n",
+    "    yslice=yslice,\n",
+    "    zslice=zslice,\n",
+    ")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "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.12.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}