{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f62c3ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "# -*- coding: utf-8 -*-\n",
    "#  Copyright 2022 United Kingdom Research and Innovation\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:    Laura Murgatroyd (UKRI-STFC)\n",
    "#                   Gemma Fardell (UKRI-STFC)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "694dafc7",
   "metadata": {},
   "source": [
    "# Exercise 02 - Preprocessing with the Core Imaging Library (CIL) \n",
    "\n",
    "### 3D laboratory micro-CT, cone-beam data of sunflower seeds in an acrylic box"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "127515ba",
   "metadata": {},
   "source": [
    "This exercise walks through the steps needed to load in a 3D cone-beam dataset of sunflower seeds in an acrylic box, acquired by laboratory micro-CT, pre-process, and reconstruct it using FDK.\n",
    "\n",
    "Notice, this uses the same sample as in [01_intro_seeds_conebeam.ipynb](./01_intro_seeds_conebeam.ipynb). However, in that notebook, the dataset file we used had already been altered to contain the centre of rotation offset. Here we use a file which has not had that applied, so we need to establish the centre of rotation offset ourselves, using CIL.\n",
    "\n",
    "Learning objectives are:\n",
    "- Load and investigate a Nikon data set.\n",
    "- Apply CIL's `TransmissionAbsorptionConverter`.\n",
    "- Apply CIL's Centre of Rotation corrector.\n",
    "- Compute FDK reconstruction using CIL.\n",
    "- Re-bin a dataset using CIL's `Binner` processor.\n",
    "\n",
    "This example requires the dataset `korn.zip` from https://zenodo.org/record/6874123#.Y0ghJUzMKUm :\n",
    "\n",
    "- https://zenodo.org/record/6874123/files/korn.zip\n",
    "\n",
    "If running locally please download the data and update the filepath in the `filename` variable below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89d2596e",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = \"/mnt/materials/SIRF/Fully3D/CIL/Korn i kasse/47209 testscan korn01.xtekct\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d48350f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from cil.io import NikonDataReader\n",
    "from cil.processors import TransmissionAbsorptionConverter, Slicer, CentreOfRotationCorrector, Binner\n",
    "from cil.recon import FDK\n",
    "from cil.utilities.display import show2D, show_geometry\n",
    "from cil.utilities.jupyter import islicer"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23520522-9dd7-4333-a370-e550db6e2008",
   "metadata": {},
   "source": [
    "Here we turn on logging for CIL's processors. This means we will get more detailed information when running the processors. This is especially useful when calculating the centre of rotation offset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f10df095-0f01-420c-a0d3-b5b6f41252e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging\n",
    "logging.basicConfig(level=logging.WARNING)\n",
    "cil_log_level = logging.getLogger('cil.processors')\n",
    "cil_log_level.setLevel(logging.INFO)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d24d898-d1b2-4496-a464-6a2e9dbf8f89",
   "metadata": {},
   "source": [
    "## Exercise A: Loading Nikon Data and looking at the Geometry"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bac0c33",
   "metadata": {},
   "source": [
    "1. Load the 3D cone-beam projection data of the seeds, using the `NikonDataReader`\n",
    "2. `print` the data to get some basic information.\n",
    "3. As well as the data itself, AcquisitionData contains geometric metadata in an AcquisitionGeometry object in the geometry field. `print` the geometry data.\n",
    "4. Use the `show_geometry` method to display the scan set up visually.\n",
    "\n",
    "*Note: This is a full 3D dataset so reading it from disk may take some time*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b30864da",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_in = ... \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6bce18d2",
   "metadata": {},
   "source": [
    "The data is loaded in as a CIL `AcquisitionData` object. How many projections does this dataset contain and how many pixels do they have? Make sure to check the axis labels."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9caa103-69a0-4696-955d-dcff26b99118",
   "metadata": {},
   "source": [
    "**Uncomment the following line and run the cell to see the solution, to run the lines you'll need to run the cell a second time**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f2d6ac1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load './snippets/02_exA.py'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15e657e6-1e14-49be-a762-02e632468be9",
   "metadata": {},
   "source": [
    "## Exercise B: Displaying the Projections with islicer"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8ab7f23",
   "metadata": {},
   "source": [
    "Use `islicer` to display the projections."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48d432ca-a8b3-4e0a-94b8-bfad3d8705dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27469287-1c26-46b9-8f3e-8e85e990c551",
   "metadata": {},
   "source": [
    "Uncomment the following line to see the solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acf157f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load './snippets/02_exB.py'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79597e26-2cdc-45fc-935c-426fe4c8a3c2",
   "metadata": {},
   "source": [
    "## Exercise C: Transmission to Absorption Conversion"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "937f501c-564c-4c2c-82d9-700388715438",
   "metadata": {},
   "source": [
    "You should have seen that the data is transmission data. We know this because the background value is 1.0. We need to apply the Beer–Lambert law to convert to the absorption data.\n",
    "\n",
    "1. Use CIL's Transmission to Absorption processor to convert the data to absorption.\n",
    "2. Use show2D to look at the central `vertical` slice of the absorption data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad8f9f5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_absorption = ..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ae5d12f-9bfd-48ca-a2c3-dfd7fb03a82d",
   "metadata": {},
   "source": [
    "Uncomment the following line to see the solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fc02a43-f80c-4f20-ad95-b82aecb471a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load './snippets/02_exC.py'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d05b90f8-6859-4914-bfb3-0e2f13688766",
   "metadata": {},
   "source": [
    "## Exercise D: Reconstructing with FDK"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4899ed4d",
   "metadata": {
    "tags": []
   },
   "source": [
    "We will use the `FDK` algorithm from CIL's recon module. FDK is filtered back-projection with special weights for cone-beam data. By default, the `recon` module uses TIGRE as a back-end. We will use \n",
    "\n",
    "1. Use reorder to ensure the data is in the correct format for tigre\n",
    "2. Create and run the `FDK` algorithm, using our `image_geometry` created below.\n",
    "3. Then show the reconstructed volume using `islicer`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a07da76f-e255-47ce-b4ee-2dfaecc6b3eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "image_geometry = data_absorption.geometry.get_ImageGeometry()\n",
    "image_geometry.voxel_num_x = 700\n",
    "image_geometry.voxel_num_y = 700\n",
    "image_geometry.voxel_num_z = 700"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "592171ef-c310-44c8-ae91-796e887be7db",
   "metadata": {},
   "outputs": [],
   "source": [
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6de2add5-c745-4bf3-a270-5671395d7548",
   "metadata": {},
   "source": [
    "Uncomment the following line to see the solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0ae0659",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# %load ./snippets/02_exD.py"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1603228-8d17-4cc0-af42-146e6d35e1ef",
   "metadata": {},
   "source": [
    "## Exercise E: Performing Centre of Rotation Correction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d03f4f89-2165-4cab-a157-9dcdb9ce66d8",
   "metadata": {},
   "source": [
    "You should notice that the above reconstruction does not look right. This edge-doubling is a classic artifact from a centre of rotation offset. \n",
    "\n",
    "In a perfectly aligned CT system the projection of the axis of rotation onto the detector is aligned with the horizontal centre of the detector. In practise it is not usually perfectly aligned. A slight offset of the centre of rotation with respect to the theoretical position used in the reconstruction will contribute to the loss of resolution; in more severe cases it will cause the severe artifacts in the reconstructed volume we see above.\n",
    "\n",
    "We can estimate the true centre of rotation offset from the acquisition data using CIL's `CentreOfRotationCorrector`. Here we will use CIL's the `image_sharpness` algorithm as it works well on cone-beam data.\n",
    "\n",
    "1. Use CIL to calculate and apply the centre of rotation correction and find out how many pixels offset it calculated\n",
    "2. Use `show_geometry` and `print` to compare the geometry before and after the correction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adb96857-8cd1-4a12-a56f-347464bae522",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_centred = ..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3777d725-00ee-48ff-9452-4cfe4e7b41b1",
   "metadata": {},
   "source": [
    "Uncomment the following line to see the solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41f3dc8e-6223-4428-868b-ccfe9f2c1d1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load ./snippets/02_exE.py"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3302677c-8613-451a-a7fd-06d08061b6f3",
   "metadata": {},
   "source": [
    "## Exercise F: Reconstruct the Centre of Rotation Corrected Data\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab8325e8-5ad4-4d50-848d-8248d3905b05",
   "metadata": {},
   "source": [
    "Now that we have applied the centre of rotation correction, perform the FDK reconstruction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "326f311e-40e9-4c16-b785-5ba3bc9a2db0",
   "metadata": {},
   "outputs": [],
   "source": [
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e656e556-2102-459a-894f-27a9505a67ff",
   "metadata": {},
   "source": [
    "Uncomment the following line to see the solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f922634a-9ce2-4026-aef6-972cc2f35106",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0dda174d-75fb-4b05-93b1-1fa38ac12b62",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load ./snippets/02_exF.py"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3492e36-c389-435e-a948-ae7c92d03348",
   "metadata": {},
   "source": [
    "## Exercise G: Binning Data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c2087a7-b3a9-443f-83a9-13b7b8287bef",
   "metadata": {},
   "source": [
    "Start again from using the dataset we read from the file (before we applied any processing).\n",
    "Re-bin the data 4x along both the horizontal and vertical axes.\n",
    "Refer to the [CIL documentation](https://tomographicimaging.github.io/CIL/nightly/processors.html?highlight=binner#cil.processors.Binner) for how to set up the `Binner` processor.\n",
    "\n",
    "Then process (transmission to absorption convert and centre of rotation correct) the binned data and reconstruct with FDK:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "418065b1-e9a3-46c9-8f70-75d942739477",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "fadb78ff",
   "metadata": {},
   "source": [
    "Uncomment the following line to see the solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c95bd28-4bb6-41ac-b366-b57bcf723ca2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load ./snippets/02_exG.py"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.6"
  },
  "vscode": {
   "interpreter": {
    "hash": "85dd3531d2b361567244a097ef0921818e021a508f13857e8138f8c4801481e0"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}