diff --git a/.gitignore b/.gitignore index d3e1110..27686f9 100644 --- a/.gitignore +++ b/.gitignore @@ -16,6 +16,7 @@ deps/src/ # Build artifacts for creating documentation generated by the Documenter package docs/build/ docs/site/ +docs/src/generated/ # File generated by Pkg, the package manager, based on a corresponding Project.toml # It records a fixed state of all packages used by the project. As such, it should not be @@ -23,6 +24,7 @@ docs/site/ # environment. Manifest.toml - -# Other +# other .DS_Store +.*.*.swp +.benchmarkci diff --git a/docs/Project.toml b/docs/Project.toml index 17f7d8c..d9c9230 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,3 +1,15 @@ [deps] +AxisArrays = "39de3d68-74b9-583c-8d2d-e117c070f3a9" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" +MIRTjim = "170b2178-6dee-4cb0-8729-b3e8b57834cc" +OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Sinograms = "02a14def-c6e6-4ab0-b2df-0ab64bc8cdd7" +Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" +UnitfulRecipes = "42071c24-d89e-48dd-8a24-8a12d9b8861f" + +[compat] +Documenter = "0.27.3" +Literate = "2" +MIRTjim = "0.11" diff --git a/docs/lit/examples/01-tomography.jl b/docs/lit/examples/01-tomography.jl new file mode 100644 index 0000000..c3eb15a --- /dev/null +++ b/docs/lit/examples/01-tomography.jl @@ -0,0 +1,287 @@ +#--------------------------------------------------------- +# # [Tomography](@id 01-tomography) +#--------------------------------------------------------- + +#= +This page gives an overview of the Julia package +[`Sinograms.jl`](https://github.com/JeffFessler/Sinograms.jl). + +This page was generated from a single Julia file: +[01-tomography.jl](@__REPO_ROOT_URL__/01-tomography.jl). +=# + +#md # In any such Julia documentation, +#md # you can access the source code +#md # using the "Edit on GitHub" link in the top right. + +#md # The corresponding notebook can be viewed in +#md # [nbviewer](http://nbviewer.jupyter.org/) here: +#md # [`01-tomography.ipynb`](@__NBVIEWER_ROOT_URL__/01-tomography.ipynb), +#md # and opened in [binder](https://mybinder.org/) here: +#md # [`01-tomography.ipynb`](@__BINDER_ROOT_URL__/01-tomography.ipynb). + + +# ### Setup + +# Packages needed here. + +using Sinograms: fbp +using MIRTjim: jim, prompt +using InteractiveUtils: versioninfo + + +# The following line is helpful when running this example.jl file as a script; +# this way it will prompt user to hit a key after each figure is displayed. + +isinteractive() ? jim(:prompt, true) : prompt(:draw); + + +#= +### Tomography + +[Tomography](https://en.wikipedia.org/wiki/Tomography) +is the process of imaging cross sections of an object +without actually slicing the object. +There are many forms of tomography; +the description here +focuses on +[X-ray computed tomography (CT scans)](https://en.wikipedia.org/wiki/CT_scan). +(See +[SPECTrecon.jl](https://github.com/JeffFessler/SPECTrecon.jl) +for a Julia package related to +[SPECT](https://en.wikipedia.org/wiki/Single-photon_emission_computed_tomography) +imaging.) + +In an +[X-ray CT imaging system](http://en.wikipedia.org/wiki/File:Ct-internals.jpg), +X-rays emitted from an X-ray source +are transmitted through an object +(e.g., a patient in medical CT) +towards a detector array. +The source and detector +[rotate rapidly](https://www.youtube.com/watch?v=2CWpZKuy-NE) +around the object. +The signal intensities recorded by the detector +are related +to the sizes and densities +of the object materials +between the source and each detector element. + + +## Radon transform + +The mathematical foundation +for 2D X-ray CT imaging +is the +[Radon transform](https://en.wikipedia.org/wiki/Radon_transform), +an +[integral transform](https://en.wikipedia.org/wiki/Integral_transform) +that maps a 2D function +``f(x,y)`` +into the collection of line integrals +through that function. +Here we describe each line +by its angle ``\phi`` +measured counter-clockwise from the ``y`` axis, +and by the radial distance ``r`` of the line from the origin. +The collection of integrals +is called a sinogram. + +Mathematically, +the Radon transform ``p(r,\phi)`` of ``f(x,y)`` is defined by +```math +p(r,\phi) = \int_{-\infty}^{\infty} +f(r \cos \phi - l \sin \phi, +r \sin \phi + l \cos \phi) \, \mathrm{d} l. +``` + +The Radon transform of the object that is 1 inside the unit disk +and 0 elsewhere is given by +```math +p_1(r,\phi) = 2 \sqrt{1 - r^2}, \ \mathrm{ for } \ |r| < 1. +``` + +By the Radon transform's translation and scaling properties, +the Radon transform +of a disk of radius ``r_0`` +centered at ``(x_0,y_0)`` +is given by +```math +p(r,\phi) = r_0 \, p_1 +\left( \frac{r - (x_0 \cos \phi + y_0 \sin \phi)}{r_0}, ϕ \right). +``` + +## Sinogram example + +Here is a display of that Radon transform +for a disk of radius ``3`` +centered at coordinate ``(5,1)``. +Note that maximum value is approximately ``6``, +the length of the longest chord +through a disk of radius ``3``. +=# + +nr = 128 +dr = 20 / nr +r = ((-(nr-1)/2):((nr-1)/2)) * dr # radial sample locations +na = 130 +ϕ = (0:(na-1))/na * π # angular samples +proj1 = r -> abs(r) < 1 ? 2 * sqrt(1 - r^2) : 0. +proj2 = (r, ϕ, x, y, r0) -> r0 * proj1(r/r0 - (x * cos(ϕ) + y * sin(ϕ))/r0) +sino = proj2.(r, ϕ', 5, 1, 3) +jimsino = (sino, title) -> jim( + r, ϕ, sino; title, aspect_ratio=:none, + xlabel = "r", ylabel = "ϕ", ylims=(0,π), yticks=([0, π], ["0", "π"]), + yflip=false, xticks = [-10, 0, 2, 5, 8, 10], + clim = (0, 6), +) +jimsino(sino, "Sinogram for one disk") + + +#= +As this figure illustrates, +the Radon transform of a unit disk +has a somewhat sinusoidal shape. +Indeed every point in the ``(x,y)`` plane +traces out a distinct sinusoid +in the sinogram. + +The mapping from the object ``f(x,y)`` +to data like the sinogram ``p(r,\phi)`` +is called the "forward problem." + + +## Image reconstruction + +[Image reconstruction](http://en.wikipedia.org/wiki/Image_reconstruction) +is the process of solving the +[inverse problem](https://en.wikipedia.org/wiki/Inverse_problem) +of recovering an estimate ``\hat{f}(x,y)`` +of the object ``f(x,y)`` +from a measured sinogram, +i.e., +from (usually noisy) samples of ``p(r,\phi)``. + + +## FBP + +A simple image reconstruction method +is called the +"filtered back-projection" (FBP) approach. +It works by filtering each row of the sinogram +with a filter, +called the +[ramp filter](https://en.wikipedia.org/wiki/Radon_transform#Radon_inversion_formula), +whose frequency response +is roughly ``|\nu|``, +where ``\nu`` is the spatial frequency variable +(units cycles / m), +followed by a +[back-projection](https://en.wikipedia.org/wiki/Radon_transform#Dual_transform) +step that is the +[adjoint](https://en.wikipedia.org/wiki/Hermitian_adjoint) +of the Radon transform. + + +### FBP example + +Here is an illustration +of using the `fbp` method +in this package +to perform image reconstruction +from the preceding sinogram. + +=# + +image = fbp(sino; dr) +(nx,ny) = size(image) +dx = dr # default +x = ((-(nx-1)/2):((nx-1)/2)) * dr +y = x +jim(x, y, image, "FBP image", + xtick=[-10, 0, 2, 5, 8, 10], + ytick=[-10, 0, -2, 1, 4, 10], +) + + +#= + +The FBP reconstructed image +looks pretty similar +to a disk of radius ``3`` +centered at ``(5,1)`` +as expected. +However, +there are quite a few ripples; +these are +[aliasing](https://en.wikipedia.org/wiki/Aliasing) +artifacts +due to the finite sampling +in ``r`` and ``\phi``. + +This is example is what is called +2D parallel-beam tomography, +because for the angles ``\phi`` are equally spaced +and for each angle +the radial samples ``r`` are also equally spaced. +This package includes +FBP reconstruction methods +for several other sinogram geometries, +including the well-known fan beam geometries +and the specialized Mojette geometry. + + +## Noise effects on FBP + +Simulating the effects of measurement noise in sinogram +leads to even worse FBP results. + +First note that a practical imaging system +has a finite field of view (FOV): +=# + +rmax = maximum(r) +fovmask = @. sqrt(abs2(x) + abs2(y)') ≤ rmax +jim(x, y, fovmask, "FOV mask") + + +# Add noise to the original sinogram: +noisy_sinogram = sino + 0.1 * randn(size(sino)) +jimsino(noisy_sinogram, "Noisy sinogram") + +# Apply FBP to the noisy sinogram: +noisy_fbp_image = fbp(noisy_sinogram; dr) +noisy_fbp_image .*= fovmask # apply FOV mask +jim(x, y, noisy_fbp_image, "Noisy FBP image"; clim=(0,1)) + + +#= +The methods in this package +(WIP) +and related methods in the +[JuliaImageRecon](https://github.com/JuliaImageRecon) +suite +are designed +to provide better reconstructions +than the simple FBP method. +In particular, +model-based image reconstruction (MBIR) methods +and methods that use suitable machine-learning approaches +can improve image quality significantly. +See this +[2020 survey paper](http://doi.org/10.1109/JPROC.2019.2936204). +=# + + +# ### Reproducibility + +# This page was generated with the following version of Julia: + +io = IOBuffer() +versioninfo(io) +split(String(take!(io)), '\n') + + +# And with the following package versions + +import Pkg; Pkg.status() diff --git a/docs/make.jl b/docs/make.jl index 97b27ab..020262c 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,29 +1,76 @@ -# push!(LOAD_PATH,"../src/") +execute = isempty(ARGS) || ARGS[1] == "run" using Sinograms using Documenter +using Literate + +# https://juliadocs.github.io/Documenter.jl/stable/man/syntax/#@example-block +ENV["GKSwstype"] = "100" +ENV["GKS_ENCODING"] = "utf-8" + +# generate examples using Literate +lit = joinpath(@__DIR__, "lit") +src = joinpath(@__DIR__, "src") +gen = joinpath(@__DIR__, "src/generated") + +base = "JeffFessler/Sinograms.jl" +repo_root_url = + "https://github.com/$base/blob/main/docs/lit/examples" +nbviewer_root_url = + "https://nbviewer.org/github/$base/tree/gh-pages/dev/generated/examples" +binder_root_url = + "https://mybinder.org/v2/gh/$base/gh-pages?filepath=dev/generated/examples" + DocMeta.setdocmeta!(Sinograms, :DocTestSetup, :(using Sinograms); recursive=true) +for (root, _, files) in walkdir(lit), file in files + splitext(file)[2] == ".jl" || continue # process .jl files only + ipath = joinpath(root, file) + opath = splitdir(replace(ipath, lit => gen))[1] + Literate.markdown(ipath, opath, documenter = execute; # run examples + repo_root_url, nbviewer_root_url, binder_root_url) + Literate.notebook(ipath, opath; execute = false, # no-run notebooks + repo_root_url, nbviewer_root_url, binder_root_url) +end + + +# Documentation structure +ismd(f) = splitext(f)[2] == ".md" +pages(folder) = + [joinpath("generated/", folder, f) for f in readdir(joinpath(gen, folder)) if ismd(f)] + +isci = get(ENV, "CI", nothing) == "true" + +format = Documenter.HTML(; + prettyurls = isci, + edit_link = "main", + canonical = "https://JeffFessler.github.io/Sinograms.jl/stable/", +# assets = String[], +) + makedocs(; modules = [Sinograms], - authors = "Jeff Fessler and contributors", - repo = "https://github.com/JeffFessler/Sinograms.jl/blob/{commit}{path}#{line}", + authors = "Jeff Fessler and contributors", sitename = "Sinograms.jl", - format = Documenter.HTML(; - prettyurls = get(ENV, "CI", "false") == "true", -# canonical = "https://JeffFessler.github.io/Sinograms.jl/stable", -# assets = String[], - ), + format, pages = [ "Home" => "index.md", + "Methods" => "methods.md", + "Examples" => pages("examples") ], ) -deploydocs(; - repo = "github.com/JeffFessler/Sinograms.jl.git", - devbranch = "main", - devurl = "dev", - versions = ["stable" => "v^", "dev" => "dev"] -# push_preview = true, -) +if isci + deploydocs(; + repo = "github.com/JeffFessler/Sinograms.jl.git", + devbranch = "main", + devurl = "dev", + versions = ["stable" => "v^", "dev" => "dev"], + forcepush = true, +# push_preview = true, + # see https://JuliaImageRecon.github.io/ImageGeoms.jl/previews/PR## + ) +else + @warn "may need to: rm -r src/generated/" +end diff --git a/docs/src/index.md b/docs/src/index.md index a43b72e..d52b4bc 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -11,19 +11,18 @@ CurrentModule = Sinograms ## Overview -todo: document Sinogram geometries, etc. - -For use with the -Michigan Image Reconstruction Toolbox -[(MIRT)](https://github.com/JeffFessler/MIRT.jl) - -## Index - -```@index -``` - -## Functions - -```@autodocs -Modules = [Sinograms] -``` +[`Sinograms.jl`](https://github.com/JeffFessler/Sinograms.jl). +is a Julia language package +with tools for tomographic image reconstruction, +including +* describing ray geometries (parallel beam, fan beam, etc.) +* performing filtered back-projection (FBP) image reconstruction +* system models for iterative reconstruction + +It is all a WIP. +See the examples + +See also: +* Michigan Image Reconstruction Toolbox + [(MIRT)](https://github.com/JeffFessler/MIRT.jl) +* [JuliaImageRecon](https://github.com/JuliaImageRecon) diff --git a/docs/src/methods.md b/docs/src/methods.md new file mode 100644 index 0000000..2bdd164 --- /dev/null +++ b/docs/src/methods.md @@ -0,0 +1,10 @@ +## Methods list + +```@index +``` + +## Methods usage + +```@autodocs +Modules = [Sinograms] +``` diff --git a/docs/lit/fbp2_example.jl b/wip/fbp2_example.jl similarity index 100% rename from docs/lit/fbp2_example.jl rename to wip/fbp2_example.jl diff --git a/wip/todo-projector.txt b/wip/todo-projector.txt new file mode 100644 index 0000000..d5bbdc0 --- /dev/null +++ b/wip/todo-projector.txt @@ -0,0 +1,452 @@ +#--------------------------------------------------------- +# # [Tomography](@id 01-tomography) +#--------------------------------------------------------- + +#= +This page gives an overview of the Julia package +[`Sinograms.jl`](https://github.com/JeffFessler/Sinograms.jl). + +This page was generated from a single Julia file: +[01-tomography.jl](@__REPO_ROOT_URL__/01-tomography.jl). +=# + +#md # In any such Julia documentation, +#md # you can access the source code +#md # using the "Edit on GitHub" link in the top right. + +#md # The corresponding notebook can be viewed in +#md # [nbviewer](http://nbviewer.jupyter.org/) here: +#md # [`01-tomography.ipynb`](@__NBVIEWER_ROOT_URL__/01-tomography.ipynb), +#md # and opened in [binder](https://mybinder.org/) here: +#md # [`01-tomography.ipynb`](@__BINDER_ROOT_URL__/01-tomography.ipynb). + + +# ### Setup + +# Packages needed here. + +using Sinograms #: todo +using MIRTjim: jim, prompt +using InteractiveUtils: versioninfo + + +# The following line is helpful when running this example.jl file as a script; +# this way it will prompt user to hit a key after each figure is displayed. + +isinteractive() ? jim(:prompt, true) : prompt(:draw); + + +#= +### Tomography + +[Tomography](https://en.wikipedia.org/wiki/Tomography) +is the process of imaging cross sections of an object +without actually slicing the object. +There are many forms of tomography; +the description here +focuses on +[X-ray computed tomography (CT scans)](https://en.wikipedia.org/wiki/CT_scan). +(See +[SPECTrecon.jl](https://github.com/JeffFessler/SPECTrecon.jl) +for a Julia package related to +[SPECT](https://en.wikipedia.org/wiki/Single-photon_emission_computed_tomography) +imaging.) + +In an +[X-ray CT imaging system](http://en.wikipedia.org/wiki/File:Ct-internals.jpg), +X-rays emitted from an X-ray source +are transmitted through an object +(e.g., a patient in medical CT) +towards a detector array. +The source and detector +[rotate rapidly]("https://www.youtube.com/watch?v=2CWpZKuy-NE) +around the object. +The signal intensities recorded by the detector +are related +to the sizes and densities +of the object materials +between the source and each detector element. + +### Radon transform + +The mathematical foundation +for 2D X-ray CT imaging +is the +[Radon transform](https://en.wikipedia.org/wiki/Radon_transform), +an +[Integral transform](https://en.wikipedia.org/wiki/Integral_transform) +that maps a 2D function +``f(x,y)`` +into the collection of line integrals +through that function. +Here we describe each line +by its angle ``\phi`` +measured counter-clockwise from the ``y`` axis, +and by the radial distance ``r`` of the line from the origin. +The collection of integrals +is called a sinogram. + +Mathematically, +the Radon transform ``p(r,\phi)`` of ``f(x,y)`` is defined by +```math +p(r,\phi) = \int_{-\infty}^{\infty} +f(r \cos \phi - l \sin \phi, +r \sin \phi + l \cos \phi) \, \mathrm{d} l. +``` + +The Radon transform of the object that is 1 inside the unit disk +and 0 elsewhere is given by +```math +p_1(r,\phi) = 2 \sqrt{1 - r^2} \mathbb{1}_{\abs{r} < 1}, +``` +where ``\mathbb{1}`` denotes the +[indicator function](https://en.wikipedia.org/wiki/Indicator_function), +i.e., +it is 1 for ``\abs{r} < 1`` and 0 elsewhere. + +By the Radon transform's translation property, +the Radon transform +of a disk of radius ``r_0`` +centered at ``(x_0,y_0)`` +is given by +```math +p(r,\phi) = r_0 p_1 +\left( \frac{r - (x_0 \cos \phi + y_0 \sin \phi)}{r_0}, ϕ \right) +``` + +### Sinogram example + +Here is an display of that Radon transform +for a disk of radius ``3`` +centered at coordinate ``(5,1)``. +Note that maximum value is approximately ``6``, +the length of the longest chord +through a disk of radius ``3``. +=# + +nr = 128 +dr = 20 / nr +r = ((-(nr-1)/2):((nr-1)/2)) * dr # radial sample locations +na = 130 +ϕ = (0:(na-1))/na * π # angular samples +proj1 = r -> abs(r) < 1 ? 2 * sqrt(1 - r^2) : 0. +proj2 = (r, ϕ, x, y, r0) -> r0 * proj1(r/r0 - (x * cos(ϕ) + y * sin(ϕ))/r0) +sino = proj2.(r, ϕ', 5, 1, 3) +jimsino = (sino, title) -> jim( + r, ϕ, sino; title, aspect_ratio=:none, + xlabel = "r", ylabel = "ϕ", ylims=(0,π), yticks=([0, π], ["0", "π"]), + yflip=false, xticks = [-10, 0, 2, 5, 8, 10], + clim = (0, 6), +) +jimsino(sino, "Sinogram for one disk") + + +#= +As the preceding figure illustrates, +the Radon transform of a unit disk +has a somewhat sinusoidal shape. +Indeed every point in the ``(x,y)`` plane +traces out a distinct sinusoid +in the sinogram. + +The mapping from the object ``f(x,y)`` +to data like the sinogram ``p(r,\phi)`` +is called the "forward problem." + + +### Image reconstruction + +[Image reconstruction](http://en.wikipedia.org/wiki/Image_reconstruction) +is the process of solving the +[inverse problem](https://en.wikipedia.org/wiki/Inverse_problem) +of recovering the object ``f(x,y)`` +from a measured sinogram, +i.e., +(usually noisy) samples of ``p(r,\phi)``. + + +### FBP + +A simple image reconstruction method +is called the +"filtered back-projection" (FBP) approach. +It works by filtering each row of the sinogram +with a filter, +called the +[ramp filter](https://en.wikipedia.org/wiki/Radon_transform#Radon_inversion_formula), +whose frequency response +is roughly ``\abs{\nu}``, +where ``\nu`` is the spatial frequency variable +(units cycles / m), +followed by a +[back-projection](https://en.wikipedia.org/wiki/Radon_transform#Dual_transform) +step that is the +[adjoint](https://en.wikipedia.org/wiki/Hermitian_adjoint) +of the Radon transform. + + +### FBP example + +Here is an illustration +of using the `fbp` method +in this package +to perform image reconstruction +from the preceding sinogram. + +=# + +image = fbp(sino; dr) +(nx,ny) = size(image) +dx = dr # default +x = ((-(nx-1)/2):((nx-1)/2)) * dr +y = x +jim(x, y, image, "FBP image", + xtick=[-10, 0, 2, 5, 8, 10], + ytick=[-10, 0, -2, 1, 4, 10], +) + + +#= + +The FBP reconstructed image +looks pretty similar +to a disk of radius ``3`` +centered at ``(5,1)`` +as expected. +However, +there are quite a few ripples; +these are +[aliasing](https://en.wikipedia.org/wiki/Aliasing) +artifacts +due to the finite sampling +in ``r`` and ``\phi``. + + +### Noise effects + +Simulating the effects of measurement noise in sinogram +leads to even worse FBP results. + +=# + +noisy_sinogram = sino + 0.1 * randn(size(sino)) +jimsino(noisy_sinogram, "Noisy sinogram") + +# +noisy_fbp_image = fbp(noisy_sinogram; dr) +rmax = maximum(r) +fovmask = @. sqrt(abs2(x) + abs2(y)') ≤ rmax +jim(x, y, fovmask) +noisy_fbp_image .*= fovmask +jim(x, y, noisy_fbp_image, "Noisy FBP image"; clim=(0,1)) +throw(0) + + +#= + +Many computational imaging methods use system models +that are too large to store explicitly +as dense matrices, +but nevertheless +are represented mathematically +by a linear mapping `A`. + +Often that linear map is thought of as a matrix, +but in imaging problems +it often is more convenient +to think of it as a more general linear operator. + +The `LinearMapsAA` package +can represent both "matrix" versions +and "operator" versions +of linear mappings. +This page illustrates both versions +in the context of single-image +[super-resolution](https://en.wikipedia.org/wiki/Super-resolution_imaging) +imaging, +where the operator `A` maps a `M × N` image +into a coarser sampled image of size `M÷2 × N÷2`. + +Here the operator `A` is akin to down-sampling, +except, rather than simple decimation, +each coarse-resolution pixel +is the average of a 2 × 2 block of pixels in the fine-resolution image. +=# + +# ### System operator (linear mapping) for down-sampling + +# Here is the "forward" function needed to model 2× down-sampling: + +down1 = (x) -> (x[1:2:end,:] + x[2:2:end,:])/2 # 1D down-sampling by 2× +down2 = (x) -> down1(down1(x)')'; # 2D down-sampling by factor of 2× + +# The `down2` function is a (bounded) linear operator +# and here is its adjoint: +down2_adj(y::AbstractMatrix{<:Number}) = kron(y, fill(0.25, (2,2))); + + +#= +Mathematically, and adjoint is a generalization +of the (Hermitian) transpose of a matrix. +For a (bounded) linear mapping `A` between +inner product space X +with inner product <.,.>_X +and inner product space Y +with inner product <.,.>_Y, +the adjoint of `A`, denoted `A'`, +is the unique bound linear mapping +that satisfies +_Y = _X +for all x ∈ X and y ∈ Y. +One can verify that the `down2_adj` function +satisfies that equality +for the usual inner product +on the space of `M × N` images. +=# + + +# ### LinearMap as an operator for super-resolution + +#= +We now pick a specific image size +and define the linear mapping +using the two functions above: +=# + +nx, ny = 200, 256 +A = LinearMapAA(down2, down2_adj, ((nx÷2)*(ny÷2), nx*ny); + idim = (nx,ny), odim = (nx,ny) .÷ 2) + +#= +The `idim` argument specifies +that the input image is of size `nx × ny` +and +the `odim` argument specifies +that the output image is of size `nx÷2 × ny÷2`. +This means that when we invoke +`y = A * x` +the input `x` must be a 2D array of size `nx × ny` +(not a 1D vector!) +and the output `y` will have size `nx÷2 × ny÷2`. +This behavior is a generalization +of what one might expect +from a conventional matrix-vector expression, +but is quite appropriate and convenient +for imaging problems. + +Here is an illustration. +We start with a 2D test image. +=# + +image = shepp_logan(ny, SheppLoganToft())[(ny-nx)÷2 .+ (1:nx),:] +jim(image, "SheppLogan") + + +# Apply the operator `A` to this image to down-sample it: + +down = A * image +jim(down, title="Down-sampled image") + + +# Apply the adjoint of `A` to that result to "up-sample" it: +up = A' * down +jim(up, title="Adjoint: A' * y") + + +# That up-sampled image does not have the same range of values +# as the original image because `A'` is an adjoint, not an inverse! + + +# ### AbstractMatrix version + +#= +Some users may prefer that the operator `A` behave more like a matrix. +We can implement approach from the same ingredients +by using `vec` and `reshape` judiciously. +The code is less elegant, +but similarly efficient +because `vec` and `reshape` are non-allocating operations. +=# + +B = LinearMapAA( + x -> vec(down2(reshape(x,nx,ny))), + y -> vec(down2_adj(reshape(y,Int(nx/2),Int(ny/2)))), + ((nx÷2)*(ny÷2), nx*ny), + ) + +#= +To apply this version to our `image` +we must first vectorize it +because the expected input is a vector in this case. +And then we have to reshape the vector output +as a 2D array to look at it. +(This is why the operator version with `idim` and `odim` is preferable.) +=# + +y = B * vec(image) # This would fail here without the `vec`! +jim(reshape(y, nx÷2, ny÷2)) # Annoying reshape needed here! + + +#= +Even though we write `y = A * x` above, +one must remember that internally `A` is not stored as a dense matrix. +It is simply a special variable type +that stores the forward function `down2` and the adjoint function `down2_adj`, +along with a few other values like `nx,ny`, +and applies those functions when needed. +Nevertheless, +we can examine elements of `A` and `B` +just like one would with any matrix, +at least for small enough examples to fit in memory. +=# + + +# Examine `A` and `A'` + +nx, ny = 8,6 +idim = (nx,ny) +odim = (nx,ny) .÷ 2 +A = LinearMapAA(down2, down2_adj, ((nx÷2)*(ny÷2), nx*ny); idim, odim) + +# Here is `A` shown as a Matrix: +jim(Matrix(A)', "A") + +# Here is `A'` shown as a Matrix: +jim(Matrix(A')', "A'") + + +#= +When defining the adjoint function of a linear mapping, +it is very important to verify +that it is correct (truly the adjoint). + +For a small problem we simply use the following test: +=# +@assert Matrix(A)' == Matrix(A') + +# For some applications we must check approximate equality like this: +@assert Matrix(A)' ≈ Matrix(A') + + +# Here is a statistical test that is suitable for large operators. +# Often one would repeat this test several times: +T = eltype(A) +x = randn(T, idim) +y = randn(T, odim) + +@assert sum((A*x) .* y) ≈ sum(x .* (A'*y)) # = + + +# ### Reproducibility + +# This page was generated with the following version of Julia: + +io = IOBuffer() +versioninfo(io) +split(String(take!(io)), '\n') + + +# And with the following package versions + +import Pkg; Pkg.status() +#using ImagePhantoms: shepp_logan, SheppLoganToft