{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# SPECTrecon 2D use\n",
    "\n",
    "This page describes how to perform 2D SPECT forward and back-projection\n",
    "using 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: SPECTplan, psf_gauss\n",
    "using SPECTrecon: project, project!, backproject, backproject!\n",
    "using MIRTjim: jim, prompt\n",
    "using ImagePhantoms: shepp_logan, SheppLoganEmis\n",
    "using LinearAlgebra: mul!\n",
    "using LinearMapsAA: LinearMapAA\n",
    "using Plots: plot, default; default(markerstrokecolor=:auto)"
   ],
   "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",
    "Real SPECT systems are inherently 3D imaging systems,\n",
    "but for the purpose of prototyping algorithms\n",
    "it can be useful to work with 2D simulations.\n",
    "\n",
    "Currently, \"2D\" here means a 3D array with `nz=1`,\n",
    "i.e., a single slice.\n",
    "The key to working with a single slice\n",
    "is that the package allows the PSFs\n",
    "to have rectangular support `px × pz`\n",
    "where `pz = 1`, i.e., no blur along the axial (z) direction.\n",
    "\n",
    "\n",
    "## Example\n",
    "\n",
    "Start with a simple 2D digital phantom."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "T = Float32\n",
    "nx,ny,nz = 128,128,1\n",
    "xtrue = T.(shepp_logan(nx, SheppLoganEmis()))\n",
    "xtrue = reshape(xtrue, nx, ny, 1) # 3D array with nz=1\n",
    "jim(xtrue, \"xtrue: SheppLoganEmis with size $(size(xtrue))\")"
   ],
   "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,pz = 11,1 # pz=1 is crucial for 2D work\n",
    "psf1 = psf_gauss( ; ny, px, pz, fwhm_start = 1, fwhm_end = 4) # (px,pz,ny)\n",
    "tmp = reshape(psf1, px, ny) / maximum(psf1) # (px,ny)\n",
    "hx = (px-1)÷2\n",
    "plot(-hx:hx, tmp[:,[1:9:end-10;end]], markershape=:o, label=\"\",\n",
    "    title = \"Depth-dependent PSF profiles\",\n",
    "    xtick = [-hx, -2, 0, 2, hx], # (-1:1) .* ((px-1)÷2),\n",
    "    ytick = [0; round.(tmp[hx+1,end] * [0.5,1], digits=2); 0.5; 1],\n",
    ")\n",
    "prompt()"
   ],
   "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 = 80\n",
    "psfs = repeat(psf1, 1, 1, 1, nview)\n",
    "size(psfs)"
   ],
   "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) # [nx,1,nview]\n",
    "sino = reshape(views, nx, nview)\n",
    "size(sino)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Display the calculated (i.e., simulated) projection views"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "jim(sino, \"Sinogram\")"
   ],
   "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": [
    "sino1 = zeros(T, nx, nview)\n",
    "sino1[nx÷2, nview÷5] = 1\n",
    "sino1[nx÷2, 1] = 1\n",
    "sino1 = reshape(sino1, nx, nz, nview)\n",
    "back1 = backproject(sino1, mumap, psfs, dy)\n",
    "jim(back1, \"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(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 `LinearMapsAA`\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": [
    "## Gram matrix impulse response"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "points = zeros(T, nx, ny, nz)\n",
    "points[nx÷2,ny÷2,1] = 1\n",
    "points[3nx÷4,ny÷4,1] = 1\n",
    "\n",
    "impulse = A' * (A * points)\n",
    "jim(impulse, \"Impulse response of A'A\")"
   ],
   "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
}