{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# 2D sinogram geometry\n",
    "\n",
    "This page describes the 2D sinogram geometries\n",
    "available in the Julia package\n",
    "[`Sinograms.jl`](https://github.com/JuliaImageRecon/Sinograms.jl)."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Setup"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Packages needed here."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Unitful: mm, °\n",
    "using Plots # these 2 must precede 'using Sinograms' for sino_plot_rays to work\n",
    "using Sinograms: SinoPar, SinoMoj, SinoFanArc, SinoFanFlat, SinoFan\n",
    "using Sinograms: sino_plot_rays, rays, angles\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 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": [
    "## 2D Sinogram geometries\n",
    "\n",
    "To perform 2D image reconstruction\n",
    "from a 2D sinogram,\n",
    "one must first describe\n",
    "how the rays in the sinogram are sampled.\n",
    "\n",
    "Mathematically,\n",
    "the Radon transform $p(r,\\phi)$ of $f(x,y)$\n",
    "is defined for all $r$ and $ϕ$\n",
    "but practical systems\n",
    "record $p(r,ϕ)$\n",
    "only for certain finite sets of samples\n",
    "of those values.\n",
    "\n",
    "\n",
    "## Parallel-beam geometry\n",
    "\n",
    "The simplest form of sinogram sampling\n",
    "is equally spaced samples\n",
    "in both $r$ and $ϕ$.\n",
    "This sampling is called a parallel-beam geometry.\n",
    "(Very few practical systems have this sampling pattern,\n",
    "but it is an easy place to start.)\n",
    "Both FBP and iterative reconstruction methods\n",
    "need to know which values\n",
    "of $r$ and $ϕ$\n",
    "are sampled.\n",
    "\n",
    "In this package,\n",
    "`SinoPar`\n",
    "is the type that describes\n",
    "a parallel-beam sinogram geometry.\n",
    "\n",
    "The built-in defaults provide helpful reminders about the usage."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "SinoPar()"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "This package supports units via the\n",
    "[Unitful.jl](https://github.com/PainterQubits/Unitful.jl)\n",
    "and using units is recommended\n",
    "(but not required).\n",
    "\n",
    "Here is an example of how to specify\n",
    "a parallel-beam geometry.\n",
    "Everything is a named keyword argument with sensible default values.\n",
    "\n",
    "* `orbit` and `orbit_start` must both be unitless (degrees)\n",
    "  or have the same units (e.g., degrees or radians).\n",
    "* detector spacing `d` (can be unitless or have units e.g., mm).\n",
    "* The projection angles $ϕ$ are equally space and given by\n",
    "  `orbit_start + (0:(nb-1))/nb * orbit`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "rg = SinoPar( ;\n",
    "    nb = 64, # number of radial samples (\"bins\")\n",
    "    na = 30, # number of angular samples\n",
    "    d = 2mm, # detector spacing\n",
    "    offset = 0.25, # quarter detector offset (unitless)\n",
    "    orbit = 180, # angular range (in degrees)\n",
    "    orbit_start = 0, # starting angle (in degrees)\n",
    ")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The struct `rg` has numerous useful properties;\n",
    "type `?SinoGeom` to see the full list.\n",
    "\n",
    "For example,\n",
    "to access the angular samples in degrees\n",
    "type `angles(rg)`"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "angles(rg)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The following function visualizes the sampling pattern."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "sino_plot_rays(rg; ylims=(0,180), yticks=(0:90:180), widen=true, title=\"Parallel\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "prompt()"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Fan-beam CT with an arc detector (3rd generation CT)\n",
    "\n",
    "For a fan-beam geometry,\n",
    "the arguments are the same as for `SinoPar`\n",
    "with the addition of specifying:\n",
    "* `dsd` distance from source to detector\n",
    "* `dod` distance from origin (isocenter) to detector\n",
    "* `dfs` distance from focal point of detector to source\n",
    "  (0 for a 3rd gen arc detector, and `Inf` for a flat detector)\n",
    "* `source_offset` for misaligned systems\n",
    "   where the ray from the source to the detector center\n",
    "   does not intersect the isocenter.\n",
    "   Not fully supported; submit an issue if you need this feature.\n",
    "\n",
    "Here is an example that corresponds to a GE Lightspeed CT system.\n",
    "These numbers are published in\n",
    "[IEEE T-MI Oct. 2006, p.1272-1283](https://doi.org/10.1109/TMI.2006.882141)."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "rg = SinoFanArc( ; nb=888, na=984,\n",
    "    d=1.0239mm, offset=1.25, dsd=949.075mm, dod=408.075mm)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here is a smaller example for plotting the rays."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "rg = SinoFanArc( ; nb=64, na=30,\n",
    "    d=20mm, offset=0.25, dsd=900mm, dod=400mm)\n",
    "sino_plot_rays(rg; ylims=(-50,400), yticks=(0:180:360), widen=true,\n",
    "    title=\"Fan-beam for arc detector\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "prompt()"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Fan-beam CT with a flat detector\n",
    "\n",
    "This geometry is the same as the arc detector\n",
    "except that `dfs=Inf`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "rg = SinoFanFlat( ; nb=64, na=30,\n",
    "    d=20mm, offset=0.25, dsd=900mm, dod=400mm)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here is its sampling plot"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "sino_plot_rays(rg; ylims=(-50,400), yticks=(0:180:360), widen=true,\n",
    "    title=\"Fan-beam for flat detector\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "prompt()"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Mojette sampling\n",
    "This is a specialized sampling geometry\n",
    "that is currently incompletely supported."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "rg = SinoMoj( ; nb=60, na=30)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here is its sampling plot"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "sino_plot_rays(rg; ylims=(0,180), yticks=(0:90:180), widen=true,\n",
    "    title=\"Mojette sampling\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here is a diagram that illustrates how the radial spacing\n",
    "is a function of the projection view angle for the Mojette geometry.\n",
    "The key aspect here\n",
    "is that in each image row\n",
    "the line intersection lengths are identical."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "plot(xlims=(-1,1) .* 3.5, ylims = (-1,1) .* 2.5, xlabel=\"x\", ylabel=\"x\")\n",
    "default(label=\"\")\n",
    "for y in -2:2\n",
    "    plot!([-3, 3], [y, y], color=:black)\n",
    "end\n",
    "for x in -3:3\n",
    "    plot!([x, x], [-2, 2], color=:black)\n",
    "end\n",
    "plot_ray!(r, ϕ) = plot!(\n",
    "    r*cos(ϕ) .+ [-1, 1] * 4 * sin(ϕ),\n",
    "    r*sin(ϕ) .+ [1, -1] * 4 * cos(ϕ),\n",
    "    color = :blue,\n",
    ")\n",
    "ia = 4 # pick an angle\n",
    "i = rays(rg) # iterator\n",
    "rs = [i[1] for i in i] # radial samples\n",
    "ϕs = [i[2] for i in i] # projection angle\n",
    "r = rs[:,ia]\n",
    "r = r[abs.(r) .< 3]\n",
    "ϕ = ϕs[1,ia]\n",
    "plot_ray!.(r, ϕ)\n",
    "plot!(aspect_ratio=1, title = \"Mojette line integrals\")"
   ],
   "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.10.3"
  },
  "kernelspec": {
   "name": "julia-1.10",
   "display_name": "Julia 1.10.3",
   "language": "julia"
  }
 },
 "nbformat": 4
}