{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# SPECTrecon overview\n",
    "\n",
    "This page gives an overview of the Julia package\n",
    "[`SPECTrecon`](https://github.com/JuliaImageRecon/SPECTrecon.jl)."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Setup"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Packages needed here."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using SPECTrecon: plan_psf, psf_gauss, SPECTplan\n",
    "using SPECTrecon: project, project!, backproject, backproject!\n",
    "using MIRTjim: jim, prompt\n",
    "using LinearAlgebra: mul!\n",
    "using LinearMapsAA: LinearMapAA\n",
    "using Plots: scatter, plot!, default; default(markerstrokecolor=:auto)\n",
    "using Plots # @animate, gif\n",
    "using InteractiveUtils: versioninfo"
   ],
   "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": [
    "## Overview\n",
    "\n",
    "To perform SPECT image reconstruction,\n",
    "one must have a model for the imaging system\n",
    "encapsulated in a forward projector and back projector.\n",
    "\n",
    "Mathematically, we write the forward projection process in SPECT\n",
    "as \"y = A * x\" where A is a \"system matrix\"\n",
    "that models the physics of the imaging system\n",
    "(including depth-dependent collimator/detector response\n",
    "and attenuation)\n",
    "and \"x\" is the current guess of the emission image.\n",
    "\n",
    "However, in code we usually cannot literally store \"A\"\n",
    "as dense matrix because it is too large.\n",
    "A typical size in SPECT is that\n",
    "the image `x` is\n",
    "`nx × ny × nz = 128 × 128 × 100`\n",
    "and the array of projection views `y` is\n",
    "`nx × nz × nview = 128 × 100 × 120`.\n",
    "So the system matrix `A` has `1536000 × 1638400` elements\n",
    "which is far to many to store,\n",
    "even accounting for some sparsity.\n",
    "\n",
    "Instead, we write functions called forward projectors\n",
    "that calculate `A * x` \"on the fly\".\n",
    "\n",
    "Similarly, the operation `A' * y`\n",
    "is called \"back projection\",\n",
    "where `A'` denotes the transpose or \"adjoint\" of `A`.\n",
    "\n",
    "\n",
    "## Example\n",
    "\n",
    "To illustrate forward and back projection,\n",
    "it is easiest to start with a simulation example\n",
    "using a digital phantom.\n",
    "The fancy way would be to use a 3D phantom from\n",
    "[ImagePhantoms](https://github.com/JuliaImageRecon/ImagePhantoms.jl),\n",
    "but instead we just use two simple cubes."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "nx,ny,nz = 128,128,80\n",
    "T = Float32\n",
    "xtrue = zeros(T, nx,ny,nz)\n",
    "xtrue[(1nx÷4):(2nx÷3), 1ny÷5:(3ny÷5), 2nz÷6:(3nz÷6)] .= 1\n",
    "xtrue[(2nx÷5):(3nx÷5), 1ny÷5:(2ny÷5), 4nz÷6:(5nz÷6)] .= 2\n",
    "\n",
    "average(x) = sum(x) / length(x)\n",
    "function mid3(x::AbstractArray{T,3}) where {T}\n",
    "    (nx,ny,nz) = size(x)\n",
    "    xy = x[:,:,ceil(Int, nz÷2)]\n",
    "    xz = x[:,ceil(Int,end/2),:]\n",
    "    zy = x[ceil(Int, nx÷2),:,:]'\n",
    "    return [xy xz; zy fill(average(x), nz, nz)]\n",
    "end\n",
    "jim(mid3(xtrue), \"Middle slices of x\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## PSF"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Create a synthetic depth-dependent PSF for a single view"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "px = 11\n",
    "psf1 = psf_gauss( ; ny, px, fwhm_end = 6)\n",
    "jim(psf1, \"PSF for each of $ny planes\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "In general the PSF can vary from view to view\n",
    "due to non-circular detector orbits.\n",
    "For simplicity, here we illustrate the case\n",
    "where the PSF is the same for every view."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "nview = 60\n",
    "psfs = repeat(psf1, 1, 1, 1, nview)\n",
    "size(psfs)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Plan the PSF modeling (see `3-psf.jl`)"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "plan = plan_psf( ; nx, nz, px)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Basic SPECT forward projection\n",
    "\n",
    "Here is a simple illustration\n",
    "of a SPECT forward projection operation.\n",
    "(This is a memory inefficient way to do it!)"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "dy = 4 # transaxial pixel size in mm\n",
    "mumap = zeros(T, size(xtrue)) # μ-map just zero for illustration here\n",
    "views = project(xtrue, mumap, psfs, dy)\n",
    "size(views)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Display the calculated (i.e., simulated) projection views"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "jim(views[:,:,1:4:end], \"Every 4th of $nview projection views\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Basic SPECT back projection\n",
    "\n",
    "This illustrates an \"unfiltered backprojection\"\n",
    "that leads to a very blurry image\n",
    "(again, with a simple memory inefficient usage).\n",
    "\n",
    "First, back-project two \"rays\"\n",
    "to illustrate the depth-dependent PSF."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "tmp = zeros(T, size(views))\n",
    "tmp[nx÷2, nz÷2, nview÷5] = 1\n",
    "tmp[nx÷2, nz÷2, 1] = 1\n",
    "tmp = backproject(tmp, mumap, psfs, dy)\n",
    "jim(mid3(tmp), \"Back-projection of two rays\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now back-project all the views of the phantom."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "back = backproject(views, mumap, psfs, dy)\n",
    "jim(mid3(back), \"Back-projection of ytrue\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Memory efficiency\n",
    "\n",
    "For iterative reconstruction,\n",
    "one must do forward and back-projection repeatedly.\n",
    "It is more efficient to pre-allocate work arrays\n",
    "for those operations,\n",
    "instead of repeatedly making system calls.\n",
    "\n",
    "Here we illustrate the memory efficient versions\n",
    "that are recommended for iterative SPECT reconstruction.\n",
    "\n",
    "First construction the SPECT plan."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "#viewangle = (0:(nview-1)) * 2π # default\n",
    "plan = SPECTplan(mumap, psfs, dy; T)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Mutating version of forward projection:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "tmp = Array{T}(undef, nx, nz, nview)\n",
    "project!(tmp, xtrue, plan)\n",
    "@assert tmp == views"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Mutating version of back-projection:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "tmp = Array{T}(undef, nx, ny, nz)\n",
    "backproject!(tmp, views, plan)\n",
    "@assert tmp ≈ back"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Using `LinearMapAA`\n",
    "\n",
    "Calling `project!` and `backproject!` repeatedly\n",
    "leads to application-specific code.\n",
    "More general code uses the fact that SPECT projection and back-projection\n",
    "are linear operations,\n",
    "so we use `LinearMapAA` to define a \"system matrix\" for these operations."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "forw! = (y,x) -> project!(y, x, plan)\n",
    "back! = (x,y) -> backproject!(x, y, plan)\n",
    "idim = (nx,ny,nz)\n",
    "odim = (nx,nz,nview)\n",
    "A = LinearMapAA(forw!, back!, (prod(odim),prod(idim)); T, odim, idim)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Simple forward and back-projection:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "@assert A * xtrue == views\n",
    "@assert A' * views ≈ back"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Mutating version:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "tmp = Array{T}(undef, nx, nz, nview)\n",
    "mul!(tmp, A, xtrue)\n",
    "@assert tmp == views\n",
    "tmp = Array{T}(undef, nx, ny, nz)\n",
    "mul!(tmp, A', views)\n",
    "@assert tmp ≈ back"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Units\n",
    "\n",
    "The pixel dimensions `deltas` can (and should!) be values with units.\n",
    "\n",
    "Here is an example with units ... (todo)\n",
    "\n",
    "using UnitfulRecipes\n",
    "using Unitful: mm"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Projection view animation"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "anim = @animate for i in 1:nview\n",
    "    ymax = maximum(views)\n",
    "    jim(views[:,:,i],\n",
    "        \"SPECT projection view $i of $nview\",\n",
    "        clim = (0, ymax),\n",
    "    )\n",
    "end\n",
    "gif(anim, \"views.gif\", fps = 8)"
   ],
   "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.11.1"
  },
  "kernelspec": {
   "name": "julia-1.11",
   "display_name": "Julia 1.11.1",
   "language": "julia"
  }
 },
 "nbformat": 4
}