{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# SPECTrecon rotation\n",
    "\n",
    "This page explains the image rotation portion of the Julia package\n",
    "[`SPECTrecon.jl`](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_rotate, imrotate!, imrotate_adj!\n",
    "using MIRTjim: jim, prompt\n",
    "using Plots: scatter, scatter!, plot!, default\n",
    "default(markerstrokecolor=:auto, markersize=3)"
   ],
   "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",
    "The first step in SPECT image forward projection\n",
    "is to rotate each slice of a 3D image volume\n",
    "to the appropriate view angle.\n",
    "In principle\n",
    "one could use any of numerous candidate interpolation methods\n",
    "for this task.\n",
    "However, because emission images are nonnegative\n",
    "and maximum likelihood methods\n",
    "for SPECT image reconstruction\n",
    "exploit that nonnegativity,\n",
    "it is desirable to use interpolators\n",
    "that preserve nonnegativity.\n",
    "This constraint rules out quadratic and higher B-splines,\n",
    "including the otherwise attractive cubic B-spline methods.\n",
    "On the other hand, nearest-neighbor interpolation\n",
    "(equivalent to 0th-order B-splines)\n",
    "does not provide adequate image quality.\n",
    "This leaves 1st-order interpolation methods\n",
    "as the most viable options.\n",
    "\n",
    "This package supports two 1st-order linear interpolators:\n",
    "* 2D bilinear interpolation,\n",
    "* a 3-pass rotation method based on 1D linear interpolation.\n",
    "\n",
    "Because image rotation is done repeatedly\n",
    "(for every slice of both the emission image and the attenuation map,\n",
    "for both projection and back-projection,\n",
    "and for multiple iterations),\n",
    "it is important for efficiency\n",
    "to use mutating methods\n",
    "rather than to repeatedly make heap allocations.\n",
    "\n",
    "Following other libraries like\n",
    "[FFTW.jl](https://github.com/JuliaMath/FFTW.jl),\n",
    "the rotation operations herein start with a `plan`\n",
    "where work arrays are preallocated\n",
    "for subsequent use.\n",
    "The `plan` is a `Vector` of `PlanRotate` objects:\n",
    "one for each thread.\n",
    "(Parallelism is across slices for a 3D image volume.)\n",
    "The number of threads defaults to `Threads.nthreads()`.\n",
    "\n",
    "\n",
    "## Example\n",
    "\n",
    "Start with a 3D image volume (just 2 slices here for simplicity)."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "T = Float32 # work with single precision to save memory\n",
    "image = zeros(T, 64, 64, 2)\n",
    "image[30:50,20:30,1] .= 1\n",
    "image[25:28,20:40,2] .= 1\n",
    "jim(image, \"Original image\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now plan the rotation\n",
    "by specifying\n",
    "* the image size `nx` (it must be square, so `ny=nx` implicitly)\n",
    "* the `Type` used for the work arrays."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "plan2 = plan_rotate(size(image, 1); T)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here are the internals for the plan for the first thread:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "plan2[1]"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "With this `plan` preallocated, now we can rotate the image volume,\n",
    "specifying the rotation angle in radians:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "result2 = similar(image) # allocate memory for the result\n",
    "imrotate!(result2, image, π/6, plan2) # mutates the first argument\n",
    "jim(result2, \"Rotated image by π/6 (2D bilinear)\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The default, shown above, uses 2D bilinear interpolation for rotation.\n",
    "That default is the recommended approach because it is faster."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here is the 3-pass 1D interpolation approach,\n",
    "included mainly for checking consistency\n",
    "with the historical ASPIRE approach used in Matlab version of MIRT."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "plan1 = plan_rotate(size(image, 1); T, method=:one)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here are the plan internals for the first thread:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "plan1[1]"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The results of rotation using 3-pass 1D interpolation look quite similar:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "result1 = similar(image)\n",
    "imrotate!(result1, image, π/6, plan1)\n",
    "jim(result1, \"Rotated image by π/6 (3-pass 1D)\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here are the difference images for comparison."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "jim(result1 - result2, \"Difference images\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Adjoint\n",
    "\n",
    "To ensure adjoint consistency between SPECT forward- and back-projection,\n",
    "there is also an adjoint routine:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "adj2 = similar(result2)\n",
    "imrotate_adj!(adj2, result2, π/6, plan2)\n",
    "jim(adj2, \"Adjoint image rotation (2D)\")\n",
    "\n",
    "adj1 = similar(result1)\n",
    "imrotate_adj!(adj1, result1, π/6, plan1)\n",
    "jim(adj1, \"Adjoint image rotation (3-pass 1D)\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The adjoint is *not* the same as the inverse\n",
    "so one does not expect the output here to match the original image!"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## LinearMap\n",
    "\n",
    "One can form a linear map corresponding to image rotation using `LinearMapAA`.\n",
    "An operator like this may be useful\n",
    "as part of a motion-compensated image reconstruction method."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using LinearMapsAA: LinearMapAA\n",
    "\n",
    "nx = 20 # small size for illustration\n",
    "r1 = plan_rotate(nx; T, nthread = 1, method=:two)[1]\n",
    "r2 = plan_rotate(nx; T, nthread = 1, method=:one)[1]\n",
    "idim = (nx,nx)\n",
    "odim = (nx,nx)\n",
    "forw! = (y,x) -> imrotate!(y, x, π/6, r1)\n",
    "back! = (x,y) -> imrotate_adj!(x, y, π/6, r1)\n",
    "A1 = LinearMapAA(forw!, back!, (prod(odim),prod(idim)); T, odim, idim)\n",
    "forw! = (y,x) -> imrotate!(y, x, π/6, r2)\n",
    "back! = (x,y) -> imrotate_adj!(x, y, π/6, r2)\n",
    "A2 = LinearMapAA(forw!, back!, (prod(odim),prod(idim)); T, odim, idim)\n",
    "\n",
    "Afull1 = Matrix(A1)\n",
    "Aadj1 = Matrix(A1')\n",
    "Afull2 = Matrix(A2)\n",
    "Aadj2 = Matrix(A2')\n",
    "jim(cat(dims=3, Afull1', Aadj1', Afull2', Aadj2'), \"Linear map for 2D rotation and its adjoint\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The following verify adjoint consistency:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "@assert Afull1' ≈ Aadj1\n",
    "@assert Afull2' ≈ Aadj2"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Applying this linear map to a 2D or 3D image performs rotation:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "image2 = zeros(nx,nx); image2[4:6, 5:13] .= 1\n",
    "jim(cat(dims=3, image2, A2 * image2), \"Rotation via linear map: 2D\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here is 3D too.\n",
    "The `A2 * image3` here uses the advanced \"operator\" feature of\n",
    "[LinearMapsAA.jl](https://github.com/JeffFessler/LinearMapsAA.jl)."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "image3 = cat(dims=3, image2, image2')\n",
    "jim(cat(dims=4, image3, A2 * image3), \"Rotation via linear map: 3D\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Examine row and column sums of linear map"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "scatter(xlabel=\"pixel index\", ylabel=\"row or col sum\")\n",
    "scatter!(vec(sum(Afull1, dims=1)), label=\"dim1 sum1\", marker=:x)\n",
    "scatter!(vec(sum(Afull1, dims=2)), label=\"dim2 sum1\", marker=:square)\n",
    "scatter!(vec(sum(Afull2, dims=1)), label=\"dim1 sum2\")\n",
    "scatter!(vec(sum(Afull2, dims=2)), label=\"dim2 sum2\")"
   ],
   "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
}