{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Tomography"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "This page gives an overview of the Julia package\n",
    "[`Sinograms.jl`](https://github.com/JeffFessler/Sinograms.jl).\n",
    "\n",
    "This page was generated from a single Julia file:\n",
    "[01-tomography.jl](https://github.com/JeffFessler/Sinograms.jl/blob/main/docs/lit/examples/01-tomography.jl)."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Setup"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Packages needed here."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Sinograms: fbp\n",
    "using MIRTjim: jim, prompt\n",
    "using InteractiveUtils: versioninfo"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The following line is helpful when running this example.jl file as a script;\n",
    "this way it will prompt user to hit a key after each figure is displayed."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "isinteractive() ? jim(:prompt, true) : prompt(:draw);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Tomography\n",
    "\n",
    "[Tomography](https://en.wikipedia.org/wiki/Tomography)\n",
    "is the process of imaging cross sections of an object\n",
    "without actually slicing the object.\n",
    "There are many forms of tomography;\n",
    "the description here\n",
    "focuses on\n",
    "[X-ray computed tomography (CT scans)](https://en.wikipedia.org/wiki/CT_scan).\n",
    "(See\n",
    "[SPECTrecon.jl](https://github.com/JeffFessler/SPECTrecon.jl)\n",
    "for a Julia package related to\n",
    "[SPECT](https://en.wikipedia.org/wiki/Single-photon_emission_computed_tomography)\n",
    "imaging.)\n",
    "\n",
    "In an\n",
    "[X-ray CT imaging system](http://en.wikipedia.org/wiki/File:Ct-internals.jpg),\n",
    "X-rays emitted from an X-ray source\n",
    "are transmitted through an object\n",
    "(e.g., a patient in medical CT)\n",
    "towards a detector array.\n",
    "The source and detector\n",
    "[rotate rapidly](https://www.youtube.com/watch?v=2CWpZKuy-NE)\n",
    "around the object.\n",
    "The signal intensities recorded by the detector\n",
    "are related\n",
    "to the sizes and densities\n",
    "of the object materials\n",
    "between the source and each detector element.\n",
    "\n",
    "\n",
    "## Radon transform\n",
    "\n",
    "The mathematical foundation\n",
    "for 2D X-ray CT imaging\n",
    "is the\n",
    "[Radon transform](https://en.wikipedia.org/wiki/Radon_transform),\n",
    "an\n",
    "[integral transform](https://en.wikipedia.org/wiki/Integral_transform)\n",
    "that maps a 2D function\n",
    "$f(x,y)$\n",
    "into the collection of line integrals\n",
    "through that function.\n",
    "Here we describe each line\n",
    "by its angle $\\phi$\n",
    "measured counter-clockwise from the $y$ axis,\n",
    "and by the radial distance $r$ of the line from the origin.\n",
    "The collection of integrals\n",
    "is called a sinogram.\n",
    "\n",
    "Mathematically,\n",
    "the Radon transform $p(r,\\phi)$ of $f(x,y)$ is defined by\n",
    "$$\n",
    "p(r,\\phi) = \\int_{-\\infty}^{\\infty}\n",
    "f(r \\cos \\phi - l \\sin \\phi,\n",
    "r \\sin \\phi + l \\cos \\phi) \\, \\mathrm{d} l.\n",
    "$$\n",
    "\n",
    "The Radon transform of the object that is 1 inside the unit disk\n",
    "and 0 elsewhere is given by\n",
    "$$\n",
    "p_1(r,\\phi) = 2 \\sqrt{1 - r^2}, \\ \\mathrm{ for } \\ |r| < 1.\n",
    "$$\n",
    "\n",
    "By the Radon transform's translation and scaling properties,\n",
    "the Radon transform\n",
    "of a disk of radius $r_0$\n",
    "centered at $(x_0,y_0)$\n",
    "is given by\n",
    "$$\n",
    "p(r,\\phi) = r_0 \\, p_1\n",
    "\\left( \\frac{r - (x_0 \\cos \\phi + y_0 \\sin \\phi)}{r_0}, ϕ \\right).\n",
    "$$\n",
    "\n",
    "## Sinogram example\n",
    "\n",
    "Here is a display of that Radon transform\n",
    "for a disk of radius $3$\n",
    "centered at coordinate $(5,1)$.\n",
    "Note that maximum value is approximately $6$,\n",
    "the length of the longest chord\n",
    "through a disk of radius $3$."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "nr = 128\n",
    "dr = 20 / nr\n",
    "r = ((-(nr-1)/2):((nr-1)/2)) * dr # radial sample locations\n",
    "na = 130\n",
    "ϕ = (0:(na-1))/na * π # angular samples\n",
    "proj1 = r -> abs(r) < 1 ? 2 * sqrt(1 - r^2) : 0.\n",
    "proj2 = (r, ϕ, x, y, r0) -> r0 * proj1(r/r0 - (x * cos(ϕ) + y * sin(ϕ))/r0)\n",
    "sino = proj2.(r, ϕ', 5, 1, 3)\n",
    "jimsino = (sino, title) -> jim(\n",
    "    r, ϕ, sino; title, aspect_ratio=:none,\n",
    "    xlabel = \"r\", ylabel = \"ϕ\", ylims=(0,π), yticks=([0, π], [\"0\", \"π\"]),\n",
    "    yflip=false, xticks = [-10, 0, 2, 5, 8, 10],\n",
    "    clim = (0, 6),\n",
    ")\n",
    "jimsino(sino, \"Sinogram for one disk\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "As this figure illustrates,\n",
    "the Radon transform of a unit disk\n",
    "has a somewhat sinusoidal shape.\n",
    "Indeed every point in the $(x,y)$ plane\n",
    "traces out a distinct sinusoid\n",
    "in the sinogram.\n",
    "\n",
    "The mapping from the object $f(x,y)$\n",
    "to data like the sinogram $p(r,\\phi)$\n",
    "is called the \"forward problem.\"\n",
    "\n",
    "\n",
    "## Image reconstruction\n",
    "\n",
    "[Image reconstruction](http://en.wikipedia.org/wiki/Image_reconstruction)\n",
    "is the process of solving the\n",
    "[inverse problem](https://en.wikipedia.org/wiki/Inverse_problem)\n",
    "of recovering an estimate $\\hat{f}(x,y)$\n",
    "of the object $f(x,y)$\n",
    "from a measured sinogram,\n",
    "i.e.,\n",
    "from (usually noisy) samples of $p(r,\\phi)$.\n",
    "\n",
    "\n",
    "## FBP\n",
    "\n",
    "A simple image reconstruction method\n",
    "is called the\n",
    "\"filtered back-projection\" (FBP) approach.\n",
    "It works by filtering each row of the sinogram\n",
    "with a filter,\n",
    "called the\n",
    "[ramp filter](https://en.wikipedia.org/wiki/Radon_transform#Radon_inversion_formula),\n",
    "whose frequency response\n",
    "is roughly $|\\nu|$,\n",
    "where $\\nu$ is the spatial frequency variable\n",
    "(units cycles / m),\n",
    "followed by a\n",
    "[back-projection](https://en.wikipedia.org/wiki/Radon_transform#Dual_transform)\n",
    "step that is the\n",
    "[adjoint](https://en.wikipedia.org/wiki/Hermitian_adjoint)\n",
    "of the Radon transform.\n",
    "\n",
    "\n",
    "### FBP example\n",
    "\n",
    "Here is an illustration\n",
    "of using the `fbp` method\n",
    "in this package\n",
    "to perform image reconstruction\n",
    "from the preceding sinogram."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "image = fbp(sino; dr)\n",
    "(nx,ny) = size(image)\n",
    "dx = dr # default\n",
    "x = ((-(nx-1)/2):((nx-1)/2)) * dr\n",
    "y = x\n",
    "jim(x, y, image, \"FBP image\",\n",
    "    xtick=[-10, 0, 2, 5, 8, 10],\n",
    "    ytick=[-10, 0, -2, 1, 4, 10],\n",
    ")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The FBP reconstructed image\n",
    "looks pretty similar\n",
    "to a disk of radius $3$\n",
    "centered at $(5,1)$\n",
    "as expected.\n",
    "However,\n",
    "there are quite a few ripples;\n",
    "these are\n",
    "[aliasing](https://en.wikipedia.org/wiki/Aliasing)\n",
    "artifacts\n",
    "due to the finite sampling\n",
    "in $r$ and $\\phi$.\n",
    "\n",
    "This is example is what is called\n",
    "2D parallel-beam tomography,\n",
    "because for the angles $\\phi$ are equally spaced\n",
    "and for each angle\n",
    "the radial samples $r$ are also equally spaced.\n",
    "This package includes\n",
    "FBP reconstruction methods\n",
    "for several other sinogram geometries,\n",
    "including the well-known fan beam geometries\n",
    "and the specialized Mojette geometry.\n",
    "\n",
    "\n",
    "## Noise effects on FBP\n",
    "\n",
    "Simulating the effects of measurement noise in sinogram\n",
    "leads to even worse FBP results.\n",
    "\n",
    "First note that a practical imaging system\n",
    "has a finite field of view (FOV):"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "rmax = maximum(r)\n",
    "fovmask = @. sqrt(abs2(x) + abs2(y)') ≤ rmax\n",
    "jim(x, y, fovmask, \"FOV mask\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Add noise to the original sinogram:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "noisy_sinogram = sino + 0.1 * randn(size(sino))\n",
    "jimsino(noisy_sinogram, \"Noisy sinogram\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Apply FBP to the noisy sinogram:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "noisy_fbp_image = fbp(noisy_sinogram; dr)\n",
    "noisy_fbp_image .*= fovmask # apply FOV mask\n",
    "jim(x, y, noisy_fbp_image, \"Noisy FBP image\"; clim=(0,1))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The methods in this package\n",
    "(WIP)\n",
    "and related methods in the\n",
    "[JuliaImageRecon](https://github.com/JuliaImageRecon)\n",
    "suite\n",
    "are designed\n",
    "to provide better reconstructions\n",
    "than the simple FBP method.\n",
    "In particular,\n",
    "model-based image reconstruction (MBIR) methods\n",
    "and methods that use suitable machine-learning approaches\n",
    "can improve image quality significantly.\n",
    "See this\n",
    "[2020 survey paper](http://doi.org/10.1109/JPROC.2019.2936204)."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Reproducibility"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "This page was generated with the following version of Julia:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "io = IOBuffer(); versioninfo(io); split(String(take!(io)), '\\n')"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "And with the following package versions"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "import Pkg; Pkg.status()"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "---\n",
    "\n",
    "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
   ],
   "metadata": {}
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.7.3"
  },
  "kernelspec": {
   "name": "julia-1.7",
   "display_name": "Julia 1.7.3",
   "language": "julia"
  }
 },
 "nbformat": 4
}