{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# LinearMapsAA overview\n",
    "\n",
    "This page illustrates the Julia package\n",
    "[`LinearMapsAA`](https://github.com/JeffFessler/LinearMapsAA.jl)."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Setup"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Packages needed here."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using LinearMapsAA\n",
    "using ImagePhantoms: shepp_logan, SheppLoganToft\n",
    "using MIRTjim: jim, prompt\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",
    "Many computational imaging methods use system models\n",
    "that are too large to store explicitly\n",
    "as dense matrices,\n",
    "but nevertheless\n",
    "are represented mathematically\n",
    "by a linear mapping `A`.\n",
    "\n",
    "Often that linear map is thought of as a matrix,\n",
    "but in imaging problems\n",
    "it often is more convenient\n",
    "to think of it as a more general linear operator.\n",
    "\n",
    "The `LinearMapsAA` package\n",
    "can represent both \"matrix\" versions\n",
    "and \"operator\" versions\n",
    "of linear mappings.\n",
    "This page illustrates both versions\n",
    "in the context of single-image\n",
    "[super-resolution](https://en.wikipedia.org/wiki/Super-resolution_imaging)\n",
    "imaging,\n",
    "where the operator `A` maps a `M × N` image\n",
    "into a coarser sampled image of size `M÷2 × N÷2`.\n",
    "\n",
    "Here the operator `A` is akin to down-sampling,\n",
    "except, rather than simple decimation,\n",
    "each coarse-resolution pixel\n",
    "is the average of a 2 × 2 block of pixels in the fine-resolution image."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## System operator (linear mapping) for down-sampling"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here is the \"forward\" function needed to model 2× down-sampling:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "down1 = (x) -> (x[1:2:end,:] + x[2:2:end,:])/2 # 1D down-sampling by 2×\n",
    "down2 = (x) -> down1(down1(x)')'; # 2D down-sampling by factor of 2×"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The `down2` function is a (bounded) linear operator\n",
    "and here is its adjoint:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "down2_adj(y::AbstractMatrix{<:Number}) = kron(y, fill(0.25, (2,2)));"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Mathematically, and adjoint is a generalization\n",
    "of the (Hermitian) transpose of a matrix.\n",
    "For a (bounded) linear mapping `A` between\n",
    "inner product space X\n",
    "with inner product <.,.>_X\n",
    "and inner product space Y\n",
    "with inner product <.,.>_Y,\n",
    "the adjoint of `A`, denoted `A'`,\n",
    "is the unique bound linear mapping\n",
    "that satisfies\n",
    "<A x, y>_Y = <x, A' y>_X\n",
    "for all x ∈ X and y ∈ Y.\n",
    "One can verify that the `down2_adj` function\n",
    "satisfies that equality\n",
    "for the usual inner product\n",
    "on the space of `M × N` images."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## LinearMap as an operator for super-resolution\n",
    "\n",
    "We now pick a specific image size\n",
    "and define the linear mapping\n",
    "using the two functions above:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "nx, ny = 200, 256\n",
    "A = LinearMapAA(down2, down2_adj, ((nx÷2)*(ny÷2), nx*ny);\n",
    "    idim = (nx,ny), odim = (nx,ny) .÷ 2)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The `idim` argument specifies\n",
    "that the input image is of size `nx × ny`\n",
    "and\n",
    "the `odim` argument specifies\n",
    "that the output image is of size `nx÷2 × ny÷2`.\n",
    "This means that when we invoke\n",
    "`y = A * x`\n",
    "the input `x` must be a 2D array of size `nx × ny`\n",
    "(not a 1D vector!)\n",
    "and the output `y` will have size `nx÷2 × ny÷2`.\n",
    "This behavior is a generalization\n",
    "of what one might expect\n",
    "from a conventional matrix-vector expression,\n",
    "but is quite appropriate and convenient\n",
    "for imaging problems.\n",
    "\n",
    "Here is an illustration.\n",
    "We start with a 2D test image."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "image = shepp_logan(ny, SheppLoganToft())[(ny-nx)÷2 .+ (1:nx),:]\n",
    "jim(image, \"SheppLogan\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Apply the operator `A` to this image to down-sample it:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "down = A * image\n",
    "jim(down, title=\"Down-sampled image\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Apply the adjoint of `A` to that result to \"up-sample\" it:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "up = A' * down\n",
    "jim(up, title=\"Adjoint: A' * y\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "That up-sampled image does not have the same range of values\n",
    "as the original image because `A'` is an adjoint, not an inverse!"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## AbstractMatrix version\n",
    "\n",
    "Some users may prefer that the operator `A` behave more like a matrix.\n",
    "We can implement approach from the same ingredients\n",
    "by using `vec` and `reshape` judiciously.\n",
    "The code is less elegant,\n",
    "but similarly efficient\n",
    "because `vec` and `reshape` are non-allocating operations."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "B = LinearMapAA(\n",
    "        x -> vec(down2(reshape(x,nx,ny))),\n",
    "        y -> vec(down2_adj(reshape(y,Int(nx/2),Int(ny/2)))),\n",
    "        ((nx÷2)*(ny÷2), nx*ny),\n",
    "    )"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "To apply this version to our `image`\n",
    "we must first vectorize it\n",
    "because the expected input is a vector in this case.\n",
    "And then we have to reshape the vector output\n",
    "as a 2D array to look at it.\n",
    "(This is why the operator version with `idim` and `odim` is preferable.)"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "y = B * vec(image) # This would fail here without the `vec`!\n",
    "jim(reshape(y, nx÷2, ny÷2)) # Annoying reshape needed here!"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Even though we write `y = A * x` above,\n",
    "one must remember that internally `A` is not stored as a dense matrix.\n",
    "It is simply a special variable type\n",
    "that stores the forward function `down2` and the adjoint function `down2_adj`,\n",
    "along with a few other values like `nx,ny`,\n",
    "and applies those functions when needed.\n",
    "Nevertheless,\n",
    "we can examine elements of `A` and `B`\n",
    "just like one would with any matrix,\n",
    "at least for small enough examples to fit in memory."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Examine `A` and `A'`"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "nx, ny = 8,6\n",
    "idim = (nx,ny)\n",
    "odim = (nx,ny) .÷ 2\n",
    "A = LinearMapAA(down2, down2_adj, ((nx÷2)*(ny÷2), nx*ny); idim, odim)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here is `A` shown as a Matrix:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "jim(Matrix(A)', \"A\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here is `A'` shown as a Matrix:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "jim(Matrix(A')', \"A'\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "When defining the adjoint function of a linear mapping,\n",
    "it is very important to verify\n",
    "that it is correct (truly the adjoint).\n",
    "\n",
    "For a small problem we simply use the following test:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "@assert Matrix(A)' == Matrix(A')"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "For some applications we must check approximate equality like this:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "@assert Matrix(A)' ≈ Matrix(A')"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here is a statistical test that is suitable for large operators.\n",
    "Often one would repeat this test several times:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "T = eltype(A)\n",
    "x = randn(T, idim)\n",
    "y = randn(T, odim)\n",
    "\n",
    "@assert sum((A*x) .* y) ≈ sum(x .* (A'*y)) # <Ax,y> = <x,A'y>"
   ],
   "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
}