{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# CBCT FDK\n",
    "\n",
    "This page describes image reconstruction\n",
    "for cone-beam computed tomography (CBCT)\n",
    "with both arc and flat detectors\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 Plots: plot, gui\n",
    "using Unitful: cm\n",
    "using Sinograms: CtFanArc, CtFanFlat # CtPar\n",
    "using Sinograms: rays, plan_fbp, Window, Hamming, fdk, ct_geom_plot3\n",
    "using ImageGeoms: ImageGeom, MaskCircle, fovs\n",
    "using ImagePhantoms: ellipsoid_parameters, ellipsoid\n",
    "using ImagePhantoms: radon, phantom\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": [
    "## CBCT projections of 3D Shepp-Logan phantom\n",
    "\n",
    "For illustration,\n",
    "we start by synthesizing\n",
    "CBCT projections\n",
    "of the 3D 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=(64,62,30), deltas = (1,1,2) .* 0.4cm)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Ellipsoid parameters for 3D Shepp-Logan phantom:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "μ = 0.1 / cm # typical linear attenuation coefficient\n",
    "params = ellipsoid_parameters( ; fovs = fovs(ig), u = (1, 1, μ))\n",
    "ob = ellipsoid(params)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "A narrow grayscale window is needed to see the soft tissue structures"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "clim = (0.95, 1.05) .* μ"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here is the ideal phantom image:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "oversample = 3\n",
    "true_image = phantom(axes(ig)..., ob, oversample)\n",
    "pt = jim(axes(ig), true_image, \"True 3D Shepp-Logan phantom image\"; clim)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Define the system geometry\n",
    "(for some explanation use `?CtGeom`):"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "p = (ns = 130, ds = 0.3cm, nt = 80, dt = 0.4cm, na = 50, dsd = 200cm, dod = 40cm)\n",
    "rg = CtFanArc( ; p...)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Examine the geometry to verify the FOV\n",
    "(this is more interesting when interacting via other Plot backends):"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "ct_geom_plot3(rg, ig)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "prompt()"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "CBCT projections\n",
    "using `Sinogram.rays` and `ImagePhantoms.radon`:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "proj_arc = radon(rays(rg), ob)\n",
    "pa = jim(axes(rg)[1:2], proj_arc ;\n",
    "    title=\"Shepp-Logan projections (arc)\", xlabel=\"s\", ylabel=\"t\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "There is no \"inverse crime\" here\n",
    "because we compute the projection views\n",
    "using the analytical phantom geometry,\n",
    "but then reconstruct\n",
    "on a discrete grid.\n",
    "\n",
    "## Image reconstruction via FBP / FDK\n",
    "\n",
    "We start with a \"plan\",\n",
    "which would save work if we were reconstructing many images.\n",
    "For illustration we include `Hamming` window."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "plan = plan_fbp(rg, ig; window = Window(Hamming(), 1.0))\n",
    "fdk_arc = fdk(plan, proj_arc)\n",
    "par = jim(axes(ig), fdk_arc, \"FDK image (arc)\"; clim)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "err_arc = fdk_arc - true_image\n",
    "elim = (-1,1) .* (0.02μ)\n",
    "pae = jim(axes(ig), err_arc, \"Error image (arc)\"; clim = elim)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Repeat with flat detector geometry"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "rg = CtFanFlat( ; p...)\n",
    "proj_flat = radon(rays(rg), ob)\n",
    "pfp = jim(axes(rg)[1:2], proj_flat ;\n",
    "    title=\"Shepp-Logan projections (flat)\", xlabel=\"s\", ylabel=\"t\")\n",
    "\n",
    "plan = plan_fbp(rg, ig; window = Window(Hamming(), 1.0))\n",
    "fdk_flat = fdk(plan, proj_flat)\n",
    "pfr = jim(axes(ig), fdk_flat, \"FDK image (flat)\"; clim)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "err_flat = fdk_flat - true_image\n",
    "pfe = jim(axes(ig), err_flat, \"Error image (flat)\"; clim = elim)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "As expected for CBCT,\n",
    "the largest errors are in the end slices."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Short scan"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "rg = CtFanFlat(:short ; p..., na = 200)\n",
    "proj_short = radon(rays(rg), ob)\n",
    "psp = jim(axes(rg)[1:2], proj_short ;\n",
    "    title=\"Shepp-Logan projections (flat,short)\", xlabel=\"s\", ylabel=\"t\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "plan_short = plan_fbp(rg, ig; window = Window(Hamming(), 1.0))\n",
    "psw = jim(axes(rg)[[1,3]], plan_short.view_weight[:,1,:];\n",
    "    title = \"View weights (including Parker)\",\n",
    "    xlabel=\"s\", ylabel=\"angle\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "fdk_short = fdk(plan_short, proj_short)\n",
    "psr = jim(axes(ig), fdk_short, \"FDK image (flat,short)\"; clim)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "err_short = fdk_short - true_image\n",
    "pse = jim(axes(ig), err_flat, \"Error image (flat,short)\"; clim = elim)"
   ],
   "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
}