{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# CBCT FDK"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "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).\n",
    "\n",
    "This page was generated from a single Julia file:\n",
    "[07-fdk.jl](https://github.com/JuliaImageRecon/Sinograms.jl/blob/main/docs/lit/examples/07-fdk.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 # these 2 must precede Sinograms for Requires to work!\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 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": [
    "### 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",
    "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",
    "cg = 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(cg, 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(cg), ob)\n",
    "jim(cg.s, cg.t, 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(cg, ig; window = Window(Hamming(), 1.0))\n",
    "fdk_arc = fdk(plan, proj_arc)\n",
    "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",
    "jim(axes(ig), err_arc, \"Error image (arc)\"; clim = (-1,1) .* (0.05μ))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Repeat with flat detector geometry"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "cg = CtFanFlat( ; p...)\n",
    "proj_flat = radon(rays(cg), ob)\n",
    "jim(cg.s, cg.t, proj_flat ;\n",
    "    title=\"Shepp-Logan projections (flat)\", xlabel=\"s\", ylabel=\"t\")\n",
    "\n",
    "plan = plan_fbp(cg, ig; window = Window(Hamming(), 1.0))\n",
    "fdk_flat = fdk(plan, proj_flat)\n",
    "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",
    "jim(axes(ig), err_flat, \"Error image (flat)\"; clim = (-1,1) .* (0.05μ))"
   ],
   "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": [
    "---\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.8.3"
  },
  "kernelspec": {
   "name": "julia-1.8",
   "display_name": "Julia 1.8.3",
   "language": "julia"
  }
 },
 "nbformat": 4
}