{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Parallel-beam tomography\n",
    "\n",
    "This page describes parallel-beam tomographic image reconstruction\n",
    "using 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 Sinograms: SinoPar, rays, plan_fbp, fbp, fbp_sino_filter\n",
    "using ImageGeoms: ImageGeom, fovs, MaskCircle\n",
    "using ImagePhantoms: SheppLogan, shepp_logan, radon, phantom\n",
    "using Unitful: mm\n",
    "using MIRTjim: jim, prompt"
   ],
   "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": [
    "## Parallel-beam sinogram of Shepp-Logan phantom\n",
    "\n",
    "For illustration,\n",
    "we start by synthesizing\n",
    "a parallel-beam sinogram\n",
    "of the Shepp-Logan phantom.\n",
    "\n",
    "For completeness,\n",
    "we use units (from Unitful),\n",
    "but units are optional."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Use `ImageGeom` to define the image geometry."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "ig = ImageGeom(MaskCircle(); dims=(128,126), deltas = (2mm,2mm) )"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Use `SinoPar` to define the sinogram geometry."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "rg = SinoPar( ; nb = 130, d = 2mm, na = 100)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Ellipse parameters for Shepp-Logan phantom:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "μ = 0.01 / mm # typical linear attenuation coefficient\n",
    "ob = shepp_logan(SheppLogan(); fovs = fovs(ig), u = (1, 1, μ))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Radon transform of Shepp-Logan phantom:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "sino = radon(rays(rg), ob)\n",
    "jim(axes(rg), sino; title=\"Shepp-Logan sinogram\", xlabel=\"r\", ylabel=\"ϕ\")\n",
    "\n",
    "# Image reconstruction via FBP"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here we start with a \"plan\",\n",
    "which would save work if we were reconstructing many images."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "plan = plan_fbp(rg, ig)\n",
    "fbp_image = fbp(plan, sino)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "A narrow color window is needed to see the soft tissue structures:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "clim = (0.9, 1.1) .* μ\n",
    "jim(axes(ig), fbp_image, \"FBP image\"; clim)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "For comparison, here is the ideal phantom image"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "true_image = phantom(axes(ig)..., ob, 2)\n",
    "jim(axes(ig)..., true_image, \"True phantom image\"; clim)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "For fun, here is the filtered sinogram:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "sino_filt = fbp_sino_filter(sino, plan.filter)\n",
    "jim(axes(rg), sino_filt; title=\"Filtered sinogram\", xlabel=\"r\", ylabel=\"ϕ\")"
   ],
   "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
}