{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# 2D Branchless Distance-driven Projection and Backprojection\n",
    "\n",
    "This page describes the 2D branchless distance-driven projection\n",
    "and backprojection method\n",
    "of\n",
    "[Basu & De Man, 2006](https://doi.org/10.1117/12.659893)\n",
    "for fan-beam geometries with a \"flat\" detector\n",
    "using the Julia package\n",
    "[`Sinograms.jl`](https://github.com/JuliaImageRecon/Sinograms.jl).\n",
    "\n",
    "This page compares the results to the `radon` transform method."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Setup"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Packages needed here."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using ImageGeoms: ImageGeom, fovs, MaskCircle\n",
    "import ImageGeoms # downsample\n",
    "using ImagePhantoms: shepp_logan, SheppLoganToft, radon, phantom\n",
    "using MIRTjim: jim, prompt\n",
    "import Plots\n",
    "using Sinograms: project_bdd, backproject_bdd\n",
    "import Sinograms # downsample, _ar, _dso\n",
    "using Sinograms: SinoFanFlat, rays, plan_fbp, fbp, Window, Hamming, sino_geom_plot!, angles\n",
    "using Unitful: mm"
   ],
   "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": [
    "## Fan-beam sinogram of SheppLoganToft phantom\n",
    "\n",
    "For illustration,\n",
    "we start by synthesizing\n",
    "a fan-beam sinogram\n",
    "of the `SheppLoganToft` phantom.\n",
    "\n",
    "For completeness,\n",
    "we use units (from `Unitful`),\n",
    "but units are optional."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Define the sinogram geometry"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "down = 4 # save time\n",
    "rg = SinoFanFlat( ; nb = 910, d = 1.0239mm, na = 360, dsd = 949mm, dod = 408mm)\n",
    "rg = Sinograms.downsample(rg, down)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Define the image geometry"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "ig = ImageGeom(MaskCircle(); dims=(512,512), deltas = (1mm,1mm))\n",
    "ig = ImageGeoms.downsample(ig, down)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Make a tuple to define imaging geometry for bdd code. todo"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "geo = (DSD = rg.dsd, DS0 = Sinograms._dso(rg), pSize = ig.deltas[1],\n",
    "    dSize = rg.d, nPix = ig.dims[1], nDet = rg.nb, angle = Sinograms._ar(rg))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Examine the geometry to verify the FOV:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "pa = jim(axes(ig), ig.mask; prompt=false)\n",
    "sino_geom_plot!(rg, ig)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "prompt()"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Ellipse parameters for SheppLoganToft phantom:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "μ = 0.01 / mm # typical linear attenuation coefficient\n",
    "ob = shepp_logan(SheppLoganToft(); fovs = fovs(ig), u = (1, 1, μ))\n",
    "testimage = phantom(axes(ig)..., ob) # Create a phantom image\n",
    "pt = jim(axes(ig), testimage) # Show the true phantom image"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Fan-beam sinograms for phantom:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "if !@isdefined(sinogramR)\n",
    "    @time sinogramR = radon(rays(rg), ob)\n",
    "    # Ensure that sinogram is not truncated\n",
    "    @assert all(==(0), sum(abs, sinogramR, dims=2)[[1,end]])\n",
    "end;\n",
    "\n",
    "if !@isdefined(sinogramB)\n",
    "    @time sinogramB = project_bdd(reverse(rot180(testimage'), dims=2), geo)\n",
    "    sinogramB = sinogramB' # todo\n",
    "end;\n",
    "\n",
    "p1 = jim(axes(rg), sinogramB; prompt=false,\n",
    " title=\"bdd_2d sinogram\", xlabel=\"r\", ylabel=\"ϕ\")\n",
    "p2 = jim(axes(rg), sinogramR; prompt=false,\n",
    " title=\"Radon test sinogram\", xlabel=\"r\", ylabel=\"ϕ\")\n",
    "p12 = jim(p1,p2)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Difference of sinogram using `bdd_2d` versus `radon` method."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "pd = jim(axes(rg), sinogramR - sinogramB;\n",
    " title=\"Difference sinogram\", xlabel=\"r\", ylabel=\"ϕ\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Backprojection with `bdd_2d`\n",
    "(The un-filtered back-projected image is not a useful reconstruction.\n",
    "See next section for image reconstruction via FBP.)"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "if !@isdefined(imageB)\n",
    "    @time imageB = backproject_bdd(sinogramB'/1mm, geo) # todo\n",
    "    imageB = rotr90(imageB) # todo\n",
    "end\n",
    "pb = jim(axes(ig), imageB;\n",
    " title=\"bdd_2d backprojection\", xlabel=\"x\", ylabel=\"y\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Image reconstruction via FBP\n",
    "Compare the reconstructed image using the `bdd_2d` sinogram\n",
    "and the `radon` sinogram.\n",
    "\n",
    "Here we start with a \"plan\"\n",
    "that 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",
    "fbp_image_b = fbp(plan, sinogramB)\n",
    "fbp_image_r = fbp(plan, sinogramR);"
   ],
   "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.04, 1.1) .* μ\n",
    "r1 = jim(axes(ig), fbp_image_b, \"FBP image with bdd_2d\"; clim)\n",
    "r2 = jim(axes(ig), fbp_image_r, \"FBP image with radon\"; clim)\n",
    "r12 = jim(r1,r2; size=(1000,300))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "For comparison, here is the difference image."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "pe = jim(axes(ig), fbp_image_b - fbp_image_r, \"Difference image\")"
   ],
   "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
}