{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# -*- coding: utf-8 -*-\n",
    "#  Copyright 2021 - 2022 United Kingdom Research and Innovation\n",
    "#  Copyright 2021 - 2022 The University of Manchester\n",
    "#\n",
    "#  Licensed under the Apache License, Version 2.0 (the \"License\");\n",
    "#  you may not use this file except in compliance with the License.\n",
    "#  You may obtain a copy of the License at\n",
    "#\n",
    "#      http://www.apache.org/licenses/LICENSE-2.0\n",
    "#\n",
    "#  Unless required by applicable law or agreed to in writing, software\n",
    "#  distributed under the License is distributed on an \"AS IS\" BASIS,\n",
    "#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
    "#  See the License for the specific language governing permissions and\n",
    "#  limitations under the License.\n",
    "#\n",
    "#   Authored by:    Gemma Fardell (UKRI-STFC)\n",
    "#                   Edoardo Pasca (UKRI-STFC)\n",
    "#                   Laura Murgatroyd (UKRI-STFC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A detailed look at CIL geometry\n",
    "CIL holds your CT data in specialised data-containers, `AcquisitionData` and `ImageData`.\n",
    "\n",
    "Each of these has an associated `geometry` which contains the meta-data describing your set-up.\n",
    "\n",
    " - `AcquisitionGeometry` describes the acquisition data and parameters\n",
    "\n",
    " - `ImageGeometry` describes the image data (i.e., the reconstruction volume)\n",
    "\n",
    "The data-readers provided by CIL (Nikon, Zeiss and diamond nexus readers) will read in your data and return you a fully configured acquisition data with the acquisition geometry already configured, however if you read in a stack of tiffs or want to tweak the parameters this is simple to create by hand."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The structure of an AcquisitionGeometry\n",
    "\n",
    "An instance of an `AcquisitionGeometry`, `ag`, holds the configuration of the system, in `config` which is subdivided in to:\n",
    " - `ag.config.system` -  The position and orientations of the `source`/`ray`, `rotation_axis` and `detector`\n",
    " - `ag.config.panel` - The number of pixels, the size of pixels, and the position of pixel 0\n",
    " - `ag.config.angles` - The number of angles, the unit of the angles (default is degrees)\n",
    " - `ag.config.channels` - The number of channels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a simple AcquisitionGeometry\n",
    "\n",
    "You can use the `AcquisitionGeometry` methods to describe circular trajectory parallel-beam or cone-beam 2D or 3D data.\n",
    "\n",
    " - `ag = AcquisitionGeometry.create_Parallel2D()`\n",
    " - `ag = AcquisitionGeometry.create_Parallel3D()`\n",
    " - `ag = AcquisitionGeometry.create_Cone2D(source_position, detector_position)`\n",
    " - `ag = AcquisitionGeometry.create_Cone3D(source_position, detector_position)`\n",
    "\n",
    "This notebook will step though each in turn and show you how to describe both simple and complex geometries with offsets and rotations.\n",
    "\n",
    "No matter which type of geometry you create you will also need to describe the panel and projection angles.\n",
    " - `ag.set_panel(num_pixels, pixel_size)`\n",
    " - `ag.set_angles(angles, angle_unit)`\n",
    "\n",
    "For multi-channel data you need to add the number of channels.\n",
    " - `ag.set_channels(num_channels)`\n",
    "\n",
    "And you will also need to describe the order your data is stored in using the relavent labels from the CIL default labels: `channel`, `angle`, `vertical` and `horizontal`\n",
    " - `ag.set_labels(['angle','vertical','horizontal'])`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A Note on CIL AcquisitionGeometry:\n",
    " - The geometry is described by a right-handed cooridinate system\n",
    " - Positive angles describe the object rotating anti-clockwise when viewed from above\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parallel geometry\n",
    "\n",
    "Parallel beams of X-rays are emitted onto 1D (single pixel row) or 2D detector array. This geometry is common for synchrotron sources.\n",
    "\n",
    "We describe the system, and then set the panel and angle data. Note that for 3D geometry we need to describe a 2D panel where `num_pixels=[X,Y]`\n",
    "\n",
    "```python\n",
    "parallel_2D_geometry = AcquisitionGeometry.create_Parallel2D()\\\n",
    "    \n",
    "                                          .set_panel(num_pixels=10)\\\n",
    "        \n",
    "                                          .set_angles(angles=range(0,180))\n",
    "\n",
    "\n",
    "parallel_3D_geometry = AcquisitionGeometry.create_Parallel3D()\\\n",
    "    \n",
    "                                          .set_panel(num_pixels=[10,10])\\\n",
    "        \n",
    "                                          .set_angles(angles=range(0,180))\n",
    "```\n",
    "Both 2D and 3D parallel-beam geometries are displayed below. Note that the detector position has been set, this is not necessary to describe and reconstruct the data, but it makes the displayed images clearer.\n",
    "\n",
    "`show_geometry()` can be used to display the configured geometry and will be used here extensively. You can also print the geometry to obtain a detailed description. If `show_geometry` is not passed an `ImageGeometry` it will show the default geometry associated with the `AcquisitionGeometry` \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An example creating a 2D parallel-beam geometry:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cil.framework import AcquisitionGeometry\n",
    "from cil.utilities.display import show_geometry\n",
    "\n",
    "ag = AcquisitionGeometry.create_Parallel2D(detector_position=[0,10])\\\n",
    "    .set_panel(num_pixels=10)\\\n",
    "    .set_angles(angles=range(0,180))\n",
    "\n",
    "show_geometry(ag)\n",
    "\n",
    "print(ag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An example creating a 3D parallel-beam geometry:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ag = AcquisitionGeometry.create_Parallel3D(detector_position=[0,10,0])\\\n",
    "    .set_panel(num_pixels=[10,10])\\\n",
    "    .set_angles(angles=range(0,180))\n",
    "    \n",
    "show_geometry(ag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fan-beam geometry\n",
    "\n",
    "A single point-like X-ray source emits a cone-beam onto a single row of detector pixels. The beam is typically collimated to imaging field of view. Collimation greatly reduce amount of scatter radiation reaching the detector. Fan-beam geometry is used when scattering has significant influence on image quality or single-slice reconstruction is sufficient.\n",
    "\n",
    "We describe the system, and then set the panel and angle data.\n",
    "\n",
    "For fan-beam data the source and detector positions are required. As default we place them along the Y-axis where the rotation-axis is on the origin. They are specified as `[x,y]` coordinates.\n",
    "\n",
    "```python\n",
    "cone_2D_geometry = AcquisitionGeometry.create_Cone2D(source_position=[0,-10],detector_position=[0,10])\\\n",
    "    \n",
    "                                      .set_panel(num_pixels=10)\\\n",
    "        \n",
    "                                      .set_angles(angles=range(0,180))\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ag = AcquisitionGeometry.create_Cone2D(source_position=[0,-10],detector_position=[0,10])\\\n",
    "    .set_panel(num_pixels=10)\\\n",
    "    .set_angles(angles=range(0,180))\n",
    "    \n",
    "show_geometry(ag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cone-beam geometry\n",
    "\n",
    "A single point-like X-ray source emits a cone-beam onto 2D detector array. Cone-beam geometry is mainly used in lab-based CT instruments.\n",
    "\n",
    "We describe the system, and then set the panel and angle data.\n",
    "\n",
    "For cone-beam data the source and detector positions are required. As default we place them along the Y-axis where the rotation-axis is on the origin and aligned in the Z-direction. They are specified as `[X,Y,Z]` coordinates.\n",
    "\n",
    "```python\n",
    "cone_3D_geometry = AcquisitionGeometry.create_Cone3D(source_position=[0,-10,0], detector_position=[0,10,0])\\\n",
    "    \n",
    "                                      .set_panel(num_pixels=[10,10])\\\n",
    "        \n",
    "                                      .set_angles(angles=range(0,180))\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ag = AcquisitionGeometry.create_Cone3D(source_position=[0,-10,0],detector_position=[0,10,0])\\\n",
    "    .set_panel(num_pixels=[10,10])\\\n",
    "    .set_angles(angles=range(0,180))\n",
    "    \n",
    "show_geometry(ag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create an offset AcquisitionGeometry\n",
    "\n",
    "It is unusual to have a perfectly aligned CT system. One of the most common offsets is the rotation-axis. If this offset is described by the `AcquisitionGeometry` then it will be accounted for in the reconstruction. This saves having to pad your data to account for this.\n",
    "\n",
    "To specify the offset, you could either add an x-component to the `source_position` and `detector_position` or you can offset the rotation axis from the origin using `rotation_axis_position`.\n",
    "\n",
    "As with the `source_position` and `detector_position` this is the `rotation_axis_position` is specified in 2D with a 2D vector `[X,Y]` or 3D with a 3D vector `[X,Y,Z]`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we offset the rotation axis by -0.5 in X by setting `rotation_axis_position=[-0.5,0]`. You can see the rotation axis position is no longer a point on the source-to-detector vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ag = AcquisitionGeometry.create_Cone2D(source_position=[0,-10],detector_position=[0,10],\n",
    "                                       rotation_axis_position=[-0.5,0])\\\n",
    "    .set_panel(num_pixels=10)\\\n",
    "    .set_angles(angles=range(0,180))\n",
    "    \n",
    "show_geometry(ag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a more complex AcquisitionGeometry\n",
    "\n",
    "We can also set up rotations in the system. These are configured with vectors describing the direction.\n",
    "\n",
    "For example a detector yaw can be described by using `detector_direction_x=[X,Y]`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ag = AcquisitionGeometry.create_Cone2D(source_position=[0,-10],detector_position=[0,10],\n",
    "                                       detector_direction_x=[0.9,0.1]\n",
    "                                       )\\\n",
    "     .set_panel(num_pixels=10)\\\n",
    "     .set_angles(angles=range(0,180))\n",
    "     \n",
    "show_geometry(ag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can set `rotation_axis_direction`, `detector_direction_x` and `detector_direction_y` by specifying a 3D directional vector `[X,Y,Z]`.\n",
    "\n",
    "For 3D datasets detector roll is commonly corrected with a dual-slice centre of rotation algorithm. You can specify `detector_direction_x` and `detector_direction_y` - ensuring they are ortogonal vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ag = AcquisitionGeometry.create_Cone3D(source_position=[0,-500,0],detector_position=[0,500,0],\n",
    "                                       detector_direction_x=[0.9,0.0,-0.1],detector_direction_y=[0.1,0,0.9]\n",
    "                                      )\\\n",
    "    .set_panel(num_pixels=[2048,2048], pixel_size = 0.2)\\\n",
    "    .set_angles(angles=range(0,180))\n",
    "    \n",
    "show_geometry(ag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In 3D datasets we can tilt the rotation axis to describe laminograpy geometry by changing `rotation_axis_direction`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ag = AcquisitionGeometry.create_Cone3D(source_position=[0,-500,0],detector_position=[0,500,0],rotation_axis_direction=[0,-1,1])\\\n",
    "    .set_panel(num_pixels=[2048,2048], pixel_size = 0.2)\\\n",
    "    .set_angles(angles=range(0,180))\n",
    "    \n",
    "show_geometry(ag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The structure of an ImageGeometry\n",
    "\n",
    "ImageGeometry holds the description of the reconstruction volume. It holds:\n",
    "\n",
    " - The number of voxels in X, Y, Z: `voxel_num_x`, `voxel_num_y`, `voxel_num_z`\n",
    " - The size of voxels in X, Y, Z: `voxel_size_x`, `voxel_size_y`, `voxel_size_z`\n",
    " - The offset of the volume from the rotation axis in voxels: `center_x`, `center_y`, `center_z`\n",
    " - The number of channels for multi-channel data\n",
    "\n",
    "You will also need to describe the order your data is stored in using the relevent labels from the CIL. The default labels are: `channel`, `vertical`, `horizontal_y` and `horizontal_x`\n",
    " - `ig.set_labels(['vertical','horizontal_y','horizontal_x'])`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a simple ImageGeometry\n",
    "\n",
    "To create a default ImageGeometry you can use:\n",
    " `ig = ag.get_ImageGeometry()`\n",
    "\n",
    "This creates an ImageGeometry with:\n",
    " - `voxel_num_x`, `voxel_num_y` equal to the number of horizontal pixels of the panel\n",
    " - `voxel_num_z` equal to the number of vertical pixels of the panel\n",
    " - `voxel_size_x`, `voxel_size_y` is given by the horizontal pixel size divided by magnification\n",
    " - `voxel_size_z` is given by the vertical pixel size divided by magnification\n",
    "\n",
    "\n",
    " You can pass a resolution argument:\n",
    " `ig = ag.get_ImageGeometry(resolution)` \n",
    "\n",
    " - `resolution=0.5` double the size of your voxels, and half the number of voxels in each dimension\n",
    " - `resolution=2` half the size of your voxels, and double the number of voxels in each dimension"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A Note on CIL ImageGeometry:\n",
    "At 0 degrees `horizontal_y` is aligned with the Y axis, and `horizontal_x` with the X axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ag = AcquisitionGeometry.create_Cone3D(source_position=[0,-500,0],detector_position=[0,500,0])\\\n",
    "    .set_panel(num_pixels=[2048,2048], pixel_size = 0.2)\\\n",
    "    .set_angles(angles=range(0,180))\n",
    "\n",
    "print(\"ImageGeometry - default\")\n",
    "ig = ag.get_ImageGeometry()\n",
    "print(ig)\n",
    "\n",
    "print(\"ImageGeometry - 0.5x resolution\")\n",
    "ig = ag.get_ImageGeometry(resolution=0.5)\n",
    "print(ig)\n",
    "\n",
    "print(\"ImageGeometry - 2x resolution\")\n",
    "ig = ag.get_ImageGeometry(resolution=2)\n",
    "print(ig)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a custom ImageGeometry\n",
    "You can create your own ImageGeometry with:\n",
    "`ig = ImageGeometry(...)`\n",
    "\n",
    "Giving you full control over the parameters.\n",
    "\n",
    "You can also change the members directly to reduce the reconstructed volume to exclude empty space.\n",
    "\n",
    "Using the previous example, we now can specify a smaller region of interest to reconstruct. We can offset the region of interest from the origin by specifiying the physical distance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ag = AcquisitionGeometry.create_Cone3D(source_position=[0,-500,0],detector_position=[0,500,0])\\\n",
    "    .set_panel(num_pixels=[2048,2048], pixel_size = 0.2)\\\n",
    "    .set_angles(angles=range(0,180))\n",
    "\n",
    "print(\"ImageGeometry - default\")\n",
    "ig = ag.get_ImageGeometry()\n",
    "show_geometry(ag, ig)\n",
    "\n",
    "print(\"ImageGeometry - RoI\")\n",
    "ig = ag.get_ImageGeometry()\n",
    "ig.voxel_num_z = 100\n",
    "show_geometry(ag, ig)\n",
    "\n",
    "print(\"ImageGeometry - Offset RoI\")\n",
    "ig = ag.get_ImageGeometry()\n",
    "ig.voxel_num_z = 200\n",
    "ig.center_z = -1024 * ig.voxel_size_z\n",
    "show_geometry(ag, ig)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also create an `ImageGeometry` directly.\n",
    "\n",
    "Here we create our ig independently of an `AcquisitionGeometry`, by first importing `ImageGeometry` from `cil.framework`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cil.framework import ImageGeometry\n",
    "\n",
    "ig = ImageGeometry(voxel_num_x=1000, voxel_num_y=1000, voxel_num_z=500, voxel_size_x=0.1, voxel_size_y=0.1, voxel_size_z=0.2 )"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.11",
   "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.11"
  },
  "vscode": {
   "interpreter": {
    "hash": "cf07678abc5cc77bc6e1a7d19b1e87ab0c29b83e7ee41c2bc72506d16d80ed44"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}