{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# ImageGeoms overview\n",
    "\n",
    "This page explains the Julia package\n",
    "[`ImageGeoms`](https://github.com/JuliaImageRecon/ImageGeoms.jl)."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Setup"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Packages needed here."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using ImageGeoms: ImageGeom, axis, axisf, axesf, downsample, grids, oversample\n",
    "import ImageGeoms # ellipse\n",
    "using AxisArrays\n",
    "using MIRTjim: jim, prompt\n",
    "using Plots: scatter, plot!, default; default(markerstrokecolor=:auto)\n",
    "using Unitful: mm, s\n",
    "using InteractiveUtils: versioninfo"
   ],
   "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": [
    "## Overview\n",
    "\n",
    "When performing tomographic image reconstruction,\n",
    "one must specify the geometry of the grid of image pixels.\n",
    "(In contrast, for image denoising and image deblurring problems,\n",
    "one works with the given discrete image\n",
    "and no physical coordinates are needed.)\n",
    "\n",
    "The key parameters of a grid of image pixels are\n",
    "* the size (dimensions) of the grid, e.g., `128 × 128`,\n",
    "* the spacing of the pixels, e.g., `1mm × 1mm`,\n",
    "* the offset of the pixels relative to the origin, e.g., `(0,0)`\n",
    "\n",
    "The data type `ImageGeom` describes such a geometry,\n",
    "for arbitrary dimensions (2D, 3D, etc.).\n",
    "\n",
    "There are several ways to construct this structure.\n",
    "The default is a `128 × 128` grid\n",
    "with pixel size $\\Delta_X = \\Delta_Y = 1$ (unitless) and zero offset:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "ig = ImageGeom()"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here is a 3D example with non-cubic voxel size:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "ig = ImageGeom( (512,512,128), (1,1,2), (0,0,0) )"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "To avoid remembering the order of the arguments,\n",
    "named keyword pairs are also supported:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "ig = ImageGeom( dims=(512,512,128), deltas=(1,1,2), offsets=(0,0,0) )"
   ],
   "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 for a video (2D+time) with 12 frames per second:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "ig = ImageGeom( dims=(640,480,1000), deltas=(1mm,1mm,(1//12)s) )"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Methods\n",
    "\n",
    "An ImageGeom object has quite a few methods;\n",
    "`axes` and `axis` are especially useful:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "ig = ImageGeom( dims=(7,8), deltas=(3,2), offsets=(0,0.5) )\n",
    "axis(ig, 2)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "or an axis of length `n` with spacing `Δ` (possibly with units)\n",
    "and (always unitless but possibly non-integer) `offset` the axis\n",
    "is a subtype of `AbstractRange` of the form\n",
    "`( (0:n-1) .- ((n - 1)/2 + offset) ) * Δ`\n",
    "\n",
    "These axes are useful for plotting:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "ig = ImageGeom( dims=(10,8), deltas=(1mm,1mm), offsets=(0.5,0.5) )"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "_ticks(x, off) = [x[1]; # hopefully helpful tick marks\n",
    "    iseven(length(x)) && iszero(off) ?\n",
    "        oneunit(eltype(x)) * [-0.5,0.5] : zero(eltype(x)); x[end]]\n",
    "showgrid = (ig) -> begin # x,y grid locations of pixel centers\n",
    "    x = axis(ig, 1)\n",
    "    y = axis(ig, 2)\n",
    "    (xg, yg) = grids(ig)\n",
    "    scatter(xg, yg; label=\"\", xlabel=\"x\", ylabel=\"y\",\n",
    "        xlims = extrema(x) .+ (ig.deltas[1] * 0.5) .* (-1,1),\n",
    "        xticks = _ticks(x, ig.offsets[1]),\n",
    "        ylims = extrema(y) .+ (ig.deltas[2] * 0.5) .* (-1,1),\n",
    "        yticks = _ticks(y, ig.offsets[2]),\n",
    "        widen = true, aspect_ratio = 1, title = \"offsets $(ig.offsets)\")\n",
    "end\n",
    "showgrid(ig)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Axes labels can have units"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "prompt()"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Offsets (unitless translation of grid)\n",
    "\n",
    "The default `offsets` are zeros,\n",
    "corresponding to symmetric sampling around origin:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "ig = ImageGeom( dims=(12,10), deltas=(1mm,1mm) )\n",
    "p = showgrid(ig)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "prompt();"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "That default for `offsets` is natural for tomography\n",
    "when considering finite pixel size:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "square = (x,y,Δ,p) -> plot!(p, label=\"\", color=:black,\n",
    "    x .+ Δ[1] * ([0,1,1,0,0] .- 0.5),\n",
    "    y .+ Δ[2] * ([0,0,1,1,0] .- 0.5),\n",
    ")\n",
    "square2 = (x,y) -> square(x, y, ig.deltas, p)\n",
    "square2.(grids(ig)...)\n",
    "plot!(p)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "prompt();"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "In that default geometry, the center `(0,0)` of the image\n",
    "is at a corner of the middle 4 pixels (for even image sizes).\n",
    "That default is typical for tomographic imaging (e.g., CT, PET, SPECT).\n",
    "One must be careful when using operations like `imrotate` or `fft`."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Odd dimensions\n",
    "If an image axis has an odd dimension,\n",
    "then each middle pixel along that axis\n",
    "is centered at 0,\n",
    "for the default `offset=0`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "igo = ImageGeom( dims=(7,6) )\n",
    "po = showgrid(igo)\n",
    "square2 = (x,y) -> square(x, y, igo.deltas, po)\n",
    "square2.(grids(igo)...)\n",
    "plot!(po)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "prompt();"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## AxisArrays"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "There is a natural connection between `ImageGeom` and `AxisArrays`.\n",
    "Note the automatic labeling of units (when relevant) on all axes by\n",
    "[MIRTjim.jim](https://github.com/JeffFessler/MIRTjim.jl)."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "ig = ImageGeom( dims=(60,48), deltas=(1.5mm,1mm) )\n",
    "za = AxisArray( ImageGeoms.ellipse(ig) * 10/mm ; x=axis(ig,1), y=axis(ig,2) )\n",
    "jim(za, \"AxisArray example\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "prompt();"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Resizing\n",
    "\n",
    "Often we have a target grid in mind but want coarser sampling for debugging.\n",
    "The `downsample` method is useful for this."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "ig = ImageGeom( dims = (512,512), deltas = (500mm,500mm) ./ 512 )\n",
    "ig_down = downsample(ig, 4)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Other times we want to avoid an \"inverse crime\" by using finer sampling\n",
    "to simulate data; use `oversample` for this."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "ig_over = oversample(ig, (2,2))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Frequency domain axes\n",
    "For related packages like\n",
    "[`ImagePhantoms`](https://github.com/JuliaImageRecon/ImagePhantoms.jl),\n",
    "there are also `axisf` and `axesf` methods\n",
    "that return the frequency axes\n",
    "associated with an FFT-based approximation\n",
    "to a Fourier transform,\n",
    "where $Δ_f Δ_X = 1/N.$"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "ig = ImageGeom( dims=(4,5), deltas=(1mm,1mm) )\n",
    "axesf(ig)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "axisf(ig, 1)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "axisf(ig, 2)"
   ],
   "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
}